In the present study, the optimal place to excavate extraction wells as the drawdown is minimized was investigated in a real aquifer. Meshless local Petrov–Galerkin (MLPG) method was used as the simulation method. The closeness of its results to the observational data compared to the finite difference solution showed the higher accuracy of this method with a root mean square error for MLPG of 0.757 m compared to the finite difference which equaled 1.197 m. Particle warm algorithm was used as the optimization model. The objective function was defined as the summation of the absolute values of difference between the groundwater level before abstraction and the groundwater level after abstraction from wells. In the Birjand aquifer, which was investigated in transient state, the value of objective function before applying the optimization model was 2.808 m, while in the optimal condition it reached 1.329 m (47% reduction in drawdown). This was investigated and observed in three piezometers. In the first piezometer, the drawdown before and after model enforcement was 0.007 m and 0.003 m, respectively. This reduction occurred in other piezometers as well.

  • In the present study, using a simulation–optimization model, the best place to bore a well in a real aquifer under transient conditions was discussed.

  • The meshless local Petrov-Galerkin method as simulation model was used.

  • The water level in a real aquifer was modeled in transient state and compared and verified with the results of finite difference 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.

Case study

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 km2 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.

Figure 1

Geographical location of Birjand aquifer and plain (Sadeghi Tabas et al. 2016).

Figure 1

Geographical location of Birjand aquifer and plain (Sadeghi Tabas et al. 2016).

Figure 2

Pascal Triangle (Liu & Gu 2001).

Figure 2

Pascal Triangle (Liu & Gu 2001).

Meshless local Petrov–Galerkin

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.

MLS approximation function

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.

If is a function of field variables in the area Ω under study, the approximation is shown at point with. MLS approximation describes locally form of the field and is expressed in Equation (1):
formula
(1)
where the number of monosyllabic constituents and is a vector of coefficients , which is defined in Equation (2):
formula
(2)
In Equation (1), P(X) is a vector of basic functions, which often contains the maximum number of monosyllables needed to achieve the minimum completeness. In one-dimensional space, a complete polynomial base of order m is expressed by Equation (3):
formula
(3)
and in two-dimensional space :
formula
(4)

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.

Figure 3

Hypothetical amplitude with support domain around the desired point.

Figure 3

Hypothetical amplitude with support domain around the desired point.

To determine the unknown coefficients, the weighting function is first defined as follows (Moharram et al. 2011):
formula
(5)
In the above equation, represents the weight function of the nodal point I and the value in brackets is the difference between the value estimated at point I and the value given at the same point. Also, n is the number of points in the support domain of the function (weight). The weight function has non-zero values in its support domain (Le Page et al. 2012). To minimize the function J, the condition of Equation (6) is examined:
formula
(6)
which eventually leads to a linear Equation (7):
formula
(7)
In Equation (7), , and are defined in Equations (8), (9) and (11), respectively:
formula
(8)
formula
(9)
formula
(10)
formula
(11)
By placing Equation (7) in Equation (1), the MLS approximations are given as Equations (12) and (13):
formula
(12)
formula
(13)
where is the approximation of the function, is the shape function and is the nodal parameter.
In other words, the shape function is expressed as follows:
formula
(14)

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.

Weight function

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.

In this study, three different weight functions–cubic spline, quartic spline, and exponential function – are evaluated. The formulas of these functions are presented, respectively, in the following equations (Liu & Gu 2005):
formula
(15)
formula
(16)
formula
(17)
In these equations, is a constant shape function and and is the influence radius of node. The value of is computed based on the Equation (18) (Liu & Gu 2005):
formula
(18)
where is a constant parameter range [2, 3] (Liu & Gu 2005). In this study, we found the best value of equals 2.1 and leads to good results based on the sensitivity analysis for this parameter.

Formulation of governed equation with MLPG

General form of the governing equation in groundwater modeling in unconfined aquifer is:
formula
(19)
where H is groundwater head, and are the hydraulic conductivity coefficients.
As the following equation:
formula
(20)
with substitution (20) in (19), we have:
formula
(21)
with assumption of isotropic condition in the Birjand aquifer ():
formula
(22)
formula
(23)
using weighted residual method:
formula
(24)
Using integration of parts, Equation (25) is achieved:
formula
(25)
As there is no constant flow boundary condition in the Birjand aquifer:
formula
(26)
with substituting (20) in (26):
formula
(27)
As we have the following equation in numerical methods, e.g. finite element and meshless methods:
formula
(28)
Substitution (28) in (27):
formula
(29)
In this study, finite difference scheme in time is used. For computing time derivative, the following forward implicit finite difference approximation is used (Mohtashami et al. 2019b, 2020b):
formula
(30)
With substitution Equation (30) in (29):
formula
(31)
formula
(32)
formula
(33)
This equation shows a set of linear equations in the form of KU = F, where (Mohtashami et al. 2021a, 2021b):
formula
(34)
formula
(35)
formula
(36)
represents force body matrix, and:
formula
(37)

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).

