A computer method based on simulated annealing to identify aquifer parameters using pumping-test data

Published online 27 April 2007 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nag.623

A computer method based on simulated annealing to identify

aquifer parameters using pumping-test data

Yen-Chen Huang, Hund-Der Yeh

and Yu-Chung Lin

SUMMARY

Conventional graphical or computer methods for identifying aquifer parameters have their own inevitable limitations. This paper proposes a computer method based on a drawdown model and a heuristic approach of simulated annealing (SA) to determine the best-fit aquifer parameters of the confined and unconfined aquifer systems. The drawdown model for the confined aquifer is the Theis solution and the unconfined aquifer is the Neuman solution. The estimated results of proposed method have better accuracy than those of the graphical methods and agree well with those of the computer methods based on the extended Kalman filter and Newton’s method. Finally, the sensitivity analyses for the control parameters of SA indicate that the proposed method is very robust and stable in parameter identification procedures. Copyrightq 2007 John Wiley & Sons, Ltd.

Received 22 November 2006; Revised 16 March 2007; Accepted 18 March 2007

KEY WORDS: groundwater; parameter identification; pumping test; simulated annealing; confined aquifer; unconfined aquifer

1. INTRODUCTION

In the past, the analysis of pumping-test data for a confined aquifer was usually made by a data-plotting and curve-fitting procedure with a type curve generated from the Theis solution[1]. The aquifer parameters, e.g. transmissivity T and storage coefficient S, are then calculated based on reading values of the match point on the graph. However, errors might be introduced during these graphical procedures.

Theis [1] obtained the solution for unsteady groundwater flow towards a well in a confined aquifer through analogy to the problem in heat conduction. Cooper and Jacob[2] developed an

∗_{Correspondence to: Hund-Der Yeh, Institute of Environmental Engineering, National Chiao Tung University,} No. 75, Bo-Ai St., Hsinchu 300, Taiwan.

†_{E-mail: hdyeh@mail.nctu.edu.tw}

approximation for the Theis equation and a data analysis method which does not require type-curve matching. The analysis of pumping-test data in an unconfined aquifer was pioneered by Boulton [3, 4], based on the theory of representing pore water drainage by an artificial index of delayed yield. Prickett[5] later developed a graphical method to determine the parameters in unconfined aquifers by employing Boulton-type curves. Cooley and Case[6] showed that Boulton’s equation is an exact solution where it describes a flow system with a rigid phreatic aquitard on top of the main aquifer and the unsaturated flow above the phreatic surface is neglected. Neuman[7, 8] presented a solution that considers the effects of elastic storage and anisotropy of aquifers on drawdown behaviour.

Neuman’s solution treated the unconfined aquifer as a compressible system and the phreatic surface as a moving boundary. His theory was extended to account for the effect of a partially pen-etrating pumping well or/and observation well in a homogeneous anisotropic unconfined aquifer. Neuman[9] also gave a graphical type-curve solution procedure to determine the aquifer param-eters. Moench [10] combined the Boulton and Neuman models for flow towards a well in an unconfined aquifer.

The aquifer parameters can also be obtained by computer methods, which usually determine parameters using the least-square approach by taking the derivative of the sum of squared errors between the observed and predicted drawdowns with respect to the parameters. The gradient-type methods are then employed to solve the nonlinear least-square equations to obtain the best-fit parameters. McElwee [11] proposed a least-squares fitting technique and sensitivity analysis to analyse the time-drawdown data for aquifer parameters. Saleem [12] proposed a nonlinear programming (NLP) technique, minimizing the sum of squares of the differences between observed and predicted drawdowns. Mania and Sucche[13] employed the least-squares approach to analyse parameters in unconfined aquifers, based on Boulton’s solution for large-time data. Sridharan

et al.[14] employed a sensitivity analysis technique based on Neuman’s model for the condition

of a fully penetrating well for identifying parameters in an unconfined aquifer.

Yeh [15] used the nonlinear least-squares and finite-difference Newton’s method (NLN) for identifying the parameters of the confined aquifer. Huang[16] used NLN to identify the unconfined aquifer parameters. The NLN approach has the advantage of high accuracy and quick convergence for reasonable initial guesses. However, those methods may yield divergent results if the guess parameter values are not very close to the target values. In addition, they may obtain poor results if improper increments were made when applying finite-difference formulae to approximate the derivative terms in the least-squared equations.

