## Abstract

Minimization of groundwater exploitation cost is examined, considering: (a) pumping from a system of wells up to a central water tank, including friction losses along the connecting pipe network and (b) amortization of network construction. Assuming that the wells are located symmetrically around the tank and directly connected to it, we derived analytically the distance between tank and wells, which minimizes the total cost. Then we compared the minimum cost of this well layout, with that of placing one well at the location of the tank and the rest symmetrically around it. Finally, we dropped any assumption on well layout, we considered that wells are connected to the tank using a minimum spanning tree and we optimized well locations and flow rates using genetic algorithms. For up to 8 wells, the resulting minimum cost is comparable to that of the symmetrical cases, even when the optimal well layout is quite different. Moreover, the analytical solution, derived for the symmetrical case, can serve to evaluate solutions achieved by sophisticated optimization techniques.

## HIGHLIGHTS

Minimizing cost of pumping groundwater from a system of wells to a tank plus amortization of pipe network construction.

Pumping cost includes friction losses along pipe network.

Analytical calculation of minimum cost for symmetrical well layouts.

Cost evaluation for symmetrical well layouts around the tank plus a well at the tank.

Minimization of cost without any well layout constraint, using genetic algorithms.

### Graphical Abstract

## LIST OF SYMBOLS

- c
Pumping cost coefficient that depends on pump efficiency and electricity cost

- c΄
Construction cost coefficient, accounting for pipes' amortization cost

- D
_{j} Pipe's j diameter

- e
Pipe roughness

- f
_{j} Friction factor

- hf
Friction losses

Cost factor, accounting for pipe network construction

Cost factor, accounting for groundwater hydraulic head level drawdown

The overall pumping cost

The overall pumping and construction cost

- K
_{trans} Cost factor, accounting for hydraulic head losses at the pipe network

- L
_{j} Length of pipe j

- P
_{coef} Electricity cost for each kWh

- q
_{i} Pumping flow rate of well i

The overall water demand

- Q
_{j} Pumping flow rate of pipe j

- R
Radius of influence of the wells

- r
_{0} Radius of the well

- r
_{k,i} Distance between wells i and k

- s
_{i} Hydraulic head level drawdown at the well i

- T
Aquifer's transmissivity

- t
_{pump} Potential pipe network lifetime

- u
Water velocity

*γ*_{w}Specific weight of water

*ε*Interest rate

*λ*The ratio between the construction cost coefficient and the pumping cost coefficient

the angle of the lines connecting wells i and j with the tank

## INTRODUCTION

Quantitative and qualitative degradation of fresh water and groundwater systems is an undisputed truth at the global scale, mainly due to the human-induced pollution or over-exploitation of water resources. At the same time, water use in the ‘Global South’ is set to increase markedly over the next few decades as a result of population growth and planned increases in irrigation (Vörösmarty *et al.* 2005). Moreover, climate change is drawing more and more attention, changing the way we approach water resources management (Mani *et al.* 2016). As acceleration of climate change is mainly attributed to human actions, the term ‘Anthropocene’ was introduced (Crutzen 2006) that is now widely used, to describe a new geological era, in which ‘humans are a significant force in the Earth System’. In the aforementioned framework, groundwater is often seen as a solution, as it responds much more slowly to meteorological conditions than surface water and, as such, provides a natural buffer against climate variability and can be seen as an opportunity to access secure water (MacDonald & Calow 2009; MacDonald *et al.* 2012). Nevertheless, to make groundwater resources massively accessible, a cost reduction approach is necessary.

Recently, there has been a significant evolution in terms of computational methods in water resources management (AWWA Engineering Modeling Applications Committee 2020), such as machine learning to predict streamflow (Parisouj *et al.* 2020), support vector machines to predict response of karst springs (Goyal *et al.* 2017) and advanced statistics to estimate groundwater contamination (Troldborg *et al.* 2012). Metaheuristics, such as genetic algorithms, are widely used for optimization purposes in difficult water resources management problems (Cheng *et al.* 2000; Ketabchi & Ataie-Ashtiani 2015; Moutsopoulos *et al.* 2017; Christelis *et al.* 2019; Katsifarakis & Kontos 2020; Seyedpour *et al.* 2021). Genetic algorithms are capable of solving optimization problems, with high complexity level. However, their use depends on availability of adequate computer resources, while there is no guarantee that they find the global optimum. Analytical solutions, when available, should be preferred, following the principle of ‘use of the simplest efficient tool’. Among others, they can serve as rules of thumbs for practitioners.

