Optimal location of pumping wells by a mesh-free numerical method

Water is the critical factor in planning for the future. The issue of water scarcity for countries such as Iran, which has an arid climate, has long been discussed, and access to water resources for drinking purposes, agricultural, and industrial uses is of particular importance. This importance is even greater in areas of Iran that are heavily dependent on groundwater. Groundwater resources in different parts of the world are used in different sectors such as agriculture, environment, irrigation, etc. To meet these demands and water needs, desirable practices were put in place to compensate for the permanent and irreversible reduction without limiting pumping from an aquifer (Akbarpour et al. 2020).

Due to the significant reduction of head in groundwater aquifers and their negative balance, optimizing water abstraction will be useful in the current condition, which first requires the knowledge of groundwater flow behavior. One way to determine groundwater flows behavior is to use mesh-dependent numerical models, such as finite element and finite difference methods to solve the governing equation of groundwater flow.

The following are the studies that have been done so far on the simulation of groundwater flow using these methods. A simulation model was designed and calibrated using MODFLOW software for unsustainable groundwater flow conditions in the Arhis-Nekour coastal aquifer in Morocco. The hydrodynamic coefficients of this aquifer were optimized.

A reservoir lumped model was proposed to simulate the long-term piezometric level process and the collection of solutes at the watershed surface. This model was used for aquifers in southern India. The simulation shows that the solute-in-water cycle may have significant negative pressure on groundwater, particularly in aquifers that are located in semi-arid areas where the main source of irrigation is groundwater. The 35-year-old forward simulation also states that long-term groundwater irrigation is not suitable in this area (Prinos et al. 2002).

To prevent the decrease of groundwater level and damage to the aquifer in the west of Egypt, and to obtain the optimal amount of pumping as well as the optimal number of wells in the region, genetic algorithm was used along with the MODFLOW model. In this study, by considering different scenarios, the groundwater status of the region was investigated. Finally, due to the scenarios and selecting the best scenario, the optimal amount of pumping from the wells in the region was determined (Moharram et al. 2011).

A groundwater flow model was developed for the Egyptian Delta Wadi. This conceptual model was a tool to analyze the hydrogeological system to understand how water flows in and out of the groundwater table. As the sampling was performed locally, the conceptual model was performed in its closest state to natural conditions. Groundwater modeling system (GMS) was provided. Then, using the finite difference method, the mathematical model of aquifer in three dimensions and the interphase estimated the data deficiencies and recorded factors such as the boundaries of aquifer using the parameter estimation code of surface geology, plant density and hydrological budget of the region as input.

After calibrating the model and achieving the desired results, considering the compatibility of the model with the natural conditions of the table, the heterogeneities in the aquifer were identified (Abdelaziz & Bakr 2012).

Another study examined a simulation model of a semi-arid aquifer in Morocco. MODFLOW model was used to simulate the desired aquifer and the values of hydraulic conductivity and storage coefficient were calculated by trial and error after manual calibration of the model (Le Page et al. 2012).

Software such as GMS using the MODFLOW package allows users to perform quantitative and qualitative modeling for porous media, including saturated and unsaturated zones, due to their study area conditions. The numerical method on which the MODFLOW package is based is the finite difference method, which divides the solution domain into rectangular grids of equal dimensions to begin its calculations. This method is not flexible on domains with irregular geometry because it is limited to rectangular grids and generates erroneous results.

In the last 5 years, researchers have used mesh-free numerical methods to model groundwater flow. In 2017, a study was conducted to model groundwater flow for the Birjand unconfined aquifer in steady and transient states, and the results were compared with the GMS solutions (Mohtashami et al. 2017a, 2017b). The higher accuracy of the mesh-free method was a proof to the fact that this method provides closer results to real condition (Mohtashami et al. 2017a, 2017b). Moreover, it is more useful to use optimization algorithms for planning and decisions based on the mesh-free model.

In a study, a combination of two genetics optimization algorithms and dynamic differential programming was used to solve the layout problem and determine the pumping rate and the optimal number of extraction wells to meet water needs. The cost of constructing a well was considered as a fixed cost and the cost of pumping was considered as time-dependent. Their results showed that applying a fixed cost to the objective function has a significant impact on determining the optimal number of wells (Hsiao & Chang 2020).