Recently, the Kalman filter has been successfully applied to the aquifer parameter and water table related estimations. Leng and Yeh[17] employed extended Kalman filter (EKF) to identify the aquifer parameters in confined and unconfined aquifer systems. The results indicated that the EKF can be utilized to analyse the drawdown data even with white noise or temporally correlated noise. However, the EKF approach still has the problems of no general guideline to assign the initial guesses and divergence also happens in some cases.

Heuristic methodologies such as simulated annealing (SA) have recently been developed in the field of stochastic optimization. SA was first proposed by Kirkpatrick et al. [18] as a method for solving combinatorial optimization problems. Zheng and Wang [19] applied SA and Tabu search to identify the parameter structure. Cunha and Sousa[20] minimized the cost in the water distribution network using SA. The solution set obtained from SA and NLP techniques for several medium-sized networks show that SA did provide a better solution in general, in comparison with that obtained with the NLP techniques. The major advantages of SA are its property of

using a descent strategy but allowing random ascent moves to avoid being trapped in a local optimum.

The objective of this study is to propose a computer method based on SA coupled with the Theis solution [1] for confined aquifers and the Neuman solution [8] for unconfined aquifers to automatically determine the best-fit aquifer parameters. In a confined aquifer, the Theis solution is used to predict the drawdown based on the guess parameters and then the SA determines the best-fit parameters to minimize the objective function which is the sum of squared errors between the observed and predicted drawdowns. On the other hand, SA combined with the Neuman solution [8] for an unconfined aquifer works in the same manner as that of a confined aquifer. Two sets of pumping-test data obtained from confined aquifers and one set from an unconfined aquifer are selected to examine the application of the proposed method. The estimated results from the proposed method are compared with those obtained from the graphical, EKF, and NLN methods. In addition, this study provides two sensitivity tests for the control parameters of the SA using data of the confined aquifer case to demonstrate the robustness and reliability when applying the proposed method for the identification of the aquifer parameters.

2. THE SIMULATED ANNEALING

The basic algorithm of SA is motivated by an analogy to the thermodynamics of annealing in solids, such as growing silicon in the form of highly ordered, defect-free crystals. In order to accomplish this, the material is annealed. It is first heated to a temperature that allows many molecules to move freely with respect to each other. After that, it is cooled slowly until the material freezes into a crystal, which is completely ordered, and thus the system is at the state of minimum energy. In other words, the molecules have high activity when the temperature is high and the crystalline configurations have various forms. If the temperature is cooled properly, the crystalline configuration is in the most stable state; thus, the minimum energy level may be naturally reached. Based on the annealing concept, SA was constructed for solving the optimization problems. During the optimization procedure, the solution, which may not be the best one, is accepted to avoid the solution becoming trapped in a local optimum.

The probability distribution of system energy at a given temperature is determined by the Boltzman probability[21]:

P(E) ∝ exp(−E/(k × Te)) (1)

where E is the system energy, k is Boltzman’s constant, Te is the temperature, and P(E) is the occurrence probability. From Equation (1), it is possible that the system might have high energy even at low temperature. Hence, the statistical distribution of energies permits the energy level of the system to escape from a local minimum. That is why the solution may not be trapped in the local optimal solution. The Boltzman probability is applied in Metropolis’ criterion [18] which takes the place
E, the difference between the objective function values of the current optimal solution and the trial solution.

As an iterative improvement method, the system starts from an initial state and is perturbed at random to a new state in the neighbourhood, for which a change of
E in the objective function

Figure 1. The flow chart of SA[21].

x is given as

x= x + (2∗D1− 1) × VM (2)

where D1 is a random number between zero and one from a uniform distribution and VM is

the step length vector. The VM can automatically be adjusted so that approximately half of all evaluations are accepted. In the minimization problem, if f(x) is smaller than f (x), then the

current solution is replaced with the trial solution. If f(x) is larger than f (x), the Metropolis criterion is then tested and a new random number D is generated between zero and one. To solve the minimization problem, the Metropolis criterion is given as[22]

PSA{accept j} = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if f( j)
f (i) exp  f(i) − f ( j)
Te  if f( j)> f (i) (3)

where PSAis the acceptance probability of the trial solution, f(i) and f ( j) are the function values

when x= xiand x= xj, and xiand xjare the current best solution and neighbourhood trial solution

