## Abstract

The main objective of this study is to propose a linked simulation–optimization approach for solving aquifer parameter estimation problems from pumping test results. In the simulation part of the proposed approach, the drawdowns at given monitoring points and times are calculated by considering the Theis and Hantush approaches for confined and leaky confined aquifers, respectively. This simulation part is then integrated to a newly proposed hybrid optimization approach, namely DE–Solver, which integrates the differential evolution (DE) algorithm and generalized reduced gradient (GRG) method of the spreadsheet Solver add-in. The performance of the proposed approach is evaluated by considering two pumping test data for confined and leaky confined aquifers. Identified results indicate that the proposed approach provides better results than those obtained by using different approaches in the literature.

## INTRODUCTION

Groundwater management is an important part of water resources planning and management. Mathematical simulation models are the essential tools to develop sustainable management plans for groundwater resources, which require knowledge about the spatial distributions of aquifer parameters. These parameter distributions are usually obtained by interpolating point observations, which are mostly determined from pumping test results. During these tests, pumping rates are kept constant and corresponding drawdowns in different monitoring locations are recorded. After obtaining time-drawdown data, parameters of the aquifer system are determined by using the manual curve-matching approaches. Although these approaches are simple to employ, some errors may occur since the accuracy of these approaches is mostly dependent on the modeller's ability (Huang *et al.* 2008). Therefore, simulation–optimization models have become popular tools to determine aquifer parameters.

The current literature includes various solution approaches integrating different simulation and optimization models for determining aquifer parameters. Although these models work in a similar manner, they have some differences in terms of the considered optimization approaches. Note that use of heuristic optimization approaches such as simulated annealing (SA) (Huang *et al.* 2008), genetic algorithm (GA) (Samuel & Jha 2003), and particle swarm optimization (PSO) (Şahin 2018) algorithms is preferred for solving the aquifer parameter estimation problems due to the non-convex solution spaces of groundwater optimization problems. However, stand-alone use of heuristic approaches may require long computation times to precisely find the global optimum solutions (Michalewicz 1992). Therefore, hybridized versions of them are employed for solving the associated optimization problems to reduce the computation times. The main philosophy of the hybrid optimization approaches is to integrate the global exploration feature of the heuristic approaches and local search ability of the deterministic optimization approaches. In these algorithms, the global exploration process starts with multiple starting points and explores the search space, and then, deterministic approaches find the optimum solution by taking the results of global exploration as their initial values (Ayvaz *et al.* 2009). This kind of an integration makes the solution of the problem more robust than either the heuristic and deterministic optimization approaches by themselves (Land 1998). Using this strategy, Ayvaz & Gurarslan (2018) proposed a solution methodology where the hybrid HS–Solver optimization approach is used. HS–Solver is a recently proposed hybrid optimization approach which integrates the heuristic harmony search (HS) algorithm and a generalized reduced gradient (GRG) approach in the spreadsheet Solver add-in. In this integration, while HS works as the global optimizer, the Solver add-in is used for local search to precisely find the global optimum solution.

As an extension of Ayvaz & Gurarslan (2018), the main objective of this study is to solve aquifer parameter estimation problems from pumping test results by proposing a new hybrid optimization approach. This optimization approach is called DE–Solver, which integrates the differential evolution (DE) algorithm and Solver add-in as the global and local optimizers, respectively. In the simulation part of the proposed approach, the flow process in the subsurface is simulated by using two analytical solution models based on the Theis and Hantush methods for confined and leaky confined aquifers. These analytical models are then integrated to DE–Solver in order to determine the aquifer parameters based on a hybrid optimization framework. The performance of the proposed DE–Solver approach is evaluated by considering two pumping test data for confined and leaky confined aquifers in the literature. The identified results indicate that the proposed hybrid DE–Solver optimization approach results in similar or better results than those obtained by using curve matching and other heuristic and hybrid optimization approaches.

## MODEL DEVELOPMENT