The problem of layout and determination of the optimal pumping rate of wells in the coastal aquifer was investigated. To simulate the influx of saline water, the numerical method of the boundary element was also used. The optimization method in this study was genetic algorithm (GA), which was combined with the simulation code to present an optimization–simulation model (Katsifarakis & Petala 2006). A simulation–optimization model was also proposed to manage groundwater resources. In this study, GMS software model is used for simulation and Meta-Heuristic Harmony Search algorithm is used for optimization. Three objective functions were investigated for two different modes, and sensitivity analysis was conducted on the parameters of the Harmony Search model (Ayvaz 2009).

In the study of an integrated model, ant colony simulation and optimization in continuous media for optimal management of coastal aquifers were presented. The objective function of their optimization model was to maximize the amount of pumping from wells in the coastal aquifer to prevent the saline seawater intrusion to them (Ketabchi & Ataie-Ashtiani 2011). In another study, using a numerical method of analytical element, groundwater flow in a real aquifer was modeled and then an optimization model based on particle swarm algorithm was proposed to maximize well pumping rate and minimize the cost of building a well system on the simulation model (Gaur et al. 2011).

Combining the support vector machine model and the particle swarm optimization (PSO) model, a simulation–optimization model was introduced. This model was used to determine the appropriate location and pumping rate in aquifer bioremediation operations. In this research, a number of selected wells were considered, and the best location as well as the optimal pumping rate were obtained. With linking an analytical method and a PSO algorithm, an optimization–simulation model was presented and an optimization model was written to maximize the water pumping rate. The aquifer used in this study was a hypothetical two-dimensional aquifer with rectangular geometry (El-Ghandour & Elsaid 2013).

A simulation–optimization model was presented to derive an optimal operation program in a hypothetically confined, heterogeneous and anisotropic aquifer with irregular geometry. The purpose of optimization in this study is to determine the number of production wells so that the amount of groundwater head drop is minimized according to the rules and costs related to water abstraction. Numerical finite-element method is used for simulation and particle clustering algorithm is used for optimization (Cyriac & Rastogi 2016).

In a study, evolutionary algorithms were used to provide an optimal management program to solve problems related to groundwater resources in coastal aquifers. Their optimization algorithms were ant colony optimization, PSO and GA (genetic algorithm), respectively. They also used the saturated–unsaturated transport model to simulate groundwater flow (Ketabchi & Ataie-Ashtiani 2015).

Multiple search method of adaptive GA (amalgam) was used to optimize the position of wells and pumping rate. In this study, the pumping rate along with three minimization objectives, namely the minimum deficit due to non-supply, changing the deficit index and minimizing the amount of depletion in the designated areas, was selected to determine an optimal solution for depletion and groundwater recharge. Their simulation model was finite difference in GMS environment. The results showed that the optimization–simulation method was able to calculate a set of optimal solutions that were portrayed on a Pareto front (Sadeghi Tabas et al. 2016).

Pathania et al. (2020) proposed a new simulation–optimization model named BIOEFGM-PSO to attain the optimal design of in-situ bioremediation system. Their model finds the optimal well locations by linking the element free Galerkin simulation method (EFGM) and PSO algorithm. They applied their model to a well-known problem and also a field problem. They achieved satisfactory results for both problems. Mohtashami et al. (2021a, 2021b) used particle filter-meshless local Petrov–Galerkin model to find the hydrodynamic coefficients of three standard aquifers. Their calibration-simulation model performance was strong and accurate due to its estimation. They recommend their model for use in all aquifers.

Researchers usually use simulation–optimization models to solve groundwater resource management problems. These models are usually used to evaluate the solution of differential equations which relate to groundwater flow that comes with an optimization model (Bear 1979; Ayvaz & Karahan 2008; Gaur et al. 2011). These optimizations can include maximum pumping usage, appropriate restrictions on permitted depletion of the aquifer, total pumping costs, and other limiting elements.