of x. Generally, the control parameter Te is the current temperature and
herein is a constant, usually taken as 1, that relates temperature to the objective function. If the random number D is smaller than PSA, the current solution would be replaced by the trial solution. Otherwise, SA

would keep on generating the trial solution within the neighbourhood of the current solution. Figure 1 displays the steps in the SA algorithm. In the first step, SA initializes the solution and sets it equal to the current optimal one. The second step is to update the current optimal solution by comparing it with the generated trial solutions within a specified boundary. If a trial solution is better than the current optimal solution or if the trial solution satisfies the Metropolis criterion, the current solution is equated to the new one, otherwise, SA continues generating other trial solutions. The temperature will be decreased by multiplying a temperature reduction factor RT when there is

no improvement to the optimum after a specified number of iterations n1are performed. Based on

Equation (3), the acceptance probability becomes small with low temperature Te. The temperature should be cooled properly to guarantee that the obtained solution is the global optimum [19]. The algorithm will be terminated when SA obtains the optimal solution or the obtained solution satisfies the stopping criteria. In general, the stopping criteria are defined initially to check if the temperature is cool at the appropriate level and then to check if the difference between the optimal objective function values and those obtained in the current iteration has reached the specified value.

3. INTEGRATING SIMULATED ANNEALING WITH OTHER SOLUTIONS

This section illustrates how SA is coupled with the Theis and Neuman solutions to, respectively, identify the confined and unconfined aquifer parameters.

3.1. Theis solution

The Theis solution describing the drawdown cone within a confined aquifer in response to the pumping as a function of radial distance and time is[1]

s= q 4TW(u) (4) and u=r 2_S 4T t (5)

where s is the drawdown, q is the pumping rate, T is the transmissivity, W(u) is the well function,

r is the distance between the pumping and the observation wells, S is the storage coefficient, and t is the time since pumping started. The well function may be expressed as

W(u) =  −0.5772157 − ln u − ∞ n=1(−1) k un n· n! (6)

The higher-order terms of u in Equation (6) may be truncated when un/(n · n!) is less than 10−7.

3.2. Neuman’s solution

Heuristic methods coupled with Neuman’s solution [8] are applied to identify the parameters of the unconfined aquifer. Given the initial guesses of the parameters, the four best-fit parameters are identified when the convergence criteria are met. The solution which describes the groundwater flow system in an unconfined aquifer developed by Neuman[8] is

s(r, z, t) = q 4T _∞ 0 4y J0(y1/2)  u0(y) + ∞  n=1un(y) dy (7)

where J0(x) is the zero-order Bessel function of the first kind,  = Kzr2/Krb2is a dimensionless

parameter, Kr and Kz are, respectively, horizontal and vertical hydraulic conductivities, b is the

thickness of aquifer, y is a dummy variable, and

u0(y) = {1 − exp[−ts(y

2_{− }2 0)]} cosh(0zD) [y2_{+ (1 + )}2 0− (y2− 20)2/] cosh(0) ×sinh[0(1 − dD)] − sinh[0(1 − lD)] (lD− dD) sinh(0) (8)

un(y) = {1 − exp[−ts(y

2_{+ }2 n)]} cos(nzD) [y2_{− (1 + )}2 n− (y2+ 2n)2/] cos(n) ×sin[n(1 − dD)] − sin[n(1 − lD)] (lD− dD) sin(n) (9)

where ts= T t/Sr2 represents the dimensionless time since pumping started, S equals Ss · b, zD= z/b is the dimensionless elevation of the observation point,  = S/Sy is a dimensionless

parameter, dD= d/b denotes the dimensionless vertical distance between the top of the perforation

in the pumping well and the initial position of the water table, and lD= l/b is the dimensionless

vertical distance between the bottom of the perforation in the pumping well and the initial position of the water table. The terms₀and_nare, respectively, the roots of the following two equations:

r0sinh(0) − (y2− 02) cosh(0) = 0, 20<y2 (10)

and

Figure 2. The flowchart of the identification procedure: (a) conventional method and (b) present method.

3.3. The objective function and error criteria

The aquifer parameters can be estimated based on the Theis solution for confined aquifers and the Neuman solution for unconfined aquifers when minimizing the sum of squared errors between the observed and predicted drawdowns. Therefore, the objective function used to replace the energy defined in Equation (1) and to be minimized is defined as