In this section, the main structure of the proposed hybrid approach is presented. For this purpose, mathematical formulations of the Theis and Hantush models are given first, and then, how to integrate them to the hybrid DE–Solver optimization approach is presented.

### Theis model

*r*and at a time

*t*,

*Q*is the pumping rate which is constant during the test,

*T*is the transmissivity of the aquifer,

*u*is a dimensionless parameter, and is the well function of

*u*, which is also known as the Theis function. The dimensionless parameter of

*u*and well function of are given as follows:where

*S*is the storage coefficient of the aquifer,

*t*is the time since the beginning of pumping, and

*y*is the variable of integration. Note that the well function in Equation (3) can also be expressed with the following series (Kresic 1997):In the proposed approach, the value of is calculated by means of Equation (4) by taking the upper limit of

*n*to be 100.

### Hantush model

*B*is the leakage factor, and is the well function for the leaky confined aquifer whose value is the function of

*u*and . Definitions of and

*B*are given as follows:where and are the thickness and hydraulic conductivity of the confining bed (aquitard). Note that values are tabulated for different

*u*and values by numerically integrating Equation (6) on MATLAB before executing the DE–Solver approach. In the optimization process, values of are determined by conducting a bilinear interpolation between these tabulated values for the calculated

*u*and values.

### Hybrid DE–Solver optimization approach

DE is an effective optimization procedure in determining the global optimum solution, and much faster than other evolutionary methods (Storn & Price 1997). Like other evolutionary algorithms, DE is also a population-based algorithm and utilizes like operators, namely crossover, mutation and selection just as given in GA. The basic computational steps of DE are as follows.

*p*be the number of candidate solutions,

*n*be the number of decision variables,

*G*be the generation index, be the vector which includes the values of the decision variables to be determined, and be the vectors which include the minimum and maximum bounds of the decision variables, respectively. Using these definitions, the initial value of the decision variable in the candidate solution at generation can be generated as follows:where represents a uniform random number within the range of . After generating the initial solutions, a new mutant vector of is generated from based on the mutation process. Note that a wide variety of mutation strategies are available to generate in the literature. One of the most widely used is the ‘DE/rand/1’ strategy, where the individuals of are randomly chosen, and one pair of solutions is selected as follows:where

*F*is a scaling factor which is used to control the strength of the mutation, and are the different integers which are randomly generated between for each candidate solution. After this process, a new solution vector of is generated based on the crossover process as follows:where is the crossover rate and is a randomly selected integer from the range of . Equation (10) states that

*n*uniform random numbers are generated and their values are compared with the value of . If or , then the component of is selected from , otherwise, it is selected from . The reason for defining the condition is to guarantee the selection of at least one value from the vector . As in the mutation process, several different crossover strategies are available and detailed information about these mutation and crossover strategies can be found in Mezura-Montes

*et al.*(2006) and Price

*et al.*(2005). After this process, candidate solutions for the next generation are determined by evaluating all members of the population as follows:where is the objective function value. As can be seen from Equation (11), the solution with the better objective function value is selected for the next generation, which is the result of the ‘survival of the fittest’ concept. Once the new generation is created, the solution sequence from Equations (8)–(11) is repeated until termination. Note that different termination criteria can be used to finalize the search process. Among them, termination after completing the specified number of generations is used in this study.

After describing the computational structure of DE, the next step is to integrate it to the spreadsheet Solver add-in for developing the new hybrid DE–Solver optimization approach. In this integration, DE first explores the entire solution space and finds potential locations where an optimum solution exists, and then, Solver precisely finds the optimum solution by using the DE results. This methodology makes finding the global optimum easier than by global and local search methods by themselves (Ayvaz *et al.* 2009; Kayhan *et al.* 2010).

