Error Analysis of Some Demand Simplifications in Hydraulic Models of Water Supply Networks

Mathematical modeling of water distribution networks makes use of simplifications aimed to optimize the development and use of the mathematical models involved. Simplified models are used systematically by water utilities, frequently with no awareness of the implications of the assumptions used. Some simplifications are derived from the various levels of granularity at which a network can be considered. This is the case of some demand simplifications, specifically, when consumptions associated with a line are equally allocated to the ends of the line. In this paper, we present examples of situations where this kind of simplification produces models that are very unrealistic. We also identify the main variables responsible for the errors. By performing some error analysis, we assess to what extent such a simplification is valid. Using this information, guidelines are provided that enable the user to establish if a given simplification is acceptable or, on the contrary, supplies information that differs substantially from reality. We also develop easy to implement formulae that enable the allocation of inner line demand to the line ends with minimal error; finally, we assess the errors associated with the simplification and locate the points of a line where maximum discrepancies occur.


Introduction
In the task of mathematical modeling of such complex structures as water distribution networks (WDNs) the use of simplifications aimed to optimize the development and use of the models is unavoidable.Such simplifications stem from the complexity of the modeled infrastructure and, at the same time, are related to the large spatial distribution typical of WDNs.These models are applied in all the areas of hydraulics-including urban hydraulics [1,2].Currently, with the generalized use of geographic information systems, models containing even hundreds of thousands of pipes and nodes are being built [3].
Extremely detailed modeling of real WDNs, even under the unrealistic hypothesis in which uncertainty can be ignored, produces a substantial amount of data and requires sophisticated computational tools and mechanisms to reliably interpret the obtained results in terms of what occurs in the system.Current computational power can be used to build hydraulic simulation models capable of providing a very detailed and accurate model but not an improved understanding of the main structure of the underlying system.Also, as such models would require very dedicated calculations, certain aspects must be considered to ensure efficient implementations.Specifically, optimization of WDNs, whether used for planning or operational purposes, often requires many iterations [4], each involving computationally expensive simulations and huge computer memory.
Let us mention just a few of the most typical simplifications: the use of one equivalent pipe to represent two or more parallel or series pipes; the removal of short pipe segments including dead ends, service connections, and hydrants; the distribution of emitter exponents without considering leakage profiles along the WDN; the assumption of friction factor values without a detailed consideration of the pipes state; model calibration assuming values of one or more of the calibration elements as fixed; the use of a single friction factor for the entire WDN or for an entire sector; or the use of a single emitter leakage exponent for the entire network.The last two simplifications, in particular, are inescapable when working with EPANET [5], a tool of water network analysis in general use worldwide.Also, some simplifications can totally undermine studies of hydraulic transient phenomena in WDNs.See, for example, [6,7].
Among the various simplifications used in the analysis of water distribution systems some are assumed in a generalized way without additional introspection.They are simplifications derived from the different levels of granularity at which a network can be considered.As mentioned before, real water distribution networks, especially those of large cities, cannot be efficiently modeled in their entirety.In approximate models some granulation reduction actions are performed such as skeletonization, grouping, pruning, and clustering.See, among others, [8][9][10][11][12][13][14][15][16][17][18].
In particular, one of these simplifications is the grouping of the consumptions associated with a line in one or both ends of the line.These points concentrate all the existing consumption points (users) within the line.The need for implementing this specific simplification is evident because a faithful consideration of demand would imply the inclusion of a large number of nodes equal to the number of consumption points.In a branched network with few consumers, this would not represent a major problem.However, problems arise in WDNs with up to 30 connections per line (e.g., a street pipeline).In a large WDN (e.g., a 500 km WDN) it would amount to considering about 150,000 nodes, which is impractical when it comes to the construction of the network model, the performance of the calculations, and the display and understanding of the results.This simplification copes with the continuity principle or conservation of mass.However, the energy aspects are completely ignored.
The methods for calculating demand load in the models are one of the least studied aspects despite the obvious importance for defining model reliability.Recently, Giustolisi et al. [19] have addressed this problem from a global perspective and have developed a matrix transformation approach that changes the classical solution of the nonlinear system of equations describing a WDN based on the conjugate gradient [20], into what they call enhanced global gradient algorithm (EGGA).Note that the solution in [20] is implemented within the widely available freeware EPANET 2, used in a generalized way by many engineers and practitioners around the world.EGGA reduces the size of the mathematical problem through a transformed topological representation of the original network model that preserves both mass and energy balance equations and improves the numerical stability of the solution procedure.According to the authors, EGGA significantly improves the model's computational efficiency without sacrificing its hydraulic accuracy.
However, this solution, although technically impeccable from a theoretical point of view, exhibits a number of practical drawbacks.Firstly, those matrix transformations are not implemented within EPANET and, as a result, are not available to the huge community of its users.Secondly, the transformations are too complex for most users of this program and, especially, too complex to be incorporated into the EPANET toolkit, given the in-depth knowledge of programming techniques that their implementation would require.Thirdly, the transformation matrices must previously be explicitly written; even though obtaining these matrices is relatively straightforward, the computational efficiency, at least in terms of memory needs, is not evident.One has to use several additional (very sparse) matrices of sizes or the order of the number of resulting pipes times the number of original pipes, and the number of resulting pipes times the number of unknown head nodes, are explicitly used.When solving real-world problems with hundreds of thousands of nodes and pipes, this may become a serious problem.Fourthly, for networks already modeled by excluding those intermediate demand nodes, the solution of matrix transformation will be useless.Fifthly, when planning and designing a new network starting from the household demand distribution, it would be perhaps desirable to start building the model by performing the simplification from the outset, in order to avoid later complications.
Compelled by those drawbacks we have addressed the problem from another, we claim, more practical point of view.This paper analyzes the possible errors from the effect of using certain types of simplifications when loading demands in models, specifically, the widespread 50% rule, which allocates half of the in-line demand to each line end.We analyze to what extent this simplification is acceptable or, on the contrary, supply information that differs substantially from reality.Also, we obtain formulae that enable to allocate inner line demand to the line ends with minimal error.Finally, a calculation of the maximum head point discrepancy associated is provided.
Our proposal involves simple, direct methods that can be easily applied by any user of any WDN analysis package, since emphasis is not placed on programming ability, but on how to make a decision about the technical aspect of simplifying the model and, thus, load the demand properly.Users, having already developed models of their networks, may revise the allocation rule used and replace it, if necessary, with the values provided by the new formulae, what will enable them to obtain more reliable results.Also, users starting the model of a new network may make an a priori decision about how to simplify the network and suitably implement the associated simplifications.
We first present a simple case that enables us to shed light onto the problem: a single line with variable distribution and granulation consumption is considered.Then an example of a real network is analyzed using the lessons learned.The paper closes with a conclusions section.