f(x) = n



i=1(Ohi− Phi)

2 ₍₁₂₎

where Ohi and Phi are, respectively, the observed and predicted drawdowns at different time steps

and n is the total number of time steps.

The SA searches for the optimal parameters depending on the objective function value. The initial guesses for SA are provided by the user; however, the SA algorithm allows the initial

guesses to be randomly given. After the initial guesses are made, the predicted drawdown can be calculated from Equation (4) or (7). Then all the possible solutions (trial solutions) will be kept and improved based on the objective function value. If the objective function value meets the specified stopping criterion, the SA process will be terminated and the optimal parameters are found. The aquifer parameter identification procedures using the conventional approach and SA are illustrated in Figure 2(a) and (b), respectively.

To assess the accuracy of the estimated parameters, two error criteria, mean error (ME) and standard error of the estimate (SEE)[15], are used to calculate the errors between the observed and predicted drawdowns in this study. The ME is defined as the average of the sum of errors between the observed and predicted drawdowns. When the ME value is equal to or very close to zero, the assumption that errors have zero mean will be satisfied. The SEE is defined as square root of the sum of squared errors between the observed and predicted drawdowns divided by the degree of freedom, which equals the number of observed data points minus the number of unknowns.

4. RESULTS AND DISCUSSION

Two sets of pumping-test data obtained from the confined aquifer are chosen to examine the application of the proposed method. The first pumping test performed in a confined aquifer with a fully penetrating well is taken from Todd [23]. The well was pumped at a uniform rate of 2500 (m3/day) and the drawdown is measured at the observation well 60 m away from the pumping well. The second data set is taken from Walton[24]. The test was conducted on 2 July 1953 at Gridley, Illinois. The constant pumping rate was 1200 (m3/day) for about 8 h and the observation well (well 1) was 251 m away from the pumping well (well 3).

The pumping test for an unconfined aquifer was done in 1965 by the Bureau of Geologic Research and Minerals of France[9, 25]. The unconfined aquifer consists of medium-grained sand with gravel in the deeper part and a clayey matrix at shallow depths. The initial saturated thickness of the aquifer equals 8.24 m. The discharge rate averages about 53 m3/h and the drawdowns are monitored at a distance of 10 m from the pumping well.

4.1. Identification of confined aquifer parameters

The upper and lower bounds of T are, respectively, 3000 and 0(m2/day) and for S are, respectively, 10−3 and 10−5 for a confined aquifer when applying the present method. The initial temperature of SA is chosen to be 10◦ and the temperature is decreased by a reduction factor of 0.75 after 100 calculations. The choice of the initial temperature is generally case by case. Nevertheless, Kirkpatrick et al.[18] gave a guideline for the initial temperature that the acceptance probability happened at the lower part of Equation (3), i.e. when the trial solution is worser than the current solution, should be larger than 80% initially. This criterion has the merit of avoiding the situation that the current solution is trapped in a local optimum at early search. The parameter identification process of SA is terminated if the absolute differences between two successive objective function values are all less than 10−6 within four iterations or the number of evaluations is greater than 106 times.

The estimated results and related errors using Todd’s data are shown in Table I. The estimated T and S obtained from SA are 1138(m2/day) and 1.93 × 10−4, respectively. The ME is−0.54 × 10−4 and the SEE value is 5.46 × 10−3. The analysed results along with their prediction errors obtained

Table I. Comparison of results when applying graphical methods, NLN, EKF, and SA, to analyse Todd’s pumping-test data obtained from a confined aquifer.

Estimated parameters Prediction errors

Method T (m2/day) S× 10−4 ME× 10−4(m) SEE1× 10−3 (m) Graphical methods[15] Theis 1110 2.06 −17.20 8.85 Cooper–Jacob 1090 1.84 −309.60 35.22 Computer methods NLN 1139 1.93 3.97 5.47 EKF 1140 1.92 −0.27 5.47 SA 1138 1.93 −0.54 5.46

from conventional graphical methods such as the Theis and Cooper–Jacob mentioned in Todd[23] are also listed in Table I. It can be seen that the estimated aquifer parameters and prediction errors obtained from SA have much better accuracy than those of the graphical methods. The estimated

T, S, and predictions errors by the NLN and EKF methods are also shown in Table I, indicating that

