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

Unsteady groundwater flow toward a fully penetrating well in a confined aquifer is represented by the Theis model, which is expressed as follows (Theis 1935): 
formula
(1)
where is the drawdown at a distance 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: 
formula
(2)
 
formula
(3)
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): 
formula
(4)
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

Unsteady groundwater flow toward a fully penetrating well in a leaky confined aquifer without aquitard storage is given by the Hantush model, which is expressed as follows (Hantush & Jacob 1955): 
formula
(5)
where 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: 
formula
(6)
 
formula
(7)
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.

Let 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: 
formula
(8)
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: 
formula
(9)
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: 
formula
(10)
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: 
formula
(11)
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.

Figure 1

Flowchart of the proposed DE–Solver optimization approach (adapted from Mandal et al. 2013).

Figure 1

Flowchart of the proposed DE–Solver optimization approach (adapted from Mandal et al. 2013).

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.

Table 1

Comparison of the identified results for each piezometer

  Piezometer Curve matchinga HSb DE HS–Solverb 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 
 P1 0.5497c 0.0529 0.052 0.0519 0.0519 
P2 0.0306c 0.0271 0.0309 0.0264 0.0266 
P3 0.0056c 0.0039 0.0037 0.0035 0.0032 
Number of function evaluations P1 N/A 1,494 816 120 (392d127 (470d
P2 N/A 6,378 4,943 130 (636d167 (516d
P3 N/A 9,196 8,911 1,858 (658d2,014 (738d
  Piezometer Curve matchinga HSb DE HS–Solverb 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 
 P1 0.5497c 0.0529 0.052 0.0519 0.0519 
P2 0.0306c 0.0271 0.0309 0.0264 0.0266 
P3 0.0056c 0.0039 0.0037 0.0035 0.0032 
Number of function evaluations P1 N/A 1,494 816 120 (392d127 (470d
P2 N/A 6,378 4,943 130 (636d167 (516d
P3 N/A 9,196 8,911 1,858 (658d2,014 (738d

aResults of Kresic (1997).

cThis value is calculated by writing the results of Kresic (1997) into the Theis model.

dThese values represent the function evaluations through Solver iterations.

Figure 2

Convergence plots for DE–Solver and HS–Solver approaches for (a) P1; (b) P2; (c) P3.

Figure 2

Convergence plots for DE–Solver and HS–Solver approaches for (a) P1; (b) P2; (c) P3.

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.

Figure 3

Convergence plots of the HS–Solver and DE–Solver approaches.

Figure 3

Convergence plots of the HS–Solver and DE–Solver approaches.

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.

Table 2

Comparison of the identified , , and values and final values and number of required iterations

Solution approach     Number of function evaluations 
Curve matchinga 257.47 1.09 640.00 3.3678b N/A 
HSc 160.48 1.25 460.03 0.0754 9,076 
DE 163.75 1.23 473.50 0.0756 7,091 
HS–Solverc 159.93 1.26 457.90 0.0753 1,053 (728d
DE–Solver 159.92 1.26 457.87 0.0753 1,122 (404d
Solution approach     Number of function evaluations 
Curve matchinga 257.47 1.09 640.00 3.3678b N/A 
HSc 160.48 1.25 460.03 0.0754 9,076 
DE 163.75 1.23 473.50 0.0756 7,091 
HS–Solverc 159.93 1.26 457.90 0.0753 1,053 (728d
DE–Solver 159.92 1.26 457.87 0.0753 1,122 (404d

aResults of Kresic (1997).

bThis value is calculated by writing the results of Kresic (1997) into the Theis model.

dThese 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

REFERENCES
Ayvaz
M. T.
&
Gurarslan
G.
2017
A new partitioning approach for nonlinear Muskingum flood routing models with lateral flow contribution
.
J. Hydrol.
553
,
142
159
.
doi:10.1016/j.jhydrol.2017.07.050
.
Ayvaz
M. T.
,
Gurarslan
G.
2018
Identification of the aquifer parameters from pumping test data by using a hybrid optimization approach
. In:
HIC 2018. 13th International Conference on Hydroinformatics
(
La Loggia
G.
,
Freni
G.
,
Puleo
V.
&
De Marchis
M.
, eds),
EasyChair
, pp.
147
154
.
Ayvaz
M. T.
,
Kayhan
A. H.
,
Ceylan
H.
&
Gurarslan
G.
2009
Hybridizing the harmony search algorithm with a spreadsheet ‘Solver’ for solving continuous engineering optimization problems
.
Eng. Optim.
41
,
1119
1144
.
doi:10.1080/03052150902926835
.
Gurarslan
G.
&
Karahan
H.
2015
Solving inverse problems of groundwater-pollution-source identification using a differential evolution algorithm
.
Hydrogeol. J.
23
(
6
),
1109
1119
.
doi:10.1007/s10040-015-1256-z
.
Hantush
M. S.
&
Jacob
C. E.
1955
Non-steady radial flow in an infinite leaky aquifer
.
Eos Trans. AGU
36
,
95
100
.
Huang
Y.-C.
,
Yeh
H.-D.
&
Lin
Y.-C.
2008
A computer method based on simulated annealing to identify aquifer parameters using pumping-test data
.
Int. J. Numer. Analyt. Methods Geomech.
32
,
235
249
.
doi:10.1002/nag.623
.
Kayhan
A. H.
,
Ceylan
H.
,
Tamer Ayvaz
M.
&
Gurarslan
G.
2010
PSOLVER: a new hybrid particle swarm optimization algorithm for solving continuous optimization problems
.
Expert Syst. Appl.
37
,
6798
6808
.
doi:10.1016/j.eswa.2010.03.046
.
Kresic
N.
1997
Quantitative Solutions in Hydrogeology and Groundwater Modeling
,
1st edn
.
CRC Press
,
Boca Raton, FL, USA
.
Land
M. W. S.
1998
Evolutionary Algorithms with Local Search for Combinatorial Optimization. PhD thesis
,
University of California
,
San Diego, CA, USA
.
Mandal
D.
,
Chatterjee
A.
&
Bhattacharjee
A. K.
2013
Design of fully digital controlled shaped beam synthesis using differential evolution algorithm
.
International Journal of Antennas and Propagation
2013
(
4
),
713680
.
doi:10.1155/2013/713680
.
Mezura-Montes
E.
,
Velazquez-Reyes
J.
&
Coello Coello
C. A.
2006
A comparative study of differential evolution variants for global optimization
. In:
Proceedings of the Genetic and Evolutionary Computation Conference (GECCO'06)
,
ACM, New York, USA
.
Michalewicz
Z.
1992
Genetic Algorithms + Data Structures = Evolution Programs
.
Springer-Verlag
,
New York, USA
.
Price
K. V.
,
Storn
R. M.
&
Lampinen
J. A.
2005
The differential evolution algorithm
. In:
Differential Evolution: A Practical Approach to Global Optimization
,
Springer
,
Berlin, Germany
, pp.
37
134
.
Samuel
M. P.
&
Jha
M. K.
2003
Estimation of aquifer parameters from pumping test data by genetic algorithm optimization technique
.
J. Irrig. Drain. Eng.
129
,
348
359
.
doi:10.1061/(ASCE)0733-9437(2003)129:5(348)
.
Storn
R.
2019
Differential Evolution (DE)
.
Available from: http://www1.icsi.berkeley.edu/∼storn/code.html (accessed 5 August 2019).