EPANET modeling of an urban groundwater field

The paper presents an EPANET model of a groundwater well field. The method used in the simulations to model the variation of the hydrodynamic levels in wells as a function of the pumped flow rate is discussed, and a comparison to previous simulations that used fixed hydrodynamic levels in the wells is performed. The case study points to a groundwater well field in Romania. The results show that the new method although requiring a more complex EPANET model, provides a solution that is closer to the actual variation of water levels in wells.


Introduction
Due to climate change and increased development of urban areas, water resources management (especially drinkable water resources) is gaining nowadays more attention from the water companies. Among the most important aspects under scrutiny are maintaining or improving the quality of the resource, improving the reliability of the resource and of the supply system, and last but not least [1], reduction of the energy needed to transport the resource were it is needed.
In order to solve the complex problems of water supply sources and water distribution networks computer modeling of these systems has gained momentum in the last years.
In this paper the focus is on modeling a well field (which is a rather common water supply source) in EPANET (which is the most used software for water distribution networks modeling).
Reports [2][3][4], papers [5][6][7] and books [8,9] are currently addressing issues of optimizing the groundwater extraction systems, including the pumping equipment. One of the problems is that the flow rate extracted from the well has to be smaller or equal to the optimum well flow rate (determined from hydro geological studies). If the extracted flow rate is bigger than the optimum flow rate, a rapid clogging of the well occurs, rendering it unusable. This is an easily controllable problem for a single well, but becomes more complex for a well field where all the extraction pipes and pumps are connected in a water system. Moreover, the operation of a well field has to be realized sometimes at partial loads either due to low water consumption, or due to revisions or damaged equipment replacement, meaning that some of the pumps are off in the system. In this case, at least intuitively, it is obvious that all duty points of the pumps in operation migrate to the right, at higher flow rates that can exceed the optimal flow rate of the wells. In order to accurately reproduce such a situation, the computer model must allow the level of the water in the well to adjust as a function of the extracted flow rate, to mimic as closely as possible the natural conditions of the aquifer. A methodology for setting up such a model in EPANET will be derived in this paper.
The methodology will be applied to the Mandruloc groundwater well field, the smallest (13 wells) of the well fields supplying the city of Arad, a medium sized city with about 160000 inhabitants, located in western Romania. A comparison to previous simulations performed in EPANET that used fixed hydrodynamic levels in the wells [10] will be performed.

Methodology
As mentioned above, EPANET is a specialized freeware for water distribution networks modeling. It operates with predefined objects that can be used in the simulations, such as nodes, reservoirs, tanks, pumps, pipes and valves. Most of the other objects that might be present in a water distribution network can be simulated:  by some property of a different object − e.g. one can set a pipe property as "check valve" (which means that the flow on the pipe is restricted to a single direction, from the initial node to the end node), or one can set a node property as "emitter coefficient" (which can be used to mimic the existence of all sorts of sprinklers, nozzles, or leakages);  by some property found in the valves category − e.g. pressure reducing valves, pressure sustaining valves, or "general purpose valves" (there are no predefined objects like turbines or wells, but in the EPANET manual [11] there is a mention that such objects can be simulated by the "general purpose valve", where a curve representing the variation of the head loss with respect to the flow can be specified).
The variation of level in an unconfined aquifer due to flow extraction at a well is presented in figure 1, where H represents the hydrostatic level, h represents the hydrodynamic level, r is the radius of the well, R is the radius of the cone of depression (drawdown cone) and Q is the flow rate extracted from the well. The equation that relates the parameters presented in figure 1 is the well known Dupuit formula [12][13][14]: where K (in m/s) is the hydraulic conductivity of the aquifer [15]. This equation can be rewritten to highlight the drawdown (level difference) H = (H − h), as: Of course, equation (2) represents the variation ) and not the other way around, but this fact is not a major concern as in EPANET the function for the General Purpose Valve (GPV) can be added as discrete pairs of values } , . Equation (2) could be easily used in spreadsheet software to yield such pairs of values, provided that all other elements in the formula are known.
Normally, the radius of the well is known. The hydrostatic level and the hydraulic conductivity are also known from hydrogeological studies, so the only problem is to determine the radius of the drawdown cone (radius of influence). There are many empirical formulas for this [1]. In this study performed on the Mandruloc well field, the Kusakin equation [9] was used: A new equation that relates the extracted flow rate with the drawdown H is obtained by inserting the equation (3) into equation (2) Just as equation (2), equation (4) In the case of the Mandruloc study, the only available data were the hydrostatic level of the aquifer, the radius of the wells, and for each well, a pair of values representing the optimal extraction flow rate Qo and the corresponding hydrodynamic level ho. In this case, equation (4) cannot be used directly (due to the lack of information about the hydraulic conductivity of the aquifer), but is valid for the optimal pair of values. Moreover, by substituting the optimal pair of values into equation (4), an implicit formula for K can be obtained: This equation can be solved iteratively, starting with an initial value for K in the right hand side of (5), and computing a resulting value. In the next iteration, the resulting value from the previous iteration becomes the initial value and a new resulting value is computed, and so on and so forth, until the relative difference between the initial value and the resulting one is below a certain threshold. This procedure was used in the case of the Mandruloc well field for all 13 wells. The threshold was set at 10 -6 and the values of the hydraulic conductivity of the aquifer were computed. The resulting values ranged between 10 -5 m/s and 410 -5 m/s (they correspond to loess). Less than 10 iterations were needed to obtain values below the threshold for all wells. Equation (4) was then used to obtain 11 pairs of } , for each well, which were in the sequel inserted as head loss curves in the EPANET model.
A similar methodology can be derived for confined aquifers, if the depth of the aquifer is known.