these three computer methods yield results with the same degree of accuracy. Figure 3 describes the estimated drawdown obtained by the present method and the pumping-test data for the observation well in the confined aquifer. It is clear that the estimated drawdowns fit the pumping-test data quite well.

Similarly, the estimated results and related prediction errors based on Walton’s data by the graphical and computer methods are shown in Table II. The estimated T and S obtained from the Theis method for test drawdown data measured at well 1 are 125.4 (m2/day) and 2.05 × 10−5, respectively, and the estimated T from Cooper–Jacob method at the pumping well is 150.3 (m2/day). Notice that the value of S was not determined because the effective radius of the pumped well was not known[24]. The T and S estimated by SA are 124.2 (m2/day) and 2.05 × 10−5, respectively, and the ME and SEE values are −1.33 × 10−3 and 3.72 × 10−2, respectively. In this case, the computer methods still provide better results than that of Theis’ method judged from the predic-tion errors. Figure 4 displays the pumping-test data and the estimated drawdowns based on the parameters obtained by the SA. Clearly, the estimated drawdown data fit the pumping-test data quite well as indicated in the figure.

4.2. Identification of unconfined aquifer parameters

The upper and lower bounds of Kr are, respectively, 10−2and 10−4(m/s), of Kzare, respectively,

10−3 and 10−5 (m/s), of S are, respectively, 5 × 10−3 and 10−5, and of Sy are, respectively,

3× 10−1and 10−2for an unconfined aquifer when applying SA. The other control parameters of SA are similar to those of the confined aquifer case given in the previous section.

The estimated results and related errors resulting from the computer methods such as NLN, EKF [17], and the present methods are listed in Table III. The estimated parameters obtained by the SA are: 2.23 × 10−3(m/s) for Kr; 1.67 × 10−5(m/s) for Kz; 1.31 × 10−3for S; and 3.83 × 10−2for Sy, respectively. In addition, Table III also lists the analysed results and the prediction errors from

the graphical methods such as the Neuman type-curve method and Neuman’s semilogarithmic method [25]. The prediction errors by using SA are generally much smaller than those by two

Figure 3. The estimated drawdown and Todd’s pumping-test data[23] obtained from a confined aquifer using SA.

Table II. Comparison of results when applying graphical methods, NLN, EKF, and SA methods to analyse Walton’s pumping-test data at well 1 obtained from a confined aquifer.

Estimated parameters Prediction errors

Method T (m2/day) S× 10−5 ME× 10−3(m) SEE× 10−2 (m) Graphical methods[24] Theis 125.4 2.00 −7.23 3.80 Cooper–Jacob 150.3 — — — Computer methods NLN 124.2 2.05 −1.33 3.72 EKF 124.2 2.04 −4.59 3.73 SA 124.2 2.05 −1.33 3.72

Note: The transmissivity estimated by Cooper–Jacob method was based on the test data measured at the observation well and the symbol — denotes that data are not available.

graphical methods, indicating that the estimated parameters of SA give a better fit to the observed drawdown data. Figure 5 displays the estimated drawdown and the pumping-test data in the unconfined aquifer. This figure also indicates that the proposed methods can optimally search the

Figure 4. The estimated drawdown and Walton’s pumping-test data at well 1[24] obtained from a confined aquifer using SA.

Table III. Comparison of results when applying graphical methods, NLN, EKF, and SA methods to analyse the pumping-test data obtained from an unconfined aquifer.

Estimated parameters Prediction errors

Kr× 10−3 Kz× 10−5 ME× 10−3 SEE× 10−3 Method (m/s) (m/s) S× 10−3 Sy× 10−2 (m) (m) Graphical methods Neuman-type curve 2.40 1.62 1.46 5.73 32.90 34.59 Neuman semilogarithmic 2.40 1.62 1.87 2.13 14.23 14.96 Computer methods NLN 2.22 1.68 1.31 3.85 0.28 8.06 EKF 2.25 1.56 0.97 4.10 1.68 8.36 SA 2.23 1.67 1.31 3.83 0.31 8.06

parameters of the unconfined aquifer. Clearly, these estimated results and related errors demonstrate that the proposed methods are much superior to the graphical methods and give the results with the same degree of accuracy when compared with those of NLN and EKF.

Figure 5. The estimated drawdown and the pumping-test data obtained from an unconfined aquifer using SA.

5. THE SENSITIVITY ANALYSIS OF SA’S CONTROL PARAMETERS