Moreover, it needs to be mentioned that minimizing water extraction cost has a direct impact in minimizing energy consumption (Kontos & Katsifarakis 2017). Reduction of energy consumption has recently become a priority in all water use sectors (Hoekstra 2009; Herath *et al.* 2014; Pérez-Sánchez *et al.* 2017; Chen *et al.* 2019). It is mainly related to the operation of water facilities, but it can be also associated with the footprint of the construction of the necessary infrastructure.

In many cases, pumping cost minimization has been examined together with quality of the water, which is provided to the end users (Arai *et al.* 2013; Mala-Jetmarova *et al.* 2017). In this paper we study cost minimization of development of groundwater resources, taking into account two major cost components: (a) operational (pumping) from a system of wells up to a central water tank, including friction losses along the pipe network that connects the wells with the tank and (b) amortization of the construction of the aforementioned network, which implies optimizing location of the wells. This approach is quite useful in planning new groundwater development schemes when selection of well locations is part of the optimization problem (Park *et al.* 2021).

As in previous works, we consider steady-state confined groundwater flow in infinite, homogeneous and isotropic aquifers. Moreover, we assume that efficiencies of all pumps are constant and equal to each other.

First, we derive the minimum of the cost analytically, based on the following additional assumption: All wells are placed symmetrically on a circle around the water tank and each of them is directly connected to it. Then we compare the cost of the fully symmetric well pattern with that of placing one well at the location of the tank and the rest symmetrically around it. Finally, we drop any assumption regarding well placement, we consider that they are connected to the tank using a minimum spanning tree (Prim 1957) and we optimize the locations of the wells and their flow rates using genetic algorithms.

### Mathematical formulation of the optimization problem

_{dr}and K

_{trans}. K

_{dr}accounts for hydraulic head level drawdown at each well due to pumping. For infinite homogeneous aquifers under steady flow conditions it is given as:

In Equation (1) n is the number of wells, q_{i} (m^{3}/s) and s_{i} (m) are the flow rate and the hydraulic head level drawdown at the well i, T (m^{2}/s) represents the aquifer's transmissivity, R (m) the radius of influence of the system of wells and r_{k,i} (m) the distance between wells k and i (r_{i,i} = r_{0}, namely it is equal to the radius of the well). Finally, c is a cost parameter, expressed in €·s/m^{4}, that depends on pump efficiency and electricity cost. It is given as c = *γ*_{w} · P_{coef} · t_{pump}, where P_{coef} (€/kWh) represents the electricity cost of each kWh, t_{pump} (s) the duration of pumping and *γ*_{w} (N/m^{3}) the specific weight of water.

In Equation (2), L_{j} (m) stands for the length of pipe j and f_{j} for the respective friction factor, which is calculated using Swamee–Jain formulas (Swamee & Jain 1976). D (m) is used to represent the pipes' diameter. Overall, q represents the pumping flow rates of the wells, while Q represents flow rates through the pipes.

In Equation (3), the term c΄(€/m) is used to take into consideration the pipes' amortization cost (€/m). To get this parameter, a regression analysis was conducted to create a relationship between pipes' diameters and pipes' market value. Then we used the interest rates (*ε*) and networks potential lifetime (t) to get a value for pipes' cost, so that comparisons can be made possible between similar costs (referring to a year).

### The analytical solution for symmetrical well placement

_{oper}(K

_{oper}= K

_{dr}+ K

_{trans}) to pump a given total flow rate q

_{tot}from n wells up to a tank through an existing pipe network is minimized, when the following conditions hold:

_{j}(j [1,n]), it is assumed, for convenience and without loss of generality, that q

_{n}depends on all other q

_{i}. Then,

*Δ*s

_{jn}is the drawdown difference between wells j and n, while

*Δ*hf

_{sj,sn}denotes the difference between the sum of head losses due to friction, along the paths of q

_{j}and of q

_{n}, through the pipe network (from the respective well to the tank). The set of flow rates q

_{j}that minimizes the cost can be found by solving a system of n equations and n unknowns with respect to q

_{j}. The n − 1 equations have the form of Equation (5), while the last one expresses the constraint regarding the well flow rates, namely:

Moreover, it has been proved that the feasible solution for an existing pipe network is unique.

In this paper, we further include amortization of the construction cost of the pipe network in the optimization problem. This approach is more appropriate when planning new groundwater resources development schemes. To derive an analytical solution, we assume that the wells will be placed at the vertices of a regular polygon, with the water tank at the centre of the respective circumcircle. Moreover, we assume that all wells are connected directly to the tank, with similar pipes. This well layout is quite reasonable for up to 6 wells, since, for equilateral hexagons the radius of the circumcircle is equal to the hexagon's side (and smaller for regular polygons with fewer sides).