EPANET model
The total design flow rate of the Mandruloc groundwater well field is about 226 l/s. Water is pumped through submersible centrifugal pumps; each vertical discharge pipe has a diameter of 150 mm and is 30 m long. The vertical discharge pipes are interconnected through pipes with diameters ranging from 200 mm to 600 mm, with a length of 235 m each. The water is supplied to a discharge reservoir through a 6.7 km long pipe with the diameter of 600 mm [10]. The EPANET model of the well field is presented in figure 2.
The EPANET numerical model matches the above geometrical description. The 13 wells of the field are numbered in figure 2 from right to left (the rightmost well is W1, while the leftmost well is W13). The wells are equipped with two types of submersible pumps, namely: Wilo-K 84-2 (wells W4, W8÷W11, W13) and Wilo-EMU K 85-2 (wells W1÷W3, W5÷W7, W12). The elevations of the nodes match the existing configuration, except for the tanks representing wells, where the initial elevation matches the predicted hydrodynamic level of water in the well working at optimal flow rate. The level of the reservoir R1 was set to 112 m, which is the hydrostatic level of the aquifer. The level of the discharge reservoir R2 is 135 m. The roughness on all pipes was set to 2 mm. Check valves were added on the discharge pipes of the pumps. The diameters of the GPV's were set to 100 mm, but this is of no consequence on the model, as long as the head losses between the reservoir and the tanks are computed with respect to the curves calculated according to the methodology presented in section 2. The Darcy-Weisbach formula was used to compute head losses on all pipes.
In order to enable tank levels to adjust in function of the extracted flow rate, an extended time simulation was performed. The total duration of the simulation was 24 hours, with a hydraulic calculation time step and report time step of 2 minutes. The quality calculation time step (that also rules the calculation of the water levels in tanks) was set to 10 seconds.

Results and discussions
In order to be able to compare the results with previous calculations obtained on a model where the water level was considered constant in the well and equal to the optimal level obtained for the optimal extracted flow [10], the previous model (fixed level) was modified here to allow an extended time simulation, with the same parameters as the new model (variable level). Results for 10 a.m. are presented in figure 3(a) for the simulations on the fixed level model, and in 3(b) for the simulations on the variable level model. In figure 3, the pipes are colored in function of the flow rate, while the nodes are colored upon the total head. It might appear from figure 3 that the flow extracted from the wells is more uniform in the variable level model (all the discharge pipes of the pumps transport a flow rate between 15 and 20 l/s). Moreover, on the GPV's supplying the tanks, the flow is different than the one extracted by the pumps and this might seem like a violation of the continuity law. These impressions are wrong. Figure 3(a) is valid for the whole simulation period (as the levels are fixed, nothing changes in the network), but this is not the case for figure 3(b), which is valid only for 10 a.m. (the differences existing between the flow rates on the GPV's and the flow rates of the pumps are constantly changing the levels in the tanks). In the first 3 hours of the simulation, the variations seem random, while after that, a certain repetitive pattern appears. The 10 a.m. moment was chosen at the positive pick of this periodic pattern, after the pattern was established.
The variation of the flow rate through the pipe supplying the discharge reservoir R2, resulted from the variable level model, is presented in figure 4. Of course, variation of level in the tanks was to be expected for the whole extended time simulation (this was the purpose of the variable level model), which would lead to variations of the flow rate trough the discharge pipe. A periodic pattern though seems to be more of the numerical nature. This could be due to the fact that EPANET adjusts the level in the tanks according to the specified quality time step and recalculates the flow rates according to the hydraulic calculation time step. In other words, with options specified in the previous section, levels in the tanks are calculated once every 10 seconds, while flow rates are calculated and reported only once in 2 minutes. This may result in an inability of the software to reach a stable solution. Tests were performed for each well operating singularly in the network (all other pumps closed) and the situation repeated itself.
Globally, as long as mean values on 24 hours are concerned, the results are still comparable with the ones obtained with fixed levels in the tanks. The total flow rate towards the discharge reservoir is of 224 l/s in the fixed level model and the average flow rate (over 24 h) on the same pipe for the variable level model is of 223 l/s. The values of the averaged flow rate extracted from the wells, and of the averaged water level in the well, resulted from the variable level model are presented in table 1, compared to the values of the flow rate and of the hydrodynamic level from the fixed level model.  The variation of the level in well 8 is presented in figure 5, while the variation of the flow rate through the corresponding pump is presented in figure 6. The biggest difference in extracted flow rate is noticed for well 13, although the difference in level is not so important.

Conclusions
In this paper, a methodology for modeling wells in an unconfined aquifer in EPANET, so that the levels of the water in the wells can adjust with respect to the extracted flow rate, was derived. The methodology was then applied to the model of a well field in Romania, and compared to previous calculation performed on a fixed level model [10]. Global results were close to the previous calculations, although differences were identified for particular wells in the field. By adopting this methodology, the EPANET model gains in complexity, but its behavior is more realistic. Further work is needed in order to assess the part of this more realistic behavior that is due to the numerical algorithm of the extended period simulations implemented in EPANET.