The use of control parameters in SA, such as the initial guess value and the temperature reduction factor, may affect the results of the parameter identification. For demonstrating robustness and reliability of SA in parameter identification, this study presents two sensitivity analyses of the control parameters in SA for the parameter identification when analysing Todd’s pumping-test data[23].

Deciding on the initial guess value and the temperature reduction factor, RT, may be difficult

in the application of SA. This section, therefore, focuses on the sensitivity analysis of these two factors. Table IV lists various values of initial guesses and the analysed results of using NLN and EKF in nine case studies[17]. Unfortunately, in four out of nine cases, the NLN and EKF fail to get convergent results. However, Table IV lists the values of initial guesses and the convergent results for those nine cases with different initial guesses when using SA in identifying confined aquifer’s parameters. In those nine cases, the estimated T ranges from 1138 to 1140(m2/day) and S ranges from 1.92 × 10−4to 1.93 × 10−4. When compared to NLN and EKF, Table IV indicates that SA allows a wider range of initial guesses and gives the same order of magnitude for the prediction error SEE. Table IV indicates that SA is successfully applied to identify aquifer parameters even when the field geological information is not available and the initial guesses are hard to decide.

Table IV. Comparison of results when applying NLN, EKF, and SA.

Initial guesses Convergence or not Estimated results Prediction errors T (m2/day) S NLN EKF SA T (m2/day) S×10−4 ME×10−4(m) SEE×10−3 (m)

700 10−3 No Yes Yes 1140 1.92 −0.27 5.47

700 10−4 Yes Yes Yes 1140 1.92 −0.27 5.47

700 10−5 Yes No Yes 1140 1.92 −0.27 5.47

1300 10−3 No Yes Yes 1140 1.92 −0.27 5.47

1300 10−4 Yes Yes Yes 1139 1.93 3.97 5.47

1300 10−5 Yes No Yes 1140 1.92 3.97 5.47

2000 10−3 No No Yes 1140 1.92 −0.27 5.47

2000 10−4 No Yes Yes 1138 1.93 −0.54 5.46

2000 10−5 Yes No Yes 1138 1.93 −0.54 5.46

Table V. The results of parameter identification using various RT values.

RT T (m2/day) S×10−4 CPU time (s)

0.50 1138.2 1.93 5.00 × 10−2 0.55 1138.2 1.93 6.62 × 10−2 0.60 1138.2 1.93 6.60 × 10−2 0.65 1138.1 1.93 6.20 × 10−2 0.70 1138.2 1.93 9.72 × 10−2 0.75 1138.2 1.93 1.12 × 10−1 0.80 1138.1 1.93 1.28 × 10−1 0.85 1138.1 1.93 2.06 × 10−1 0.90 1138.2 1.93 2.69 × 10−1 Statistical summary Mean 1138.2 1.93 × 10−4 Standard deviation 3.06 × 10−3 3.60 × 10−9

Table V shows the estimated parameters and their mean and standard deviation when using the

RT value varying from 0.50 to 0.90 with 0.05 increments. The identified parameter T ranges from

1138.1 to 1138.2 (m2/day) and the parameter S is 1.93 × 10−4for all cases. The mean of T and

S is 1138.2 (m2/day) and 1.93 × 10−4, respectively, which is very close to those estimated by various methods as shown in Table I. The standard deviations of T and S are 3.06 × 10−3 and 3.60 × 10−9, respectively, which are very small when compared to their mean values, indicating that the identified results are independent of RT values. In other words, the influence of choosing

various values of RT on the results of the parameter identification is negligible. In addition, Table

V also shows that the CPU time increases with the value of RT because the cooling schedule slows

down when RT gets larger. The results of the sensitivity analysis indicate that one can apply the

6. SUMMARY AND CONCLUSIONS

A novel method is developed based on simulated annealing (SA) integrated with aquifer draw-down equations to identify aquifer parameters of confined and unconfined aquifer systems. The Theis solution combined with SA can optimally determine the aquifer transmissivity and storage coefficient for a confined aquifer if the assumptions of Theis model are satisfied. Likewise, the Neuman solution can also be employed with SA to estimate the horizontal and vertical hydraulic conductivities, storage coefficient, and specific yield for an unconfined aquifer if the assumptions of Neuman model are met. Two sets of pumping-test data in the confined aquifers and one set in the unconfined aquifer are utilized to demonstrate the application of the proposed method in parameter identification. The results show that the present method can determine the aquifer parameters with very good accuracy. The identified results and related prediction errors indicate that the proposed method is superior to the graphical methods and gives results with the same degree of accuracy when compared with those of NLN and EKF.