_{in}and Δhf

_{si,sn}are equal to zero. Then Equation (5) holds, and for a given circumcircle radius, this flow rate distribution is optimal, as shown in Equation (7):

To include construction cost in the minimization problem, though, we have to consider the pipe network length, namely the radius of the circumcircle L, as variable. As shown in Figure 1, when the distance between the tank and the wells increases, K_{dr} decreases, whereas K_{trans} and K_{constr} increase. We seek the critical length value that minimizes the overall cost K_{total}. A simple example follows.

Example 1. Two wells are symmetrically placed with respect to the tank. The cost items are calculated for many distances L between the tank and the wells, for the following data set: T = 0.002 m^{2}/s, q_{tot} = 0.1 m^{3}/s, f = 0.0235, D = 0.3 m, R = 2,000 m, c = 5,156.136 €·s/m^{4}, c’ = 2.8415(€/m) and r_{0} = 0.2 m. As shown in Figure 1, the total cost is minimized when L = 305 m.

#### Critical point and absolute minimum

Based on the aforementioned assumptions, the characteristics of all wells and pipes are the same, namely . Moreover, the flow rate through each pipe equals the flow rate of the respective well, namely .

_{total}, its derivative with respect to L should be calculated. K

_{constr}is proportional to L. Moreover, the friction head loss hf can be expressed using a length-independent parameter J, as: hf = J*L. K

_{dr}depends on L, too, since it depends on the distances between wells, and the respective derivative should be calculated. We have, then:

_{i}and angles

*θ*

_{i}. Using this transformation and the cosine rule, the distance r

_{i,j}between any two wells i and j can be expressed as follows:

In Equation (12), u represents the velocity of the water (m/s). Using the previous relationships, one can get L_{critical} = 305 m in the Example 1, which is the actual length that minimizes the exploitation cost.

As the second derivative is always positive (for *n* > 1), the critical length corresponds to a minimum and is unique.

#### Solution limitations

In all calculations involving K_{dr} we assume that pumping at any well affects drawdown at all others. In practical applications, this assumption should be checked, comparing well distances with the radius of influence R. If R is exceeded, any increase in K_{constr} does not result in further K_{dr} decrease. This procedure is illustrated in Figure 2.

In Equation (14) r_{i,j} is the distance between wells i and j, L the distance between the wells and the tank and *φ*_{ij} is the angle of the lines connecting wells i and j with the tank. If R< r_{ij} there is no interference between the wells and the respective term in Equation (1) should be set to 0. Such a check can be easily incorporated in the drawdown calculation procedure. In practice, such occasions may arise when K_{trans} and K_{constr} are too small compared to K_{dr} and the calculation of L_{critical} is of restricted value.

### Placing one well at the location of the tank

Another reasonable planning approach is to place one well at the location of the tank and the rest symmetrically around it, at a distance L, calculated by Equation (11), as for the fully symmetrical case. Compared with that case, K_{constr} will be always smaller, but K_{dr} may be larger. To find the minimum cost we calculate the optimal flowrate distribution, as described by Nagkoulis & Katsifarakis (2020) and briefly outlined in the previous section, in order to find the respective K_{dr}. Comparison for different data sets appear in a following section.

### Cost minimization using genetic algorithms

In this section we use genetic algorithms to minimize K_{total} for given total well flow rate and number of wells. We drop the assumptions on well layouts, made in the previous sections, and we allow free choice of their locations, inside a square area around the tank. The side of the square is equal to 2R. Moreover, the wells are connected to the tank using a minimum spanning tree. Our purpose is to evaluate the quality of the solutions obtained for the symmetrical well settings.

The aim of the genetic algorithm is to minimize K_{total}, namely its objective function is described by Equation (4). If the number of wells is n, each chromosome (solution of the optimization problem) represents a combination of n well flowrates q_{i} and 2n well coordinates (x_{i} and y_{i}). In the evaluation process, the Prim's algorithm is used to connect the wells to other wells or to the water tank, minimizing the overall length of the network and therefore the construction cost. The minimum spanning tree uses the wells and the tank as vertices and the pipes as lines. It needs to be mentioned that this network layout does not necessarily minimize the friction losses (for given locations of the wells and pumping flow rates), as it might result in high flow rates in some pipes. Genetic algorithms reprioritize hydraulic drawdowns, friction losses or construction cost, based on their relative magnitude, and they may find the global cost minimum. Nevertheless, they offer no guarantee that the optimum is reached.