There are some studies in the field of groundwater simulation with the meshless methods. Mategaonkar & Eldho (2011) used the polynomial point collocation (PPCM) method to simulate groundwater flow in three aquifers. Two of them were one-dimensional, and one of them was two-dimensional. Their achieved results showed the higher accuracy of PPCM than mesh-based methods, e.g. finite element and finite difference. Thomas et al. (2014) engaged the point collocation method (PCM) for modeling groundwater flow in two standard aquifers. Their simulation process was carried out under unsteady state. They indicated that PCM has more satisfactory results than finite-element solutions. Mohtashami et al. (2017b) used the meshless local Petrov–Galerkin (MLPG) method to simulate groundwater flow in an aquifer. They compared their results with the finite difference method and found their model more accurate than finite difference. Boddula & Eldho (2018) developed a simulation–optimization model to investigate three different scenarios, including maximization of the total pumping volume, minimization of the pumping cost with respect to the given demand, and minimization of cost to establish a new pumping well system. Their simulation and optimization models were MLPG method and modified artificial bee colony, respectively.

In this study, by adding PSO model to the MLPG groundwater simulation model, a simulation–optimization model is presented to find the appropriate and optimal location of extraction wells so that the drawdown in the aquifer is minimized. To do this, in the unconfined Birjand aquifer, groundwater is first simulated and then the optimization model is added to it.

The Birjand aquifer study area, with approximate coordinates 32°34′00″N to 33°08′00″N, and 58°41′00″E to 59°44′00″N, is located in South Khorasan province in Iran. This aquifer, which is of unconfined type with an approximate area of 269 km² and an average saturation thickness of 30 m, is marked in red in Figure 1 (Sadeghi Tabas et al. 2016). This figure also shows the geographical location of the aquifer in Iran. The annual precipitation of the aquifer is 170 mm and it has nine inflow pathways and one outflow pathways, as shown in Figure 1. Based on the Domarton climate index this region is counted as an arid region and, due to the high depth of groundwater, evaporation value is zero. (Hamraz et al. 2015) Maximum and minimum elevations are 2,736 and 1,167 m, respectively. The slope is deep in the east part of the plain and is more gentle toward the west. The bedrock of the aquifer is mostly hard rock, such as sandstone, conglomerate, and tuff formations. Most of the area in the east is young gravel fans and clay flat extends from the east to the center of the aquifer. The area in the west consists of young gravel fans and gravel plains with a saline flat. Low-density pasture is the major land-use land-cover category in the study region.

The MLPG method is one of the real mesh-free methods because meshing is not required at any stage of analysis, including the approximation of field variable and the numerical integration process. Atluri & Zhu (1998) first proposed this method which is known as the weak form, in 1998. Two functions are used in this numerical method: approximation function and weight function. Moving least squares (MLS) approximation function helps MLPG to compute the shape function denoted ⁠, and cubic spline weight function calculates the value of weights for each node. The following sections explained these two functions in detail. Gaussian integration is used to solve the integral equations (Liu & Gu 2005).

Due to its independency to meshing, this method removes drawbacks and errors generated from gridding process. Meshes in finite difference method limits to rectangular grids. This indicates the weakness of the finite-element method in domains with irregular geometry. In the finite-element method, the computational cost for gridding is high, which makes it inappropriate for domains with large areas and complicated boundary conditions. However, in the MLPG method, all these drawbacks are removed.

This function was first used in 1992 in the field of solid mechanics (Atluri & Zhu 1998). In this function, the field variables are calculated by multiplying polynomials into dimensionless coefficients. Two main features of this function have increased its use: (1) the function approximated by this method is continuous in the problem domain; (2) it has the ability to make approximations to the desired order of consistency. It was used in this study to construct shape function.

In general, the vector P(X) is based on the Pascal triangle (Figure 2; Le Page et al. 2012).

The vector of coefficients in Equation (1) is determined using the function values in a set of nodes located in the support domain of the point x. The number of nodes used locally to estimate the value of the function at point x is determined by the support domain.

For example, in Figure 3, the number of points located in the support domain is four.

To enforce essential boundary conditions (known as constant head boundaries) in this study, the penalty approach was used. As the MLS approximation function does not satisfy the Kronecker Delta property, the penalty technique was engaged. This function provides MLS approximation function to impose essential boundary conditions without deficiency (Liu & Gu 2001). The function, which is developed by Atluri & Zhu (1998), is efficient and does not need any other additional unknown variables.