Solver is a powerful gradient-based optimization add-in and is available in many commercial spreadsheet products such as Lotus^{®} and Microsoft Excel^{®}. It can solve nonlinear optimization problems based on the GRG approach. Note that the DE and GRG optimization processes are integrated by developing Visual Basic for Applications (VBA) codes on the background of the Excel^{®} spreadsheet platform. In this integration, the execution of Solver is controlled by using a small probability value of just as given in Ayvaz & Gurarslan (2018). The flowchart of the proposed DE–Solver optimization approach can be seen in Figure 1. As can be seen, once the new solution is found by DE, this value is taken as the initial solution for Solver. After completing the local search process, the final decision variable and objective function values are transferred to DE again. The required solution parameters of the proposed DE–Solver optimization approach are the number of candidate solutions , scaling factor , crossover rate , and probability of . After defining these parameters and executing DE–Solver, the problem of aquifer parameter estimation can be solved by using the proposed approach. Note that the sum of squared errors calculated between the observed (pumping test data) and simulated drawdown values is used as the objective function in the proposed approach. Decision variables of the optimization model are *T* and *S* for the confined and *T*, *S*, and *B* for the leaky confined aquifers.

## NUMERICAL APPLICATIONS

The performance of the proposed hybrid DE–Solver approach is evaluated by considering two pumping test data for confined and leaky confined aquifers, respectively. For both tests, the identified aquifer parameters are compared with the ones which are determined by using the curve-matching approach in Kresic (1997) and hybrid HS–Solver optimization approach in Ayvaz & Gurarslan (2018). The solution parameters are taken as: , , , and . Note that these values are selected based on previous experience and the recommendations of Prof. Storn (developer of the DE optimization algorithm) for solving different engineering optimization problems. On the DE algorithm web page (Storn 2019), Prof. Storn states that the parameter *p* can be used as ten times the number of variables. Empirically, a value of *p* greater than 40 does not improve the solution's quality. A selection of 0.80 for both *F* and can be a good choice to solve various optimization problems (Elçi & Ayvaz 2014; Gurarslan & Karahan 2015; Ayvaz & Gurarslan 2017). However, the values of both parameters can be adjusted dynamically in the range of [0.50, 1.00] for problems with noisy objective functions. Regarding the parameter of , there is no specific recommendation to determine it. This parameter is used to call the Solver from the DE optimization approach. Use of higher values makes the local search by DE more dominant and use of smaller values makes the global search by Solver more dominant. Since Ayvaz & Gurarslan (2018) solved both problems by considering 10,000 iterations, the number of DE generations is set to 500 to obtain comparable results.

### Example 1: confined aquifer

A 24-hour pumping test is conducted to determine *T* and *S* in a confined aquifer. The pumping rate of a single well *Q* is kept constant at 0.008 and the corresponding drawdowns at piezometers P1, P2, and P3 are recorded. The distance of each piezometer to the well location is m, m, and m. By applying the proposed approach for each piezometer, the results given in Figure 2 are obtained. As can be seen, the search process is started from the same initial solutions since the same random number seeds are considered for both DE and DE–Solver approaches. For each piezometer, a significant improvement is observed in the values with inclusion of local search by Solver. For these outcomes, Table 1 compares the model identification results with the ones obtained by using different solution approaches. It is clearly seen that the identified *T* and *S* values are similar to the ones obtained from the curve-matching approach for P2 and P3 although there is a small difference in both for P1. The reason for this difference is associated with the distance between the pumping well and P1. Since P1 is the closest piezometer to the well location, the recorded drawdown data of P1 are much less curved than for P2 and P3. Therefore, the matching of P1 and Theis curve is more arbitrary, as expected (Kresic 1997). Comparison of the pure heuristic and hybrid approach results indicates that inclusion of Solver in both HS and DE approaches results in a small change in the identified aquifer parameters. For these results, Table 1 also compares the model results in terms of the and function evaluation numbers. As can be seen, the calculated for P1 in the curve-matching approach is approximately ten times bigger than in the DE and DE–Solver approaches. The reason for this outcome is also associated with the less curved nature of P1. From the results, it is clear that inclusion of Solver as a local optimizer slightly improves the objective function values in the DE–Solver approach. Comparison of DE–Solver and HS–Solver results indicates that DE–Solver finds the same value as HS–Solver for P1. Although the result from HS–Solver is better for P2, a vice versa result is obtained for P3. Note that Table 1 also includes a comparison in terms of the number of required function evaluations. It can be seen that DE requires fewer function evaluations than HS. For the hybrid approach results, use of Solver as a local optimizer significantly reduces the number of function evaluations especially for P2 and P3. These function evaluation numbers are slightly higher than those obtained by using the HS–Solver approach.