Line with Associated Demand
Let us first consider the case of a single line associated with some internal consumption under steady state condition.Such a line is representative of the simplest installation (a line between two nodes) with a given inflow rate.The characteristics of the line are (i) length: , (ii) diameter: , (iii) upstream head (boundary condition at the upstream node):  0 , (iv) friction factor: , with its associated line resistance:  = 8/( 2  5 ), Various demand scenarios of consumption in the line may be tested.Such scenarios are associated with two characteristics: (i) total demand in the line with regard to inflow, (ii) specific distribution of the demand along the line.
Let us assume that the flow consumed within the line (total in-line demand) represents a percentage of the line inflow.If this fraction is represented by   , 0 <   ≤ 1, the actual demand in the line is given by the expression (1)

Uniformly Distributed Demand along the Line.
To start with the study, we will assume that the actual demand of the line is uniformly distributed into a number  of equally spaced interior points (nodes);  can take a value ranging from 1 (in the case of a line with a single demand node in the middle) to a large integer number (in the case of an equally distributed demand throughout the line).Observe that we do not consider any demand at the end nodes, since we are only interested in the line demand associated with the interior nodes.
Figure 1 shows various distributions of piezometric head corresponding to values of   equal to 1, 0.8, 0.6, and 0.4, for a set of values of , ,  0 , , and  in .To build Figure 1 we have used the specific values  = 500 m,  = 300 mm,  0 = 50 m,  = 0.018, and  in = 0.25 m 3 /s.As mentioned before, the demand has been equally distributed among  equally spaced interior nodes.Specifically, in Figure 1,  takes the values 1, 3, 7, 11, 19, and "infinity." The "infinity" case represents a uniform continuous demand.The various curves in Figure 1 have straightforward interpretation.
For the polygonal hydraulic grade lines (HGLs) made out of segments between consumption points, the calculations correspond to the usual hydraulic calculation of losses.The polygonals start at the boundary condition (0,  0 ), the other vertices being the  + 1 points as follows: ) ,  = 1, . . .,  + 1. ( Let us call   (  , ; ) these HGLs, the subindex "" standing for "real" distribution of piezometric head along the line, as real demands are used to calculate (2).
The calculation for the ideal HGL corresponding to a uniform continuous demand, which is used here as the limit for  → ∞ of the discrete uniform distribution of demands, is performed by integrating the differential loss along the line.The value   / is the (constant) demand per unit length and  is the distance to the upstream node.By integrating, and using (1), the piezometric head corresponding to this continuous loss is given by which corresponds to the upper curve, a cubic, in each of the graphs in Figure 1.
It becomes clear that the greatest discrepancies occur for values of   close to 1 (e.g., when a high percentage of the inflow is consumed along the line).
As mentioned before, these "real" HGLs in Figure 1 have been calculated according to the demand distribution at the various inner points in the line.However, models of large WDNs do not usually take intermediate demands into account; in contrast, the demand of each line is allocated to the end nodes of the line, the 50% rule being generally used.

Allocation of In-Line Demand to the Line Ends: Is the 50%
Rule Adequate?Let    be the factor that allocates a part of the line distributed demand,   , to its upstream end.Thus, the demand assigned to this upstream node is  0 =      .As a result,   =  in −  0 is the flowrate through the line.Note that, using (1), Then, the calculated head value   for a given value of The HGL obtained is, thus, a straight line that connects the point (0,  0 ) with the point (,   (   ; )).This last value corresponds to the calculated head at , the downstream node.
In Figure 2, dashed lines have been added to the first two charts of Figure 1.These new HGLs have been calculated to give the same piezometric head at the downstream node as the line corresponding to a demand concentrated in the middle point,  = 1, for the case   = 1 (left chart), and as the line corresponding to a continuous demand for   = 0.8 (right chart).These lines have been obtained by allocating a fraction of the interior line demand to the upstream node and the rest to the downstream node.Analogous dashed lines can be obtained for other combinations of   and .If the allocated fractions to the line ends are different, the corresponding (straight) lines (the lines given by the numerical model with lumped demands at the ends of the line) will also be different.
Two problems arise at this point.
(a) Firstly, it would be desirable to know the best allocation of the total in-line demand to the line ends, that is to say, to know the value    that solves the following problem: Minimize We start by discussing (7), and then address (8) in Section 2.4.
The solution of ( 7) is here constrained by the nature of the problem: we have to adhere to the fact that one or more lines (pipes) may be connected to the downstream end of the considered line.The connected lines need the correct piezometric head at -upstream end for themto suitably perform their respective calculations.It means that the piezometric head at , given by   and   , must coincide.That is to say, (7) reduces, in our case, to Solve   (   ; ) =   (  , * ; ) for    .
By solving (9) the following expressions for    are obtained.
(i) Case of continuous demand: (ii) Case of demand equally distributed among  equally distributed nodes: In Table 1 and in two-dimensional Figure 3, values for (10) and (11) for   values between 0.2 and 1 and for  varying along the previously used values, namely, 1, 3, 7, 11, 19, and infinity, are presented.Note that values 1 to 6 on the frontal axis of Figure 3 symbolize, as shown in Table 1, the values  = 1, 3, 7, 11, 19, and infinity, respectively.
The following facts are remarkable.
(i) These values are independent of the problem data, namely,  0 ,  in , , , and , and depend only on   and , in the case of (11); that is to say, they depend on the magnitude and the pattern of the distributed demand.This is a very remarkable result, since the present study thus becomes nondimensional and, as a result, completely general.(ii) The values of    range from approximately 0.3 to 0.5.This means that about 30%-50% of the total demand must be allocated to the upstream node and the remainder to the downstream node.
(iii) The lowest    values correspond to the most awkward cases: the rate of demand is close to or equals the total inflow in the line, and the demand is highly concentrated at a few points (upper right corner of the table, right front of the figure).
(iv) The highest    values, closer to 50%, correspond to the less problematic cases, meaning little total distributed demand in relation to the total inflow in the line, and widely distributed demand (lower left corner of the table, bottom left of the figure).This value approaches 50% as the rate of inflow consumed in the line approaches zero.(Observe that for both (10) and (11) lim   → 0  (⋅)   (  ) = 0.5.)(v) It is also worth noting that lim  → ∞     =  ∞   for all   ; that is to say, the monotonic sequence of continuous functions {    } ∞ =1 converges to the continuous function  ∞   uniformly in [0, 1] (Dini theorem; see, e.g., [21,22]).This has a direct interpretation as the continuous demand is the limit of a uniformly distributed demand among an increasing number of equally distributed points on the pipe.
As mentioned before, it is common practice in mathematical modeling of WSN engineering to distribute the line flow into two parts: 50% for the upstream node and 50% for the downstream node, which approximately coincides with what is observed in Table 1 and Figure 3, except for cases where the demand line is highly concentrated and represents a large percentage of the total flow through the line.
As a result of what has been presented so far, it can be said that, for uniformly distributed in-line demand, the usual 50% rule seems, in principle, an appropriate solution provided that the inner demand of the line is small compared with the pipe inflow and that such a demand is widely distributed.However, equal demand allocation to the end nodes of the line may produce important discrepancies, since the study highlights the need for other assignments in certain cases.

Arbitrary Demand along the Line.
To state the problem in its more general form, let us now consider a demand distribution on the line whose accumulated demand is given by a function () =   (), where () is the accumulated demand ratio, a function increasing monotonically from 0 to 1.While   is calculated as in ( 6),   is calculated by integrating the loss Δ = −( in − ()) 2 Δ through the line [0, ] as follows: Observe that this function is monotonically decreasing and concave upwards.
Using (12) and the expression (6) for   in , written, using (5), as the equation   =   in  may be rewritten as By substituting one gets from where the following expression is readily obtained: The general solution for the problem at hand, when an arbitrary demand through the line is considered, may only be solved after having a specific expression for ().As a consequence, for arbitrary demands we will restrict ourselves to the case of discrete demands on a finite number of points of the pipe, as happens in real life.which concentrates the whole demand   at  1 , where (⋅) is the well-known Dirac delta.In this case, ( 12) is written as where  =  1 / is the fraction of the pipe where the concentrated demand is located from the origin.Equating again   =   at  gives As could be expected, these values not only depend on   , as in the case of uniform demand, but also strongly depend on .In Figure 4 we have plotted these values as a function of  for various instances of   , namely, 1.0, 0.8, 0.6, 0.4, and 0.2.
As expected, the worst cases occur once more for withdrawals representing a large percentage of the inflow to the line.For small in-line demands (see, e.g., the curve for   = 0.2; also calculate the limit of (20) for   approaching to 0), the demand should be allocated to the end nodes almost linearly proportional to the relative distance of the demand point to the downstream end, as is completely natural.In contrast, this rule does not apply to large in-line demands, as their corresponding curves show, by becoming less and less linear.As an extreme case, let us consider the case of   = 1.Various HGLs have been plotted in Figure 5 by varying the location of the withdrawal point, specifically for values  = 0.1, 0.25, 0.5, and 0.9.The (straight dashed) line, corresponding to the 50% allocation rule, has also been represented.We can observe that only for  = 0.25 does the dashed line match the correct piezometric head at the downstream end (also, observe that the lower curve in Figure 4, corresponding to   = 1, contains the point (,    ) = (0.25, 0.5)).In the other cases, disagreements are not only important at the downstream end but all along the line.
In the general case,   is calculated by By denoting this expression can be written as Then, equating again   =   at  gives ) . ( This expression can be easily calculated using, for example, a worksheet, as in Figure 6. In this figure, we consider a demand distribution at the (inner) points given by the values of  (0.246, 0.338, . . ., 0.954), the demand values given by the values (0.171, 0.084, . . ., 0.329) of , representing demand fractions at those points, according to (22).In the worksheet we can also read, besides the specific variable values used for the calculations, the value of   = 0.8 used, meaning that a demand of   = 0.8 ⋅  in = 0.8 ⋅ 0.25 = 0.2 m 3 /s is extracted in the line.By using formula (24), implemented in the cell below    , we obtain the rate of demand that must be allocated to the upstream end to get the correct piezometric head at the downstream end.The graph in Figure 6 represents this situation.It can be clearly observed that the use of (24) provides a calculated HGL for the considered example (mid line) that perfectly matches the right end of the line representing the real HGL (lower polygonal).On the other hand, the application of the 50% rule (upper HGL) produces unacceptable errors.

Maximum Head Discrepancy When
Using the Proposed Formula.As mentioned before, allocation of in-line demands to the end nodes of a line may be of great interest in order to reduce the size of the mathematical model of a WDN.In the previous section, we have given formulae to obtain allocation values that zero the piezometric head error at the downstream end of the line, a mandatory condition for correct calculation on the line(s) connected to this end.However, this reduction of the model size is at the price of making some piezometric head errors at the inner points of the line.The engineer analyzing a WDN should be aware of the magnitude of those discrepancies, in order to have a better representation and understanding of the problem.In this section we answer this question by solving problem (8).
25 500 0.018 0.3 0.61205079 100 In the (ideal) case of a continuous demand along the line, it is easy to show that (8) reduces to find the point  where the derivatives with respect to  of Let us observe that   (  ; ) is a differentiable function since () is continuous and that   (   ; ) is a linear function of .Thus, it is easily seen that the maximum discrepancy is located at the point  given by Example 4. In the case of uniform continuous demand the solution of ( 8) is given by obtained using () = /, as in Example 1, and (10).
In the (real) case of a discrete demand along the line, the maximum discrepancy occurs at one of the points   , since the real HGL is a decreasing concave upwards polygonal (see Figure 6).Then, the problem reduces to identify the first  0 for which the next section of the polygonal has a slope equal to or lower than the slope of   (   ; ) in absolute values (if equal, all the points between  0 and  0+1 will provide the maximum since the mentioned section of the polygonal and   (   ; ) will run parallel between both points).So we have to solve the following problem.
Find the first point  0 such that This inequality may be simplified and written, in terms of   , as which gives in terms of   =   /  .That is to say, the problem may be rewritten as follows: Find the first point  0 such that Example 5.In the case of uniform demand on  equally spaced points, this gives the following: Find the first point  0 such that Observe that lim  → ∞ (  0 /) ≥ lim  → ∞     =  ∞   , since     converges to  ∞   .We also have that lim  → ∞ (  0−1 /) ≤ lim  → ∞     =  ∞   .But obviously these two limits are the same.As a result, one has from (28) that lim Observe that   0 is the point where the polygonals (2) exhibit the biggest head discrepancy with the calculated head given by ( 9) corresponding to the allocation factor     .Also,  ∞ is the point of maximum head discrepancy between the HGL corresponding to the uniform continuous demand and the calculated head corresponding to the allocation factor  ∞   .As we could expect, lim Expression (32) does not need any additional extra work in the general case if the calculations have been organized as shown in Figure 6.In effect, in this figure the first value of the accumulated rated demand that exceeds    = 0.261 is ∑ 0 =1   = 0.272, which corresponds to the third inner demand point  0 = , with  = 0.604.This maximum is   (0.261;  0 ) −   (0.8;  0 ) = 92.767− 90.763 = 2.003. (36)

Illustrative Example on a Real-World Network
This section is aimed to present some significant results obtained when applying the lessons learned in the previous section to a practical situation.The assessment is carried out in a real water supply network from Tegucigalpa City (the capital of Honduras).The network, which is intended to be a district metered area, is supplied by one main pipe connected to one of the water tanks administered by the local water authority (SANAA).The network has 203 nodes and 211 pipes.Figure 7 shows the network model in EPANET.To illustrate the aspects previously studied, three scenarios have been put forward.For each scenario, intermediate demand nodes (mid nodes) are placed in the main pipe.In each case, the demand values loaded in the intermediates nodes, the distribution along the line, or the pipe diameter vary as it is indicated as follows.
Case A. As it may be seen in Figure 8, an internal demand representing 60% of the total inflow in the main line is fixed in two nodes placed along the line.The first one, loaded with 25% of the total internal demand, is located at 50% of the total length.The second intermediate node is located at 80% of the total length and the remaining 75% of the internal demand is allocated to it.
Case B. As shown in Figure 9, one single node is placed along the main line, at 50% of its total length.As it can be expected, 100% of the total internal demand is allocated to it.In this case the internal demand represents 80% of the total inflow.
Case C.This case is fairly similar to Case B. However, in this case, the diameter of the main line is reduced from 600 mm to 450 mm.

Result (Case A)
. Figure 8 clearly reflects the characteristics and the results obtained for Case A. First, a hydraulic simulation was conducted with the internal demand loaded in the mid nodes.For that (real) case, the pressure in node B (downstream node) is 42.35 water column meters (wcm).When the internal demand is reallocated to the upstream and downstream nodes following the 50% rule, the result varies up to 12.2 wcm in comparison with the initial case.Through the application of formula (24), a more suited demand allocation rule is found (22.45% for the upstream node and 77.55% for the downstream node).When this rule is applied, the resulting pressure value is, as expected, exactly the same as in the (real) case with the internal demand located in the mid nodes.

Result (Case B and Case C).
Given the fact that between Cases B (Figure 9) and C (Figure 10) there is only one difference, namely, the pipe diameter, the allocation rule in both cases is the same; 34.86% of the internal demand should be allocated to the upstream node and the rest to the downstream node.Nevertheless, the differences between the initial pressure values and the pressure obtained when applying the 50% allocation rule are larger in Case C than in Case B, since in Case C the head loss is higher than that in Case B. If Case A is included in the same comparison (higher head loss than Cases B and C), the importance of a demand allocation based on the proposed methodology over a 50%-50% becomes evident.
Finally, to better visualize the negative impact of using the common 50% demand distribution rule instead of the suited rule (24) developed in this paper, resulting pressure values obtained in four nodes of the example network (see Figure 11) are compared for Case A: first, when the demand is allocated using the 50% rule, and then, when the suited rule, given by formula (24), is implemented.The results of such comparison may be seen in Table 2.

Conclusions
The complexity of the interaction of all the input data in a mathematical model makes it impossible to include all the data accurately.The reality is very complex [2], and the use of simplifications in order to make the model feasible is therefore inevitable.The use of such tools should be, in any case, accompanied by a clear understanding of the consequences of such assumptions.It is obvious that the lower level of simplification corresponds with more complex tools.In this sense, software packages available in the market devoted to the analysis, design, and, in general, the simulation of the various states in a WDN should be used with the necessary caution.
This research focuses on the study of the influence that the concentration of a distributed demand in a line on the line ends represents in modeling steady state conditions in

Example 2 .
Let us start by considering a single demand withdrawn at a specific point of the line.That is to say, let us consider a demand distribution given by  () =  ( −  1 )

Figure 4 :
Figure 4: Demand allocation fraction to the upstream node depending on the location of a single withdrawal in the line.

Figure 5 :
Figure 5: HGLs depending on the location of a single withdrawal in the line and HGL (dashed line) corresponding to the 50% allocation rule.

Figure 6 :
Figure 6: Calculation of    for a nonuniform distribution of demand within a line and comparison among hydraulic grade lines.

Figure 7 : 2 Q
Figure 7: Water supply system used for the practical example.

Figure 8 :
Figure 8: Details and results for Case A.

Figure 9 :
Figure 9: Details and results for Case B.

Figure 10 :Figure 11 :
Figure 10: Details and results for Case C.

Table 1 :
(9)ues of    as a function of   and  when solving(9).
(2)ure 1: Hydraulic grade lines in one single line for various uniform distributions of demand.forcertainfunctionalnorm‖⋅ ‖, where   (   ; ) is calculated by(6), according to the lumped demand allocation to both line ends, and   (  , * ; ) accounts for the real demand distribution, either calculated by(2)or (4) (or any other more general expression corresponding to not uniformly distributed demand, which we address later).The asterisk denotes other parameters, such as  in (2), which may appear in the expression of   .(b) Secondly, after having made a suitable line-end allocation decision, it is of interest to know the actual distribution of piezometric head errors on the line and, specifically, to identify at what point or points the maximum head discrepancy occurs: Maximize         (   ; ) −   (  , * ; )      .