Choice of weight function is important in the performance of approximation function (MLS). The engaged weight function must have at least three following properties (Liu & Gu 2005): (1) the value of weight function must be positive; (2) this value must be zero outside the computational domain; and (3) W(x) must be decreased monotonically from the interest point x.

The first and second term in the right side of Equation (37) represent the concentrated flow rate (wells) and the distributed flow rate (precipitation or evaporation).

The particle swarm algorithm is an optimization technique based on laws of probability, first presented by Russell Eberhart, a computer scientist, and James Kennedy, a social psychologist, in 1995 (Kennedy & Eberhart 1995). PSO is based on the behavior of collective components such as the movement of fish and birds in nature. This algorithm searches the space of the objective function with each particle or component. The motion of each particle has a definite and probable component. The algorithm of this method first starts with a set of components, each of which is an answer to the problem and are randomly selected.

The particle swarm algorithm steps are as follows:

• 1.

  By considering a population of components at random (solutions), the velocity values as well as the location of the components are determined. Then, after determined iteration, the search begins to find the optimal condition.

• 2.

  During each iteration, the competency value of each component is calculated.

• 3.

  If the calculated competency value for the particle is the best value, the location of the component is stored as ⁠.

• 4.

  In each iteration, among all the components, the component with the best value is selected and the component location is saved as ⁠.

• 5.

  Velocities are corrected based on the position of and ⁠.

• 6.

  The position of the components is fixed.

• 7.

  If the termination condition is met, such as maximum iteration or minimum error, the operation terminates. Otherwise, it will be repeated from step 2.

These steps are known as standard PSOs and apply only to continuous problems, which can also be applied to the problem in this study.

Boundary conditions in groundwater are divided into two categories: constant head and constant flow. The first is the Dirichlet boundary, where the groundwater level remains constant over time and the second is the Neumann boundary condition. At boundary points with a constant head, the groundwater level for the boundary is determined and its value remains constant along the boundary. These boundaries include nine inflow pathways and one outflow pathway in the form of drainage.

The most important factor in the decrease of groundwater level in the aquifer is abstraction through production wells. There are 190 wells in this area, which are shown in Figure 4. There are also ten observation wells that stand out in red from the rest of the wells.

Due to the region's arid climate and low rainfall, the same small amount of rainfall is considered as the amount of recharge in the aquifer. The amount of rainfall is 0.0000132 m per day; this value was extracted based on the data of Birjand plain rain-gage station in 2012–2013. The amount of water abstracted per unit of time from wells is considered as the amount of discharge in the model.

Based on the results of pumping tests in this aquifer, the hydraulic conductivity coefficients and specific yield values range [1.11, 75.09] and [0.01, 0.37], respectively. These values depend on the properties of the soil in the aquifer. These ranges are considered in Sadeghi Tabas et al. (2016) research, which calibrated the values using generalized likelihood uncertainty estimation (GLUE) method. In their study, this aquifer has been divided into several zones, each of which shows the numerical areas as hydraulic conductivity and specific yield. Figure 5 shows these calibrated values.

The geometry of the Birjand aquifer was modeled using MATLAB software, and the nodal points were evenly distributed so that the distance between two consecutive points in horizontal and vertical directions was 500 m (Mohtashami et al. 2019a, 2020a). The distance was selected based on previous studies in the Birjand plain aquifer Sadeghi Tabas et al. (2016). The four layers of information for constructing the groundwater flow model in the Birjand aquifer were boundary conditions, discharge rates of wells, specific yield, and hydraulic conductivity values.

Three error indices were calculated. MLPG has a lower error than the finite element method, in other words the MLPG method was more accurate than the MODFLOW solution. The results in Table 1 were computed when the weight function was cubic spline. For quartic spline and exponential weight functions, our model indicates results with meaningful errors, such that the root mean square error (RMSE) is 3.32 and 5.21 m, respectively.

This study aims to find the optimal location of extraction wells with the aim that the aquifer has the least drawdown, and a number of existing wells were removed to have better and more understandable analysis of the results. The selection of removed wells was completely random and the number of removed wells was 94. The remaining wells are shown in Figure 6.