. | Piezometer . | Curve matching^{a}
. | HS^{b}
. | DE . | HS–Solver^{b}
. | DE–Solver . |
---|---|---|---|---|---|---|

P1 | 119.23 | 127.87 | 130.27 | 128.74 | 128.66 | |

P2 | 130.46 | 129.6 | 127.75 | 128.74 | 128.96 | |

P3 | 129.6 | 129.6 | 130.27 | 131.33 | 130.93 | |

P1 | 8.7 | 5.69 | 5.1 | 5.46 | 5.46 | |

P2 | 4.7 | 4.81 | 5.2 | 4.92 | 4.91 | |

P3 | 5.2 | 5.1 | 5.1 | 4.96 | 5 | |

P1 | 0.5497^{c} | 0.0529 | 0.052 | 0.0519 | 0.0519 | |

P2 | 0.0306^{c} | 0.0271 | 0.0309 | 0.0264 | 0.0266 | |

P3 | 0.0056^{c} | 0.0039 | 0.0037 | 0.0035 | 0.0032 | |

Number of function evaluations | P1 | N/A | 1,494 | 816 | 120 (392^{d}) | 127 (470^{d}) |

P2 | N/A | 6,378 | 4,943 | 130 (636^{d}) | 167 (516^{d}) | |

P3 | N/A | 9,196 | 8,911 | 1,858 (658^{d}) | 2,014 (738^{d}) |

. | Piezometer . | Curve matching^{a}
. | HS^{b}
. | DE . | HS–Solver^{b}
. | DE–Solver . |
---|---|---|---|---|---|---|

P1 | 119.23 | 127.87 | 130.27 | 128.74 | 128.66 | |

P2 | 130.46 | 129.6 | 127.75 | 128.74 | 128.96 | |

P3 | 129.6 | 129.6 | 130.27 | 131.33 | 130.93 | |

P1 | 8.7 | 5.69 | 5.1 | 5.46 | 5.46 | |

P2 | 4.7 | 4.81 | 5.2 | 4.92 | 4.91 | |

P3 | 5.2 | 5.1 | 5.1 | 4.96 | 5 | |

P1 | 0.5497^{c} | 0.0529 | 0.052 | 0.0519 | 0.0519 | |

P2 | 0.0306^{c} | 0.0271 | 0.0309 | 0.0264 | 0.0266 | |

P3 | 0.0056^{c} | 0.0039 | 0.0037 | 0.0035 | 0.0032 | |

Number of function evaluations | P1 | N/A | 1,494 | 816 | 120 (392^{d}) | 127 (470^{d}) |

P2 | N/A | 6,378 | 4,943 | 130 (636^{d}) | 167 (516^{d}) | |

P3 | N/A | 9,196 | 8,911 | 1,858 (658^{d}) | 2,014 (738^{d}) |

^{a}Results of Kresic (1997).

^{b}Results of Ayvaz & Gurarslan (2018).

^{c}This value is calculated by writing the results of Kresic (1997) into the Theis model.

^{d}These values represent the function evaluations through Solver iterations.

### Example 2: leaky confined aquifer

A pumping test is conducted at a fully penetrating well in a leaky confined aquifer which is overlain by an aquitard and an unconfined aquifer. The value of *Q* is kept at 0.012 and the corresponding drawdowns are recorded at a piezometer, which is away from the well. By using these recorded data, the parameters of *T*, *S*, and *B* are determined by using the DE–Solver approach. For this analysis, Figure 3 compares the convergence plots with the ones obtained by HS–Solver.