Our genetic algorithm code is rather simple and includes the following genetic operators: (a) tournament selection including elitism, (b) one-point crossover and (c) mutation. The respective parameters are: (a) number of generations: 2,000, (b) population size: 70, (c) crossover probability: 0.6, (d) mutation probability: 0.012

## RESULTS AND COMPARISONS

First we present some results of the genetic algorithm code, to offer some insight on the cost magnitude and its dependence on the number of the wells. The following data have been used:

The radius of influence of the wells is R = 2,000 m, and the wells' diameter r = 0.2 m. The total flow rate, to be distributed to the wells is q_{tot} = 0.1 m^{3}/s. The friction losses coefficient (f) is calculated using Swamee-Jain formula and a pipe roughness e = 0.5·10^{−3}. Considering the pipes' diameter equal to D = 0.3 m, an interest rate of 5% and supposing that the network will operate for 20 years, we calculated the construction cost parameter for pipes c΄ = 2.8€/m. Similarly, the pumping coefficient c = *γ*_{w} · P_{coef} · t_{pump} = 5,156.136 €·s/m^{4}, using P_{coef} = 0.06€/kwh and considering a duration of pumping equal to a year. Finally, the transmissivity varies from 10^{−3} to 10^{−2} (m^{2}/s), and the number of wells from 1 to 8.

Results for one of the examined scenarios are shown in Figure 3. It can be seen that the genetic algorithm leads to lower costs when many wells are used. It can also be seen that the cost has stabilized before reaching the final generation, which means that increasing the number of generations and the population size would not affect the results.

Then we investigated the two ‘symmetrical placement’ solutions. The sums of the respective operational and annual amortization of construction costs have been calculated and compared with the cost found using genetic algorithms, for 80 scenarios (10 transmissivity scenarios and 8 number of wells scenarios). The difference between the cost results is presented in Figure 4(a) and 4(b), respectively, while Tables 1 and 2 summarize the results for transmissivity values T = 10^{−3} m^{2}/s and T = 10^{−2} m^{2}/s, respectively. The investigation process is depicted in Figure 5.

Wells . | Symmetric . | Central Well . | GA . | L critical . |
---|---|---|---|---|

1 | 75,582 | – | 75,592 | 0 |

2 | 43,916 | – | 43,909 | 611.01 |

3 | 34,171 | 34,066 | 33,679 | 608.34 |

4 | 30,391 | 29,957 | 29,149 | 528.99 |

5 | 28,739 | 28,187 | 26,357 | 456.45 |

6 | 28,022 | 27,435 | 24,674 | 398.22 |

7 | 27,770 | 27,185 | 23,610 | 352.00 |

8 | 27,770 | 27,201 | 22,542 | 314.88 |

Wells . | Symmetric . | Central Well . | GA . | L critical . |
---|---|---|---|---|

1 | 75,582 | – | 75,592 | 0 |

2 | 43,916 | – | 43,909 | 611.01 |

3 | 34,171 | 34,066 | 33,679 | 608.34 |

4 | 30,391 | 29,957 | 29,149 | 528.99 |

5 | 28,739 | 28,187 | 26,357 | 456.45 |

6 | 28,022 | 27,435 | 24,674 | 398.22 |

7 | 27,770 | 27,185 | 23,610 | 352.00 |

8 | 27,770 | 27,201 | 22,542 | 314.88 |

Wells . | Symmetric . | Central Well . | GA . | L critical . |
---|---|---|---|---|

1 | 7,558 | – | 7,558 | 0 |

2 | 5,336 | – | 5,334 | 61.10 |

3 | 4,677 | 4,665 | 4,630 | 60.83 |

4 | 4,456 | 4,411 | 4,332 | 52.89 |

5 | 4,386 | 4,328 | 4,174 | 45.64 |

6 | 4,377 | 4,315 | 4,114 | 39.82 |

7 | 4,397 | 4,335 | 3,975 | 35.20 |

8 | 4,430 | 4,371 | 4,130 | 31.48 |

Wells . | Symmetric . | Central Well . | GA . | L critical . |
---|---|---|---|---|

1 | 7,558 | – | 7,558 | 0 |

2 | 5,336 | – | 5,334 | 61.10 |

3 | 4,677 | 4,665 | 4,630 | 60.83 |

4 | 4,456 | 4,411 | 4,332 | 52.89 |

5 | 4,386 | 4,328 | 4,174 | 45.64 |

6 | 4,377 | 4,315 | 4,114 | 39.82 |

7 | 4,397 | 4,335 | 3,975 | 35.20 |

8 | 4,430 | 4,371 | 4,130 | 31.48 |