The analysed results based on SA with various control parameters are compared and discussed. In the confined case, SA is shown to have a wide range of initial guess values and the accuracy of estimated results is as good as those of other computer methods such as EKF and NLN. In addition, the temperature reduction factor does not seem to affect the results of the parameter identification. These analyses demonstrate that the proposed method is robust and reliable even if the user is not experienced in using SA.

ACKNOWLEDGEMENTS

This study was partly supported by the Taiwan National Science Council under the grant NSC92-2211-E-009-008. The authors sincerely thank three anonymous reviewers for constructive comments and suggested revisions.

REFERENCES

1. Theis CV. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground-water storage. EOS Transactions of the AGU 1935; 16:519–524.

2. Cooper Jr HH, Jacob CE. A generalized graphical method for evaluating formation constants and summarizing well field history. Transactions, American Geophysical Union 1946; 27(IV):526–534.

3. Boulton NS. Unsteady radial flow to a pumped well allowing for delayed yield from storage. International

Association of Scientific Hydrology Publications 1954; 37:472–477.

4. Boulton NS. Analysis of data from non-equilibrium pumping tests allowing for delayed yield from storage.

Proceedings of the Institution of Civil Engineers, vol. 26, 1963; 469–482.

5. Prickett TA. Type-curve solution to aquifer tests under water-table conditions. Ground Water 1965; 3:5–14. 6. Cooley RL, Case CM. Effect of a water-table aquitard on drawdown in an underlying pumped aquifer. Water

Resources Research 1973; 9(2):434–447.

7. Neuman SP. Theory of flow in unconfined aquifers considering delayed response of the water table. Water

Resources Research 1972; 8:1031–1044.

8. Neuman SP. Effects of partial penetration on flow in unconfined aquifers considering delayed aquifer response.

Water Resources Research 1974; 10:303–312.

9. Neuman SP. Analysis of pumping test data from anisotropic unconfined aquifers considering delayed gravity response. Water Resources Research 1975; 11:329–342.

10. Moench AF. Combining the Neuman and Boulton models for flow to a well in an unconfined aquifer. Ground

Water 1995; 33:378–384.

11. McElwee CD. Theis parameter evaluation from pump tests by sensitivity analysis. Ground Water 1980; 18:56–60. 12. Saleem ZA. A computer method for pumping test analysis. Ground Water 1970; 13:21–24.

13. Mania J, Sucche M. Automatic analysis of pumping test data—application to the Boulton and Hantush hypothesis.

Journal of Hydrology 1978; 37:185–194.

14. Sridharan K, Ramaswamy R, Lakshmana Rao NS. Identification of parameters in unconfined aquifers. Journal

of Hydrology 1985; 79:73–81.

15. Yeh HD. Theis’ solution by nonlinear least-squares and finite-difference Newton’s method. Ground Water 1987; 25:710–715.

16. Huang WC. A numerical method for determining parameters of unconfined aquifer. MS Thesis, National Chiao Tung University, Hsinchu, Taiwan, 1996.

17. Leng CH, Yeh HD. Aquifer parameter identification using the extended Kalman filter. Water Resources Research 2003; 39(3):1062. DOI: 10.1029/2001WR000840

18. Kirkpatrick S, Gelatt Jr CD, Vecchi MP. Optimization by simulated annealing. Science 1983; 220(4598):671–680. 19. Zheng C, Wang P. Parameter structure identification using Tabu search and simulated annealing. Advances in

Water Resources 1996; 19(4):215–224.

20. Cunha MDC, Sousa J. Water distribution network design optimization: simulated annealing approach. Journal of

Water Resources Planning and Management (ASCE) 1999; 125(4):215–221.

21. Pham DT, Karaboga D. Intelligent Optimisation Techniques. Springer: Great Britain, 2000.

22. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. Journal of Chemical Physics 1953; 21(6):1087–1092.

23. Todd DK. Ground Water Hydrology (2nd edn). Wiley: New York, 1980. 24. Walton WC. Groundwater Resource Evaluation. McGraw-Hill: New York, 1970.