As can be seen, both DE and DE–Solver start the search process from the same initial solution and proceed together in early iterations. After Solver gets the initial solution provided by DE, the calculated value is significantly improved without requiring many more iterations, just as given in HS–Solver. For these outcomes, the identified results are compared in Table 2 in terms of the aquifer parameters, values, and number of required iterations. It can be seen that the identified *T*, *S*, and *B* values of DE and DE–Solver agree well with the ones obtained by HS and HS–Solver although there are some differences from the results of the manual curve-matching approach. This difference is associated with the inaccuracy of the manual curve-matching approach such that the final value of the curve-matching approach is approximately 45 times bigger than that of DE–Solver. For these results, it can be observed that inclusion of Solver as a local search engine does not significantly improve the final values. When the number of function evaluations are compared, it can be seen that DE–Solver requires 1,526 function evaluations while the same result is obtained by using HS–Solver in 1,781 function evaluations. For the stand-alone DE and HS, the results given in Table 2 are obtained in 7,091 and 9,076 function evaluations, respectively. These results indicate that inclusion of the Solver as a local optimizer significantly reduces the number of function evaluations compared with the stand-alone heuristic approaches.

Solution approach . | . | . | . | . | Number of function evaluations . |
---|---|---|---|---|---|

Curve matching^{a} | 257.47 | 1.09 | 640.00 | 3.3678^{b} | N/A |

HS^{c} | 160.48 | 1.25 | 460.03 | 0.0754 | 9,076 |

DE | 163.75 | 1.23 | 473.50 | 0.0756 | 7,091 |

HS–Solver^{c} | 159.93 | 1.26 | 457.90 | 0.0753 | 1,053 (728^{d}) |

DE–Solver | 159.92 | 1.26 | 457.87 | 0.0753 | 1,122 (404^{d}) |

Solution approach . | . | . | . | . | Number of function evaluations . |
---|---|---|---|---|---|

Curve matching^{a} | 257.47 | 1.09 | 640.00 | 3.3678^{b} | N/A |

HS^{c} | 160.48 | 1.25 | 460.03 | 0.0754 | 9,076 |

DE | 163.75 | 1.23 | 473.50 | 0.0756 | 7,091 |

HS–Solver^{c} | 159.93 | 1.26 | 457.90 | 0.0753 | 1,053 (728^{d}) |

DE–Solver | 159.92 | 1.26 | 457.87 | 0.0753 | 1,122 (404^{d}) |

^{a}Results of Kresic (1997).

^{b}This value is calculated by writing the results of Kresic (1997) into the Theis model.

^{c}Results of Ayvaz & Gurarslan (2018).

^{d}These values represent the function evaluations through Solver iterations.

## CONCLUSIONS

In this study, the hybrid DE–Solver approach is first proposed and applied to the solution of aquifer parameter estimation problems by using the results of pumping tests. In the proposed approach, drawdown values at a given distance and monitoring times are calculated by means of the Theis and Hantush models for confined and leaky confined aquifers, respectively. These models are then used in conjunction with the DE–Solver approach to determine the given aquifer parameters. Evaluation of the proposed approach is conducted by solving two pumping test examples. The identified results indicate that use of a simulation–optimization approach outperforms the inadequacy of the manual curve-matching approach. Furthermore, use of a local search approach in a heuristic optimization approach significantly reduces the number of required iterations to find an optimum solution. Note that the proposed approach requires homogeneous and isotropic aquifer parameters since the problem is solved by means of the Theis and Hantush models. However, the proposed approach can be applied to problems with heterogeneous and anisotropic aquifer parameters by including the numerical solution of the governing equations in a radial coordinate system.

## ACKNOWLEDGEMENTS

A part of this paper was presented orally at the 13th International Conference on Hydroinformatics (HIC2018). The authors would like to thank two anonymous reviewers and the editor for their contributions to our paper.

## REFERENCES

*PhD thesis*

*Differential Evolution (DE)*