From Tables 1 and 2, it can be seen that the differences are small. For up to 4 wells, they are smaller than 5%. If we examine the whole sample, using a regression analysis model, we get an extra 2% increase in average cost difference with each new well. Moreover, high T values generally result in lower differences. Comparing the two symmetrical layouts to each other, we can also see that placing one well at the location of the tank leads to smaller cost. In absolute values, the cost difference increases up to 6 wells and then it decreases slowly. The geometrical explanation is that the side of a regular polygon with up to 5 vertices is larger than the radius of the respective circumcircle. However, such an improvement should be considered minor, as it remains less than 1% in most cases.

The difference in cost between the results of the GA code and those of the symmetrical solutions, is also small. Nevertheless, the optimal well layout reached by the GA code, is quite different from the symmetrical one, as shown in Figure 6. To investigate this discrepancy, we have reduced the pipe diameter from D = 0.3 m to D = 0.1 m, using T = T = 10^{−3} and 4 wells; namely, we have increased substantially the magnitude of K_{trans} and decreased slightly the magnitude of K_{constr}, compared to the previous case.

Figure 7 represents the results obtained for 4 wells. In this figure we can see that the optimal layout reached by the GAs is similar to the symmetrical one, when the diameter decreases. Still the cost is slightly smaller. Overall, increasing the influence of K_{trans} results in layouts in which wells are directly connected to the tank, which reduces the path of the water and therefore the friction losses. Cost results are summarized in Table 3.

Diameter . | Symmetric . | Central well . | GA . | L critical . |
---|---|---|---|---|

0.3 | 30,391 | 29,957 | 29,149 | 528.99 |

0.1 | 42,542 | 41,614 | 40,231 | 73.45 |

Diameter . | Symmetric . | Central well . | GA . | L critical . |
---|---|---|---|---|

0.3 | 30,391 | 29,957 | 29,149 | 528.99 |

0.1 | 42,542 | 41,614 | 40,231 | 73.45 |

Based on the above, the following approach is proposed: calculate the distance L between the tank and the wells, using Equation (11), which leads to the cost minimum, for the fully symmetrical well layout. Use this value for the ‘symmetrical’ layout, with one well placed at the location of the tank. These simple calculations lead to a cost-efficient solution.

To achieve even better results, one can resort to genetic algorithms. In this case, the results of the previous approach can serve to check the quality of the new ones.

## CONCLUSIONS

We have studied minimization of the sum two cost components, which are important, when planning optimal development of groundwater resources: (a) operational (pumping) from a system of wells up to a central water tank, including friction losses along the pipe network that connects the wells with the tank and (b) amortization of the construction of the aforementioned network.

We have derived an optimal solution analytically for any number of wells, under the following, rather restrictive assumptions: (a) the wells are placed at the vertices of a regular polygon, with the water tank at the centre of the respective circumcircle and (b) all wells are connected directly to the tank, with similar pipes. Under these assumptions, the optimal solution entails equal well flow rates and the cost minimization problem is reduced to calculation of the optimal distance between the tank and the wells. We have derived the analytical formula given by Equation (11) to calculate this distance.

Then, we have investigated another ‘symmetrical’ well layout, with one well placed next to the tank. We have concluded that it leads to slightly smaller cost.

Finally, we have investigated the more general cost minimization problem, with no restrictions on well layout, using genetic algorithms. The locations and the pumping rates of the wells are considered as variables, and minimum spanning trees are created to achieve the shortest pipe network, connecting the wells to other ones or to the tank. Comparing the results for the optimal symmetrical well layout, with those of the genetic algorithms (for 80 scenarios) we have concluded that for up to 4 wells, the difference in cost is smaller than 5% and that it generally increases with the number of wells (around 2% for each additional well). An interesting remark is that cost differences are small, even if the respective well layouts are quite different.

For practical purposes, the following approach is proposed: Calculate the distance L between the tank and the wells, using the formula for the fully symmetrical well layout. Use this value for the ‘symmetrical’ layout, with one well placed at the location of the tank, to achieve a cost-efficient solution. To achieve even better results, one can resort to genetic algorithms. In this case, the results of the previous approach can serve to check the quality of the new ones.

Finally, the proposed analytical solution, which corresponds to a local optimum of the general problem, can be used in connection with advanced optimization algorithms, to reduce computational cost. Last, but not least, it can be used to evaluate solutions achieved by sophisticated optimization techniques, which, nevertheless, might be trapped to local optima.

## ACKNOWLEDGEMENT

All computational work has taken place at the facilities of the Aristotle University of Thessaloniki, Greece. Its support is gratefully acknowledged

## CONFLICT OF INTEREST

None

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## REFERENCES

**21**, 2968-2976.

**34**, 4113–4131