The groundwater level with 96 wells in the aquifer was modeled and the results for piezometers 1, 5, and 10 were taken in the form of outflow diagrams, which are shown in Figure 7.

Figure 7 shows the groundwater level over a 1-year period, with monthly time steps for piezometers 1, 5, and 10 of the aquifer. According to the graphs, they all have a decreasing trend, such that the amount of drawdown increases. The diagram for piezometer 1 shows a drop of 0.007 m. Moreover, during a modeling period, the groundwater level decreases with a gentle slope. In the diagram for piezometer 5, the drop in groundwater level reaches about 0.6 m, and the high value of this drop is due to the proximity of this observation well to the pumping wells, which reduces the water level in this area. In the groundwater level diagram related to piezometer 10, due to the low density of wells in that area, the groundwater level drop is about 0.1 m, which decreased significantly in the third month so that the chart suddenly descends and then it continues its decreasing trend with a gentle slope.

The optimal location of extraction wells was determined. For this purpose, the PSO model was connected to the simulation model and implemented. The time taken to run this model each time was about 6 h, which was done on a home computer with the following specifications CPU i5/5300 U/8 Gb RAM. Figure 8 shows the optimal location of production wells. The conditions that are observed in the optimization model are: (1) the wells are not placed on top of each other, and (2) the wells do not fall on the border.

Figure 8 shows the optimal location of extraction wells so that the drop in the aquifer is minimized. The wells have moved to areas of the aquifer that are close to constant head boundaries. The reason for this is the high groundwater level in these areas over time. Therefore, by boring wells in these areas, there is no high drop in aquifer piezometers. The concentration of wells in the southwestern part of the aquifer is also an important point, and may have been due to the region's bedrock conditions and topography.

The sensitivity analysis of the effective parameters in the PSO model was examined and different values of C1 and C2 coefficients are evaluated. Table 2 shows the different values of the C1 and C2 coefficients in calculating the total groundwater level difference at the beginning and end of the modeling cycle in piezometers.

According to Table 2, in all cases the value of the objective function is less than when only the simulation model is executed. The appropriate value of the coefficients C1 and C2, such that the objective function is at its minimum, is related to row 1 of Table 2, where C1 and C2 are equal to 1.5 and 2, respectively. Figure 7 is based on row 2 values.

Figure 9 shows groundwater level in piezometers 1, 5, and 10 after applying the optimization model.

As can be seen in Figure 9, the groundwater level drop at the beginning and end of the modeling period in the piezometers has decreased significantly. This means that the optimization model has worked properly and has found the optimal location of the wells in a way that causes the least amount of drop. In the piezometer diagram in Figure 7, the drop in water level at the beginning and end of the cycle is 0.002 m, while in Figure 6 it is 0.007 m. The drop in groundwater level in piezometer 5 in Figure 7 is 0.04 m, and this value is 0.6 m in Figure 6. Therefore, the optimization model has correctly reduced the groundwater level drop.

The accurate prediction of groundwater behavior helps the decision makers to manage this aquifer more precisely and have pervasive knowledge about the quantity, volume and groundwater head of the aquifer for now and the future. This model also helps the managers prevent the formation of land subsidence due to extraction of wells.

In the present study, a simulation–optimization model is presented to determine the optimal position of well excavation so that the drawdown is minimized. The MLPG method as simulation model was used. The groundwater level in a real aquifer was modeled in transient state and compared and verified with the results of finite difference numerical method. Then, by adding the PSO model to the simulation model, the optimal well location was found. In the Birjand aquifer, with the implementation of the optimization-simulation model, the groundwater level drop was associated with a 47% decrease; in other words, the total groundwater level drop decreased from 2.808 m to 1.329 m. The drawdown value was computed for three piezometers before and after enforcement of the optimization model. The results showed a significant reduction. For instance, in the first piezometer, the drawdown was 0.007 and 0.003 m before and after model enforcement, respectively. Furthermore, the sensitivity analysis of optimization algorithm parameters was performed, and the best values of C1 and C2 were calculated, which were 1.5 and 2, respectively. Finally, the effect of optimization model on water level drop in piezometers was investigated.

The authors thank Ali Mohtashami for his insightful comments during the preparation of this paper.

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