PSO algorithm

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 algorithm starts by considering a group of random answers, then searches for the optimal solution in the problem space by updating the position and velocity of the particles. Each particle is identified by two values, x and v, which represent position and velocity. At each step of the set motion, each component or particle is updated through the two best values. The first value is the best answer in terms of competence for each particle separately – this value is called. The second value is the best value among all the particles. This value is called . The new velocity and position of the particles are updated at each stage by the following equation (Stacey et al. 2003):
formula
(38)
formula
(39)
where and are cognitive coefficients that are determined according to the higher priority or whether it moves more toward its own experiences or toward the experiences of the model person. is usually considered insert equal or less than symbol 4 in implementations. and are also random numbers obtained from a uniform distribution between 0 and 1. is known as the inertia coefficient and is ignored. It should be noted that w is in the range of 0.4–0.9. Also, is the k-th location of the particles and is their previous location.

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.

The objective function

In this study, after the pumping wells have taken the position at the aquifer, the groundwater level at all nodes in the aquifer is calculated. Then, having the previous level of groundwater head, its difference with the calculated level can be estimated. In this study, the objective function is defined as a way to minimize the difference value between the calculated groundwater level and the initial level in the aquifer. Equation (40) states the objective function (Akbarpour et al. 2020):
formula
(40)
where Z denotes the objective function, N Node is the number of nodes distributed in the aquifer, H stands for the groundwater level before the abstraction of pumping wells, and shows the head of groundwater after abstraction. The considered constraints for this problem are: (1) the wells must be in different locations; (2) groundwater level in each node must have a non-zero value; and (3) the wells cannot fall on the boundaries and edges.

Boundary conditions

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.

Extraction and observation wells

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.

Figure 4

Extraction wells and observation wells.

Figure 4

Extraction wells and observation wells.

Recharge parameter

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.

Aquifer hydrodynamic coefficient

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.

Figure 5

Calibrated values of hydraulic conductivity (k) and specific yield coefficient (Sy) (Sadeghi Tabas et al. 2016).

Figure 5

Calibrated values of hydraulic conductivity (k) and specific yield coefficient (Sy) (Sadeghi Tabas et al. 2016).

Creating conceptual model

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.

This aquifer has already been modeled, calibrated, and validated by Mohtashami et al. (2017a, 2017b) in both steady and transient state. They found that this model performs more accurately than traditional methods such as finite difference. Table 1 shows the three indices computed from the following equations:
formula
(41)
formula
(42)
formula
(43)
where and indicates the observed and simulated heads, respectively, while m is the number of monthly time steps and n is the number of piezometers.
Table 1

Calculation of mean error, absolute mean, and RMSE in transient state

Finite difference method (Farpour et al. 2018)MLPG method
0.159 −0.08 Mean error (meter) 
1.434 0.573 Absolute mean error (meter) 
1.197 0.757 Root mean square error (meter) 
Finite difference method (Farpour et al. 2018)MLPG method
0.159 −0.08 Mean error (meter) 
1.434 0.573 Absolute mean error (meter) 
1.197 0.757 Root mean square error (meter) 

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.

Figure 6

Extraction wells (96) in the aquifer.

Figure 6

Extraction wells (96) in the aquifer.

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

Groundwater level in piezometers 1, 5, and 10.

Figure 7

Groundwater level in piezometers 1, 5, and 10.

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

Optimal location of production wells in Birjand aquifer.

Figure 8

Optimal location of production wells in Birjand aquifer.

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.

Table 2

Sensitivity analysis of C1 and C2 coefficients

fC2C1RowfC2C1Row
1.764 1.604 
1.381 1.5 1.329 1.5 
1.965 0.5 1.571 
fC2C1RowfC2C1Row
1.764 1.604 
1.381 1.5 1.329 1.5 
1.965 0.5 1.571 

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.

Figure 9

Water level obtained in piezometers after applying the optimization model

Figure 9

Water level obtained in piezometers after applying the optimization model

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.

Abdelaziz
R.
&
Bakr
M. I.
2012
Nverse modeling of groundwater flow of Delta Wadi El-Arish
.
Water Resource and Protection
4
(
7
),
432
438
.
Akbarpour
A.
,
Zeynali
M. J.
&
Nazeri Tahroudi
M.
2020
Locating optimal position of pumping wells in aquifer using meta-heuristic algorithms and finite element method
.
Water Resources Management
34
,
21
34
.
Bear
J.
1979
Hydraulics of Groundwater
.
The University of Michigan
,
McGraw-Hill International Book Co
.
Cyriac
R.
&
Rastogi
A. K.
2016
Optimization of pumping policy using coupled finite element-particle swarm optimization modelling
.
ISH Journal of Hydraulic Engineering
22
(
1
),
88
99
.
El-Ghandour
H. A.
&
Elsaid
A.
2013
Groundwater management using a new coupled model of flow analytical solution and particle swarm optimization
.
International Journal of Water Resources and Environmental Engineering
5
(
1
),
1
11
.
Farpour
A.
,
Ramazani
Y.
&
Akbarpour
A.
2018
Numerical simulation of chromium changes trend in aquifer of birjand plain
.
Iranian Journal of Irrigation and Drainage
12
(
5
),
1203
1216
.
Hamraz
B. S.
,
Akbarpour
A.
,
Pourreza Bilondi
M.
&
Sadeghi Tabas
S.
2015
On the assessment of ground water parameter uncertainty over an arid aquifer
.
Arabian Journal of Geosciences
8
,
10759
10773
.
Hsiao
C. T.
&
Chang
L. C.
2020
Dynamic optimal groundwater management with inclusion of fixed costs
.
Water Resources Planning and Management
128
(
1
),
57
65
.
Katsifarakis
K. L.
&
Petala
Z.
2006
Combining genetic algorithms and boundary elements to optimize coastal aquifers management
.
Hydrology
1–2
(
327
),
200
207
.
Kennedy
J.
&
Eberhart
R.
1995
Particle swarm optimization
. In:
Proceeding of IEEE Intenational Conference on Neural Networks
.
Ketabchi
H.
&
Ataie-Ashtiani
B.
2011
Development of combined ant colony optimization algorithm and numerical simulation for optimal management of coastal aquifers
.
Iranian Water Resources Research
7
(
1
),
1
12
.
Le Page
M.
,
Berjamy
B.
,
Fakir
Y.
,
Bourgin
F.
,
Jarlan
L.
,
Abourida
A.
,
Benrhanem
M.
,
Jacob
G.
,
Huber
M.
,
Sghrer
F.
,
Simonneaux
V.
&
Chehbouni
G.
2012
An integrated DSS groundwater managment based on remote sensing. The case of semi-arid aquifer in Morocco
.
Water Resources Management
26
,
3209
3230
.
Liu
G. R.
&
Gu
Y. T.
2001
A point interpolation method for two-dimensional solid
.
International Journal of Numerical Methods in Engineering
50
,
937
951
.
Liu
G. R.
&
Gu
Y. T.
2005
An Introduction to Meshfree Methods and Their Programming
.
Springer
,
Singapore
.
Mategaonkar
M.
&
Eldho
T. I.
2011
Simulation of groundwater flow in unconfined aquifer using meshfree point collocation method 2011
.
Engineering Analysis with Boundary Elements
35
(
4
),
700
707
.
Moharram
S. H.
,
Gad
M. I.
,
Saafan
T. A.
&
Khalaf Allah
S.
2011
Optimal groundwater management using genetic algorithm in El-Farafra oasis, Western Desert, Egypt
.
Water Resources Management
26
,
927
948
.
Mohtashami
A.
,
Akbarpour
A.
&
Mollazadeh
M.
2017b
Modeling of groundwater flow in unconfined aquifer in steady state with meshless local Petrov–Galerkin
.
Modares Mechanical Engineering
17
(
2
),
393
403
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2019a
Determination the capture zone of wells by using meshless local Petrov-Galerkin numerical model in confined aquifer in unsteady state (case study: Birjand Aquifer)
.
Iranian Journal of Ecohydrology
6
(
1
),
239
55
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2019b
Prediction of groundwater fluctuations using meshless local petrov-galerkin numerical method in a field aquifer (Birjand Aquifer)
.
International Journal of Numerical Methods in Civil Engineering
3
(
4
),
33
41
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2020a
Computation of groundwater balance using numerical MLPG method (case study: birjand unconfined aquifer)
.
Iranian Journal of Irrigation & Drainage
14
(
4
),
1460
1474
.
Mohtashami
A.
,
Hashemi Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2021a
Estimation of parameters in groundwater modeling by particle filter linked to the meshless local Petrov-Galerkin numerical method
.
Journal of Hydraulic Structures
7
(
1
),
16
37
.
Mohtashami
A.
,
Monfared
S. A.
,
Azizyan
G.
&
Akbarpour
A.
2021b
Usage of particle filter for exact estimation of constant head boundaries in unconfined aquifer
.
AmirKabir Journal of Civil Engineerin
53
(
11
),
1
21
.
Prinos
S.
,
Lietz
A.
&
Irvin
R.
2002
Design of A Real-Time Groundwater Level Monitoring Network and Portrayal of Hydrologic Data in Southern Floria
.
US Geological Survey Report
,
Tallahassee, FL
.
Sadeghi Tabas
S.
,
Samadi
S. Z.
,
Akbarpour
A.
&
Pourreza Bilondi
M.
2016
Sustainable groundwater modeling using single-and multi-objective optimization algorithms
.
Hydroinformatics
18
(
5
),
1
18
.
Stacey
A.
,
Jancic
M.
&
Grundy
I.
2003
Particle swarm optimization with mutation. In:
The 2003 Congress on Evolutionary Computation. CEC'03. 2: 1425–1430, IEEE
.
Thomas
A.
,
Eldho
T. I.
&
Rastogi
A. K.
2014
A comparative study of point collocation based meshfree and finite element methods for groundwater flow simulation
.
ISH Journal of Hydraulic Engineering
20
(
1
),
65
74
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).