Published online 23 January 2008 in Wiley InterScience (www.interscience.wiley.com) DOI: 10.1002/hyp.6875
Identifying groundwater pumping source information using
simulated annealing
Yu-Chung Lin and Hund-Der Yeh*
Institute of Environmental Engineering, National Chiao Tung University, No.75, Po-Ai Street, Hsinchu City 300, Taiwan
Abstract:
This study develops an approach referred to as SA-MF using simulated annealing (SA) and a three-dimensional groundwater flow model (MODFLOW-2000) to determine the pumping source information, namely the pumping source location, the pumping rate and the pumping period. Eight scenarios, each with various cases for aquifers with homogeneous or heterogeneous conductivities, are used to test the applicability and accuracy of SA-MF in determining the pumping source information. The results show that at least five measured heads should be used to analyse the pumping source identification problem and at least three observation wells are required to effectively determine the pumping source information. The SA-MF gives good estimated results in both synthetic and real problems even if the measured heads contain measurement errors. Copyright 2008 John Wiley & Sons, Ltd.
KEY WORDS groundwater; inverse problem; heuristic approach; pumping source identification; simulated annealing Received 13 November 2006; Accepted 18 July 2007
INTRODUCTION
Groundwater is an ideal freshwater resource for human usage. The merits of using groundwater include con-venient availability near the point of use, excellent quality and relatively low cost of development (Todd and Mays, 2005). The groundwater resource reveals its important significance as there is a shortage of sur-face freshwater for domestic and industrial uses as well as for irrigating agriculture which is the main con-sumer of groundwater. Effective and appropriate man-agement strategy of groundwater usage can provide fresh groundwater for human uses and environmental functions (Lo´aiciga, 2004). In Taiwan, the groundwater resource is not allowed to be pumped for private use without obtain-ing a governmental permit as the groundwater resource is a public property. Illegal pumping may have adverse effects on the groundwater flow system. Effects may include seawater intrusion into coastal aquifers, land subsidence, upwelling of poor-quality groundwater into fresh-water aquifers, decreased base flow and spring flow to support riparian and aquatic ecosystems, and aquifer dewatering and increased pumping cost (Lo´aiciga, 2004). Moreover, illegal pumping may also disturb the quality of the groundwater in a groundwater contamination site, change the groundwater flow pattern or agitate the capture zone at an ongoing remediation site. There is therefore a need to develop a methodology for the determination of unknown pumping source information from a practical viewpoint.
* Correspondence to: Hund-Der Yeh, Institute of Environmental Engi-neering, National Chiao Tung University, No.75, Po-Ai Street, Hsinchu City 300, Taiwan.
E-mail: hdyeh@mail.nctu.edu.tw
Bear (1979) indicated that the groundwater manage-ment problem could be solved by incorporating a simu-lation model with an optimization program. The ground-water management problems could be classified into two categories (Todd and Mays, 2005). The first is the hydraulic management problems that are aimed at man-aging pumping and recharge. The second is the policy evaluation problems that consider the economics of water allocations. Traditionally, linear programming, nonlin-ear programming or a Newton-type approach was often employed to solve the optimization problem. It is diffi-cult to find optimal solutions for large-scale combinatorial optimization problems using gradient-type approaches or to resolve parameter estimation problems with complex search spaces (Cai et al., 2001). In general, two major difficulties obstruct inexperienced users. The first is that the gradient-type approaches might converge to a local solution nearest to the starting point if an improper ini-tial guess was made (Cai et al., 2001). In addition, the gradient-type approaches need to compute the deriva-tives with respect to the decision variables and iteratively update the solution to obtain the optimal solution. The second is that those derivatives may be difficult to cal-culate analytically or numerically in the highly nonlinear and non-convex optimization problems (Shieh and Per-alta, 2005).
Stochastic optimization methods such as genetic algo-rithms (GAs) and simulated annealing (SA) are global optimization solvers and have been widely used recently. Differing from the gradient-type approaches, these global optimization methods iteratively renew the trial solu-tion by the objective funcsolu-tion value (OFV) in determin-ing the nearby global optimum. The first advantage of using global optimization methods is that the users don’t
require much experience in providing the initial guesses for solving nonlinear problems. The initial guesses can be arbitrarily given or even generated by a random number generator, and still achieve optimal results. The second advantage is that these global optimization methods do not need to calculate the derivates to renew the trial solu-tion.
SA is a random search algorithm that allows, at least in theory or in probability, to obtain the global optimum of a function in a given domain. The philosophy of SA is to imitate the physical annealing process of heating up a solid and then cooling it down slowly until it crystallizes. The temperature should be reduced properly to minimize the energy state and crystallize a regular structure. Comparatively, if the cooling process occurs too fast, the crystallized structure becomes irregular and its energy does not stay in minimum. Romeo and Sangiovanni-Vincentelli (1991) used homogeneous and inhomogeneous Markov chain theories to prove that SA can converge to global optimal solutions.
SA evolved from the descent search method, which has the drawback of obtaining the solution in a local optimum (Rayward-Smith et al., 1996). Basically, SA starts with a single initial guess named as the current optimum. Next, a trial solution neighboured to the current optimum is randomly generated and its corresponding OFV is then calculated. Finally, the nearby global optimum is obtained by improving the current optimum iteratively based on the OFV. To avoid obtaining a local optimum, SA uses the Metropolis mechanism to control which ascent moves, i.e. inferior solution, could be accepted as a current optimum (Rayward-Smith et al., 1996). The Metropolis mechanism is a criterion based on Boltzman’s probability for the difference between the OFVs of current optimum and trial solution at a given state (Pham and Karaboga, 2000). Therefore, the SA has a property of going descent but allowing random ascent to avoid a possible trap in a local optimum, which prevents the SA from having the same problem as the descent method.
SA was successfully applied in a wide range of opti-mization applications such as capacity extension for pipe network system (Cunha and Sousa, 1999; Monem and Namdarian, 2005), parameter calibration and identifica-tion problems (Zheng and Wang, 1996; Cooper et al., 1997; Li et al., 1999; Gou and Zheng, 2005; Lin and Yeh, 2005; Yeh et al., 2007), groundwater management problems (Doutherty and Marryott, 1991; Marryott et al., 1993), groundwater remediation system problems (Kuo et al., 1992; Marryott, 1996; Rizzo and Dougherty, 1996; Tung et al., 2003; Shieh and Peralta, 2005) and source identification of groundwater contamination problems (Aral and Guan, 1996; Aral et al., 2001; mar and Sayeed, 2005, 2006; Sayeed and Mahinthaku-mar, 2005). Shieh and Peralta (2005) mentioned that the straightforward formulation, no requirement for comput-ing derivatives and easy implementation with ground-water simulation models, are the main advantages of employing SA to solve the optimization problem.
The aim of this work is to propose an approach referred to as SA-MF based on simulated annealing (SA) and a three-dimensional groundwater flow model (MODFLOW-2000) to determine the pumping source information, namely the pumping source location, pump-ing rate and pumppump-ing period, in a groundwater flow sys-tem. Three synthetic problems are designed to test the performance of the proposed approach. The problems include a hypothetical pumping well which is assumed at a specific location and withdraws the groundwater at a specified flow rate for a period of time. The tempo-ral distribution of water level (referred to as measured hydraulic head or simply measured head) at the obser-vation wells can be either observed at the field for real problems or simulated by MODFLOW-2000 for testing the performance of SA-MF. In the estimation of ing source information, the trial solutions for the pump-ing source location, pumppump-ing rate and pumppump-ing time are generated by SA. The trial solutions are used as input to MODFLOW-2000 to generate the simulated heads. The pumping source information can then be determined by SA-MF when minimizing the sum of square errors between the measured heads and simulated heads at the observation wells.
This study is organized into six sections. Following this introduction, brief illustrations for the groundwater flow model, SA and a detailed description for SA-MF are given in the second section. Sections 3– 5 contain the problems of a homogeneous aquifer with four scenarios, heterogeneous aquifers with three scenarios and a real aquifer with one scenario to test the performance of SA-MF and related requirement in the estimation of pumping source information. Finally, we summarize the findings and offer concluding remarks in the sixth section.
METHODOLOGY Groundwater Flow Model
The three-dimensional equation describing the ground-water flow can be expressed as (Harbaugh et al., 2000)
∂ ∂xi Kij ∂h ∂xi CW D Ss ∂h ∂t 1
where h is the hydraulic head (l), Kij is the hydraulic
conductivity tensor (l T1), S
sis the specific storage (l1),
W is the volumetric flux per unit volume representing sources and/or sink (negative for flow out and positive for flow in (l T1)), xi are the Cartesian coordinates, and
tis time (T).
Equation (1), together with the flow and/or head con-ditions at the boundaries of an aquifer system and ini-tial head condition, constitutes the mathematical model for a groundwater flow system. The computer model MODFLOW-2000 developed based on Equation (1) can be used to simulate the head distribution for a given aquifer system. MODFLOW-2000, an extended version of MODFLOW, can simulate steady state and transient
flow in the flow system with an aquifer which is con-fined, unconfined or a combination of these (Harbaugh et al., 2000).
Simulated annealing
Similar to Newton’s methods, the SA starts with an initial guess xkand calculates the objective function value
(OFV) fx which is treated as the current optimum. Then one trial solution x0
k is generated near the current
optimal solution to update the optimal solution based on the OFV, fx0
k. In minimization problems, if fx0k <
fxk the current optimal solution takes place with
trial solution x0
k. If fx0k ½ fxk, one has to use the
Metropolis criterion to determine the acceptance of the trial solution. The Metropolis criterion is given as (Pham and Karaboga, 2000): PSAD 1, fx0 k fxk exp fxk fx0k Te , fx0 k > fxk 2 where PSAis the acceptance probability of the trial
solu-tion and Te, a control parameter, is the current tempera-ture. If PSA is greater than the random number generated
from a uniform distribution and distributed between 0 and 1, then the trial solution is accepted as the new cur-rent optimal. Otherwise, the trial solution is rejected. The algorithm will be iteratively applied and terminated when the obtained solution satisfies the stopping criterion. In general, two stopping criteria are used in SA. One is to check whether the total objective function evaluations exceed the specified maximum number of iterations. This criterion is applied to avoid the waste of large computer time if the optimal solution is too far to obtain. The sec-ond is to check whether the current optimal OFV is small enough to be considered a minimum.
SA-MF implementation
Figure 1 depicts a flowchart describing the eleven steps of SA-MF. The first step is to initialize random initial solutions within the specified solution domains. The second step is to generate the simulated heads by MODFLOW-2000 with the estimated solutions. The objective function is defined as
Minimize f D 1 nwellðnstep nstep mD1 nwell lD1 hlm,simhlm,mes2 3
where nstep is the number of the measured heads, nwellis
the number of observation wells, hlm,simis the simulated
head at the mth time step at the lth observation well, hlm,mes is the measured head at the mth time step and the
lth observation well. Note that the initial solutions are considered as the current optimal solutions.
The third step is to generate a new trial solution for source location x0
k given as
x0kDxkCInteger[2 ð RD11 ð VMk] 4
and the trial solutions for pumping rate and pumping period are generated as
x0kDxkC2 ð RD21 ð VMk 5
where RD1 and RD2 are random numbers ranging from
0– 1 and produced from a uniformly distributed random generator and VMk is a step length vector. The
vari-able VMk is first defined as the length of the difference
between the upper and lower bounds and then automati-cally adjusted periodiautomati-cally so that half of the evaluations could be accepted as new optima in the next exploration domain. The goal here is to sample the function widely. The acceptance of the trial solutions below 0Ð4 indi-cates that the global optimum does not locate within the exploration domain. Comparatively, a high percentage of points accepted for xi means that the relevant element
of VMk is enlarged. If a new trial solution does not fall
within the solution domain, the SA-MF has to generate another. The fourth step is to simulate the hydraulic head distribution by MODFLOW-2000 with the trial solutions and calculate its corresponding OFV. With the OFV, the trial solution is checked to determine if whether or not this is a new optimum in the fifth step. If the OFV sat-isfies Metropolis criterion, the current optimal solution is replaced by the trial solution in the sixth step. Other-wise, the algorithm will continue generating a new trial solution in the seventh and eighth steps.
After NS steps through all considered variables, the step length vector VMk is adjusted so that 50% of all
moves are accepted in the ninth step. In other words, after N ð NS function evaluations, each element of VMk
is adjusted so that half of all function evaluations are accepted. Note that N represents the number of consid-ered variables and NS represents the number of steps at a specific temperature. After NT times through the above loop, the temperature Te is reduced in the tenth step. Accordingly, after N ð NT ð NS function evaluations, the temperature is decreased by the temperature reduc-tion factor RTe even if no improvement in the optimum
takes place. The new temperature is then
Te0DRTeðTe. 6
Note that the value of RTeis generally chosen to be less
than one (Pham and Karaboga, 2000). The temperature should be reduced properly to guarantee the obtained solution is the global optimal solution (Zheng and Wang, 1996). In the eleventh step, the identification process will be terminated when the obtained solution satisfies the stopping criterion, which is used to check whether the absolute value of the difference between the two OFVs obtained at two consecutive temperatures is less than 106, four times in succession.
APPLICATION TO A SYNTHETIC HOMOGENEOUS AQUIFER
Four unknown variables, i.e. the source location of x and y coordinates, pumping rate and pumping period,
Figure 1. Flowchart depicting SA-MF processes are considered for solving the problem of pumping
source identification. Note that SA-MF determines the source location at the centre of the finite difference grid. Four scenarios are considered in this section. Scenario 1 with eight cases is designed to verify the minimum number of measured heads for effectively solving the pumping source identification problem. The purpose of scenario 2, with 30 cases, is to examine the influence of measurement error on the estimated results. Scenario 3 with eight cases is used to find an effective way to increase the accuracy of the source identification with high measurement error level. Finally, the fourth scenario with 23 cases is to determine the required
number of observation wells to identify the pumping source.
Site description
A hypothetical site is used to illustrate the identifica-tion procedure for a pumping source in a homogenous and isotropic confined aquifer. The aquifer length and width are both 1000 m and the aquifer thickness is 20 m. The top and bottom elevations, hydraulic conductivity, porosity and storage coefficient are 10Ð0 and 10Ð0 m, 5 ð 106 m s1, 0Ð1 and 1 ð 104, respectively. The
finite difference grids are block-centred and the related boundary conditions for the aquifer system are shown in
Figure 2. Plan view of hypothetical aquifer site
Figure 2. The upper and lower boundaries are specified as no-flux and the left and right boundaries are specified as constant-head. The elevations of the hydraulic head at left and right boundaries are 30Ð0 m and 20Ð0 m, respec-tively. The grid width and length are both 10 m; thus, the number of finite difference meshes is 100 ð 100. It is assumed that the pumping well P1 is located at (495 m, 495 m) and the pumping rate and pumping period are 200Ð0 cubic metres per day (m3 day1) and 60Ð0 min,
respectively. The measured heads at 15 observation wells indicated in Figure 2 are simulated using MODFLOW-2000 and listed in Table I. Note that the drawdown at a distance of 500 m away from P1 is smaller than 0Ð001 m when the pumping rate is 200 m3 day1 and the
pump-ing period is less than 170 min based on the calculation using Theis solution (Todd and Mays, 2005).
All the nodal points are considered as suspected pumping source locations except those meshes containing the observation wells. The upper bounds for the pumping rate and pumping period are set as 1000Ð0 m3 day1 and 120Ð0 min, respectively, and the lower bounds are both zero. The NS, NT, initial temperature and RTeare chosen
as 20, 10, 5 and 0Ð8, respectively. All initial guesses in the case studies are produced by a random number generator. Scenario 1: Number of measured heads
The purpose of this scenario is to investigate the required number of measured heads in determining the pumping source information, namely the source location, and the pumping rates and periods, at a homogeneous and isotropic aquifer site. Mathematically, at least four measured drawdowns are needed to solve four unknowns. Eight cases, with different combinations of observation wells as indicated in Table II, are designed to study the influence of the number of measured drawdowns on the results of pumping source identification.
Table I. Location and observed hydraulic head for homogeneous and isotropic confined aquifers
Observation Location Measured heads(m)
wells
Initial Pumping period state T D30 (min) T D60 (min) N-1 (465 m, 495 m) 25Ð041 22Ð621 21Ð275 N-2 (395 m, 495 m) 25Ð040 24Ð769 24Ð388 NW-1 (395 m, 395 m) 26Ð031 25Ð945 25Ð781 NW-2 (295 m, 295 m) 27Ð023 27Ð02 27Ð013 W-1 (495 m, 425 m) 25Ð734 25Ð086 24Ð408 W-2 (495 m, 395 m) 26Ð031 25Ð759 25Ð378 SW-1 (595 m, 395 m) 26Ð031 25Ð945 25Ð781 SW-2 (695 m, 295 m) 27Ð022 27Ð020 27Ð012 S-1 (585 m, 495 m) 25Ð040 24Ð680 24Ð216 S-2 (595 m, 495 m) 25Ð040 24Ð769 24Ð388 SE-1 (595 m, 595 m) 24Ð051 23Ð965 23Ð801 E-1 (495 m, 535 m) 24Ð645 22Ð958 21Ð809 E-2 (495 m, 595 m) 24Ð051 23Ð779 23Ð398 E-3 (495 m, 695 m) 23Ð062 23Ð043 22Ð995 NE-1 (395 m, 595 m) 24Ð050 23Ð964 23Ð800
Four measured heads are used to determine the pump-ing source information in cases 1Ð1 to 1Ð4 and five mea-sured heads are employed in cases 1Ð5 to 1Ð8. The anal-ysed results using four observations with different alloca-tion of the observaalloca-tion wells are shown in Table II. The pumping locations are correctly determined at (495 m, 495 m) in cases 1Ð1, 1Ð2 and 1Ð4. However, the estimated location of (485 m, 505 m) in case 1Ð3 is incorrect. In cases 1Ð1, 1Ð2, and 1Ð4, the estimated pumping rates are 200Ð99, 216Ð22 and 198Ð92 m3 day1, respectively, and the estimated pumping periods are 54Ð22, 43Ð98 and 59Ð57 min, respectively. The greatest relative errors are 8Ð11 % in the estimated pumping rate and 26Ð7 % in the estimated pumping period in case 1Ð2.
When the number of measured heads is increased to five, the source locations for cases 1Ð5 to 1Ð8 are all correctly predicted and the predicted pumping rates and pumping periods are all close to the correct solutions. The greatest relative errors are 1Ð89% in the estimated pumping rate in case 1Ð7 and 2Ð3% in the estimated pumping period in case 1Ð5. Obviously, one more mea-sured head is helpful in determining pumping source information. Thus, the number of measured heads used to solve the pumping source identification problem is suggested to be at least one more than the number of unknowns.
Scenario 2: Measurement error
Scenario 2 considers that the measured head has measurement error and the disturbed measured head is expressed as
h0i,obsDhi,obsð1 C Er ð RD3 7
where h0
i,obs is the disturbed measured head, Er is the
Table II. Required number of measured heads to determine pumping source information (pumping source is located at (495 m, 495 m) and the pumping rate and pumping period are 200 m3 d1and 60 min, respectively)
Case Observation wells Source location Pumping rate
(m3 day1) Pumping period (min) 1– 11 N-1, W-1, S-1, E-1 (495 m, 495 m) 200Ð99 54Ð22 1– 21 N-2, W-2, S-2, E-2 (495 m, 495 m) 216Ð22 43Ð98 1– 31 NW-1, SW-1, SE-1, NE-1 (485 m, 505 m) 226Ð90 72Ð42 1– 41 NW-2, SW-2, E-3, E-1 (495 m, 495 m) 198Ð92 59Ð57 1– 52 N-1, W-1, S-1, S-2, E-1 (495 m, 495 m) 200Ð21 58Ð62 1– 62 N-1, N-2, W-2, S-2, E-2 (495 m, 495 m) 198Ð73 60Ð02
1– 72 NW-1, SW-1, SE-1, E-2, NE-1 (495 m, 495 m) 196Ð22 61Ð26
1– 82 NW-2, SW-2, E-2, E-3, E-1 (495 m, 495 m) 198Ð99 60Ð38
1Number of measured heads D 4 2Numbers of measured heads D 5
normal deviate generated by the routine RNNOF of IMSL (IMSL, 2003). Three different values of Er, i.e. 1, 2 and 3 cm, are chosen in this scenario and ten cases for each level of measurement error are studied. In cases 2Ð1 to 2Ð30, the wells N-1, W-1, S-1, S-2, and E-1 depicted in Figure 2 are chosen as the observation wells. Table III shows that the source location is all correctly identified when Er D 1 and 2 cm in cases 2Ð1 to 2Ð20. The average of estimated pumping rate and pumping period are 199Ð42 m3 day1 and 60Ð41 min,
respectively. The greatest relative error is 1Ð24% in the estimated pumping rate and 3Ð72 % in the estimated pumping period in case 2Ð2 as Er D 1 cm. In addition, the average of the estimated pumping rate and pumping period are 199Ð24 m3 day1and 60Ð64 min, respectively. The greatest relative error is 2Ð84% in the estimated pumping rate and 11Ð55 % in the estimated pumping period in case 2Ð20 as Er D 2 cm. However, when the measurement error level is increased to 3 cm, only six out of ten cases obtain the correct source location as indicated in Table III. The average of the estimated pumping rate and pumping period in those six cases are 200Ð09 m3 day1 and 60Ð20 min, respectively. The greatest relative error is 3Ð65 % in the estimated pumping rate and 9Ð28% in estimated pumping period in case 2Ð30 as Er D 3 cm.
Scenario 3: Performance improvement
Based on the results in scenario 2, the estimated pump-ing source location may be incorrect and the pumppump-ing rate and pumping period may not be accurate if the level of measurement error is increased. Accordingly, cases 3Ð1 to 3Ð8 are designed to find an effective way to increase the accuracy of the identified result as the mea-surement error level is high. Cases 3Ð1 to 3Ð4, cases 3Ð5 and 3Ð6, and cases 3Ð7 and 3Ð8 use five, six and seven measured heads, respectively, to determine the pump-ing information. Note that all the measured heads are taken from different observation wells and are considered to have random measurement errors with Er D 3 cm. The estimated pumping source location, pumping rate and pumping period are correct in cases 3Ð1 and 3Ð2. However, the accuracy of the estimated results decreases
Table III. Results for measurement with error (pumping source is located at (495 m, 495 m) and pumping rate and period are
200 m3 day1 and 60 min, respectively)
Case Source location Pumping rate (m3day1 ) Pumping period (min) 2–1 (495 m, 495 m) 198Ð42 60Ð54 2–2 (495 m, 495 m) 197Ð53 62Ð23 2–3 (495 m, 495 m) 200Ð12 60Ð08 2–4 (495 m, 495 m) 199Ð43 60Ð12 2–5 (495 m, 495 m) 200Ð63 59Ð93 2–6 (495 m, 495 m) 198Ð12 61Ð12 2–7 (495 m, 495 m) 198Ð96 60Ð43 2–8 (495 m, 495 m) 200Ð72 59Ð21 2–9 (495 m, 495 m) 198Ð43 61Ð73 2–10 (495 m, 495 m) 201Ð87 58Ð75 2–11 (495 m, 495 m) 197Ð32 61Ð21 2–12 (495 m, 495 m) 200Ð88 59Ð83 2–13 (495 m, 495 m) 201Ð32 58Ð76 2–14 (495 m, 495 m) 199Ð08 60Ð19 2–15 (495 m, 495 m) 203Ð43 56Ð42 2–16 (495 m, 495 m) 196Ð43 62Ð48 2–17 (495 m, 495 m) 197Ð83 62Ð13 2–18 (495 m, 495 m) 203Ð53 56Ð83 2–19 (495 m, 495 m) 198Ð21 61Ð63 2–20 (495 m, 495 m) 194Ð32 66Ð93 2–21 (495 m, 495 m) 198Ð50 61Ð42 2–22 (505 m, 505 m) 205Ð51 54Ð32 2–23 (505 m, 505 m) 199Ð23 60Ð45 2–24 (495 m, 495 m) 202Ð67 58Ð48 2–25 (495 m, 505 m) 200Ð68 60Ð43 2–26 (495 m, 495 m) 204Ð13 55Ð46 2–27 (495 m, 495 m) 193Ð98 66Ð43 2–28 (495 m, 495 m) 194Ð15 64Ð98 2–29 (505 m, 495 m) 203Ð12 62Ð30 2–30 (495 m, 495 m) 207Ð13 54Ð43
Note: Maximum error for cases 2– 1 to 2– 10, 2– 11 to 2– 20 and 2– 21 to 2– 30 is 1, 2 and 3 cm, respectively.
with the increasing average distance from the observa-tion wells, as indicated in cases 3Ð3 and 3Ð4. However, if the number of the measured heads is increased to six or seven, the pumping source location, pumping rate and pumping period are all correctly estimated in cases 3Ð5 to 3Ð8 indicated in Table IV. Accordingly, it is found that the best way to obtain good estimated results is to employ more measured head data in the pumping source
Table IV. Required number of measured heads to determine the pumping source information when the maximum measurement error is 3 cm (pumping source is located at (495 m, 495 m) and pumping rate and period are 200 m3 day1and 60 min, respectively)
Case Observation wells Ave. dist.
(m)
Source location Pumping rate (m3 day1) Pumping period (min) 3–1 N-1, W-1, S-1, S-2, E-1 66Ð0 (495 m, 495 m) 198Ð50 61Ð42 3–2 N-1, N-2, W-2, S-2, E-2 86Ð0 (495 m, 495 m) 197Ð26 62Ð38 3–3 NW-1, SW-1, SW-2, SE-1, NE-1 169Ð7 (505 m, 485 m) 195Ð18 64Ð29
3–4 NW-2, SW-2, E-1, E-2, E-3 181Ð1 (455 m, 545 m) 208Ð31 50Ð41
3–5 NW-1, SW-1, S-2,SE-1, NE-1, E-2 120Ð5 (495 m, 495 m) 198Ð76 61Ð18
3–6 NW-2, SW-2, S-2 E-1, E-2, E-3, 150Ð9 (495 m, 495 m) 202Ð42 58Ð31
3–7 N-2, NW-1, SW-1, S-2,SE-1, NE-1, E-2 123Ð7 (495 m, 495 m) 199Ð15 60Ð73
3–8 N-2, NW-2, SW-2, S-2, E-1, E-2,E-3 157Ð9 (495 m, 495 m) 202Ð42 57Ð43
Table V. Number of observation wells required to determine pumping source information (pumping source is located at (495 m, 495 m) and pumping rate and period are 200 m3 day1and 60 min, respectively)
Case Observation wells Source location Pumping
rate (m3 day1) Pumping period (min) 4–11 N-1, W-1, S-1, E-1 (495 m, 495 m) 202Ð18 57Ð43 4–21 N-2, W-2, S-2, E-2 (495 m, 495 m) 198Ð66 61Ð29 4–31 NW-2, SW-2, E-3, E-1 (495 m, 495 m) 198Ð38 60Ð81 4–42 W-1, S-1, E-2 (495 m, 495 m) 198Ð64 61Ð32 4–52 N-1, W-1, S-1 (495 m, 495 m) 198Ð62 61Ð63 4–62 N-1, E-2, S-1 (495 m, 495 m) 198Ð58 61Ð58 4–72 N-1, W-1, E-2 (495 m, 495 m) 198Ð45 61Ð42 4–82 N-2, S-2, E-2 (495 m, 495 m) 197Ð88 60Ð92 4–92 N-2, W-1, S-2 (495 m, 495 m) 198Ð45 61Ð35 4–102 W-1, S-2, E-2 (495 m, 495 m) 198Ð48 61Ð72 4–112 N-2, W-1, E-2 (495 m, 495 m) 197Ð85 60Ð67 4–123 W-1, E-2 (565 m, 385 m) 216Ð79 50Ð48 4–133 S-1, E-2 (495 m, 495 m) 198Ð36 60Ð93 4–143 N-1, E-2 (495 m, 495 m) 197Ð17 62Ð04 4–153 W-1, S-1 (425 m, 415 m) 199Ð51 59Ð87 4–163 N-1, W-1 (495 m, 495 m) 198Ð12 61Ð42 4–173 N-1, S-1 (495 m, 495 m) 199Ð63 60Ð59 4–183 N-2, S-2 (495 m, 475 m) 203Ð43 58Ð46 4–193 S-2, E-2 (495 m, 495 m) 198Ð18 61Ð32 4–203 W-1, S-2 (565 m, 385 m) 216Ð67 50Ð14 4–213 N-2, E-2 (495 m, 495 m) 198Ð42 60Ð68 4–223 N-2, W-2 (495 m, 495 m) 198Ð32 61Ð32 4–233 W-1, E-2 (495 m, 495 m) 197Ð43 62Ð43
1Number of observation wells D 4 2Number of observation wells D 3 3Number of observation wells D 2
identification when the level of measurement error is high.
Scenario 4: Number of observation wells
In scenario 4, 23 cases are considered to examine the effect of using the observations measured at two different times on the determination of the pumping source information. The results are listed in Table V. In cases 4Ð1 to 4Ð3, eight observations measured at four wells are used to determine the pumping information. However, six observations measured at three wells in cases 4Ð4 to 4Ð11 and four observations measured at two wells in cases 4Ð12 to 4Ð23 are utilized. As the number of observation wells is four or three, the pumping source location, pumping rate and pumping period are
correctly determined. However, when the number of observation wells is two, eight out of twelve cases yield the correct results of pumping source information. Obviously, the number of observation wells should be at least three to determine the pumping source information.
HETEROGENEOUS AQUIFER
In this section, we estimate the pumping source informa-tion for problems with heterogeneous formainforma-tions repre-senting possible real-world aquifers. We consider three scenarios, i.e. scenarios 5 to 7, with different degrees
of variation in random hydraulic conductivity fields and conditioning conductivity data.
Site description
Field aquifer tests such as a slug test or pumping test are usually used to determine hydrogeologic parameters, i.e. hydraulic conductivity and storativity. These aquifer parameters obtained at specific locations can be used as the conditioning data. The heterogeneous aquifers may be characterized by three statistical parameters: the mean hydraulic conductivity K, the variance in hydraulic conductivity and the correlation length. The heterogeneous formations are assumed to have random conductivity fields which are spatially correlated and log-normally distributed. The correlation length is 100 m and the mean of the logarithm of hydraulic conductivity ln K is 12Ð206 m s1. The mean and standard deviation of ln K are denoted as y and y, respectively, where y D
ln K. Three different values of y, representing different
levels of aquifer heterogeneity, are considered: 0Ð5, 1Ð0 and 1Ð5 m s1. The dimensions of the aquifer are the same as those given previously, i.e. the length and width of the aquifer are both 1000 m and the thickness is 20 m. The locations of pumping source and observation wells are as depicted in Figure 2.
For scenarios 5 to 7, monitoring wells are N-1, W-1, S-1 and E-1. Table VI lists the hydraulic conductivities, which are assumed to be obtained from the aquifer test and adopted as the conditioning data for the heteroge-neous sites. All the random conductivity fields are gener-ated by the program SASIM of the geostatistical software, GSLIB (Deutsch and Journel, 1998). SASIM can pro-duce spatially correlated random conductivity fields by preserving the known mean and variance of ln K and known conductivities at specific locations. Therefore, the SASIM is chosen to generate random hydraulic conduc-tivity fields for heterogeneous aquifers with those three statistical parameters. Then MODFLOW-2000 is used with the generated conductivity field and under the same
Table VI. The hydraulic conductivity for the generated random conductivity field and the corresponding measured heads Scenario y (m s1) Observation well Hydraulic conductivity Measured head (m) (ð106 m s1) T D 30 min T D 60 min 5 0Ð5 N-1 5Ð40 22Ð847 21Ð506 W-1 2Ð41 25Ð148 24Ð499 S-1 6Ð79 24Ð700 24Ð256 E-1 7Ð07 22Ð751 21Ð483 6 1Ð0 N-1 15Ð70 23Ð589 22Ð579 W-1 1Ð39 25Ð330 24Ð808 S-1 4Ð64 24Ð689 24Ð247 E-1 6Ð18 21Ð978 20Ð602 7 1Ð5 N-1 10Ð30 23Ð790 23Ð017 W-1 0Ð59 25Ð526 25Ð236 S-1 5Ð58 24Ð790 24Ð483 E-1 18Ð30 22Ð228 20Ð996
hydrogeological conditions described previously to pro-duce the measured heads. The generated measured heads by MODFLOW-2000 with the pumping rate and period of 200Ð0 m3 day1and 60Ð0 min, respectively, are listed
in Table VI.
Monte Carlo simulation
Field conductivities, obtained from a slug test or pumping test, are generally available only at a few locations. The Monte Carlo simulation is employed to test the applicability of SA-MF. SASIM is used to produce 50 realizations of random conductivity fields with the same conditioning conductivity data and statistical parameters used in the process of generating the measured heads. The SA-MF is then employed to estimate the pumping source information for an aquifer with a realization, i.e. a conditioning random conductivity field. Note that in scenarios 5 to 7, the conductivity fields used to generate the simulated heads in identifying the pumping source information are different from the conductivity field used in generating the measured heads.
The estimated results show that the locations of pumping source are all correct in scenarios 5 to 7. In scenario 5, the average estimated pumping rate and pumping period are 197Ð11 m3 day1 and 61Ð11 min, respectively, when the value of y is 0Ð5 m s1.
The standard deviation of the estimated pumping rate and pumping period are 5Ð53 m3 day1 and 4Ð48 min,
respectively. The greatest relative error is 7Ð09% in the estimated pumping rate and 22Ð69 % in pumping period. In scenario 6, the average estimated pumping rate and pumping period are 193Ð74 m3 day1 and 63Ð37 min,
respectively, and the standard deviation of the estimated pumping rate and pumping period are 9Ð39 m3 day1
and 8Ð63 min, respectively, when y D1Ð0 m s1. The
greatest relative error is 11Ð06% in the estimated pumping rate and 38Ð01 % in the pumping period. If y is increased to 1Ð5 m s1, the average estimated
pumping rate and pumping period are 188Ð02 m3 day1
and 65Ð23 min, respectively. The standard deviation of the estimated pumping rate and pumping period are 10Ð59 m3 day1and 7Ð88 min, respectively. The greatest relative error is 16Ð07% in the estimated pumping rate and 38Ð52 % in the pumping period. The accuracy of the estimated pumping rate and pumping period obviously decreases with increasing y.
ACTUAL FIELD CASE
In the last scenario (scenario 8) the data obtained from a field pumping test is used to determine the pumping source information. The field pumping test was performed in 1987 at Yenliao, a site of the fourth nuclear power plant, located at the northeast corner of Taiwan. The dimension of the groundwater flow system is chosen as 200 m in both length and width. One pumping well, denoted as P, and four observation wells, denoted as A to D, were installed for the aquifer test, as shown in
Figure 3. Plan view of Yenliao site
Figure 3. The pumping well is located at (97Ð5 m, 97Ð5 m). Observation distances from the pumping well are A: 5 m, B: 10 m, C: 20 m and D: 20 m. Wells A to D are located at (92Ð5 m, 97Ð5 m), (87Ð5 m, 97Ð5 m), (77Ð5 m, 97Ð5 m), and (92Ð5 m, 110Ð0 m), respectively. According to the drill log obtained at P, the thickness of the aquifer is about 16 m and the formation composed of medium sand with a small amount of shale.
The surface topography is flat and the gradient is about 0Ð01. The water level is at a depth of 1 m under the ground surface. The groundwater flows toward the northeast, the upper and lower boundaries are assumed as no-flow and the left and right boundaries are considered as constant-head. The elevations of the hydraulic head at the left and right boundaries are 17 m and 15 m, respectively. The grid width and length are both 2Ð5 m; thus, the number of finite difference mesh is 80 ð 80. The aquifer is assumed as isotropic, the pumping well and observation wells are fully penetrated through the aquifer. The conductivity field is divided into two zones based on the results of pumping test data analysis. The estimated transmissivity and storage coefficient are 3Ð90 m2 day1 and 3Ð73 ð 104, respectively, in zone
1 and 20Ð51 m2 day1 and 3Ð83 ð 103, respectively,
in zone 2 as indicated in Figure 3. The pumping rate was 10 l min1 and the drawdown data were taken at the pumping times of 10 and 20 min. The drawdowns at wells A to D were 0Ð529, 0Ð130, 0Ð038 and 0Ð001 m, respectively, at pumping time of 10 min and 0Ð774, 0Ð250, 0Ð101 and 0Ð005 m, respectively, at pumping time of 20 min.
All the nodal points inside the problem domain are considered as the suspected pumping source locations except those nodes containing the observation wells. The upper bounds for the pumping rate and pumping period are set as 100Ð0 l min1and 60Ð0 min, respectively, and
Table VII. Identified results for the Yenliao site (pumping source is located at (97Ð5 m, 97Ð5 m) and the pumping rate and period
are 10 l min1 and 20 min, respectively)
Case Observation wells
Source location Pumping rate (l min1) Pumping period (min) 8–1 A, B, C (97Ð5 m, 97Ð5 m) 10Ð38 19Ð29 8–2 A, B, C, D (97Ð5 m, 97Ð5 m) 10Ð73 18Ð87 8–3 B, C, D (97Ð5 m, 97Ð5 m) 10Ð34 19Ð21 8–4 A, B, D (97Ð5 m, 97Ð5 m) 10Ð30 19Ð17
the lower bounds are both given as zero. The NS, NT, initial temperature and RTe are chosen as 20, 10, 5 and
0Ð8, respectively. All initial estimates are produced by a random number generator.
Four cases are used to test the performance of SA-MF in this scenario. Table VII lists the observation wells used in these four cases. The estimated source locations are all correct at (97Ð5 m, 97Ð5 m) in cases 8Ð1 to 8Ð4. The average estimated pumping rate and pumping period are 10Ð43 l min1 and 19Ð14 min, respectively. The greatest relative error is 7Ð3 % in the estimated pumping rate and 5Ð65% in the pumping period in case 8Ð2. Therefore, the proposed SA-MF has been shown to be able to estimate the pumping source information for a real field problem.
CONCLUSIONS
This study proposes a new approach, SA-MF, based on a 3D groundwater flow, MODFLOW-2000, and SA for solving the pumping source identification problem. The SA is first employed to generate the trial solu-tions randomly, i.e. pumping source location, pumping rate and pumping period and MODFLOW-2000 is then employed to simulate the hydraulic heads at the obser-vation wells. Eight scenarios were designed to test the accuracy and applicability of SA-MF in the determina-tion of pumping source informadetermina-tion on the homogeneous, heterogeneous aquifer systems. Four major conclusions were drawn.
First, the proposed approach can be used to deter-mine the pumping source information for problems in either synthetic heterogeneous aquifers or a real aquifer system and the estimations generally give good results. Second, SA-MF can provide good estimated results even if the initial guesses are generated by a random number generator. Third, at least five measured heads should be used to analyse the pumping source identification prob-lem and at least three observation wells are required for the effective determination of the pumping source loca-tion, the pumping rate and the pumping period. Fourth, if the level of measurement error is high, increasing the number of measured heads is helpful in determining the pumping source information correctly. In other words, the key for the successful determination of the pump-ing source information is to use more observed hydraulic heads.
ACKNOWLEDGEMENTS
This study was partly supported by the Taiwan National Science Council under the grant NSC94-2211-E-009-012. The authors thank two anonymous reviewers for their very detailed and insightful comments.
REFERENCES
Aral MM, Guan J. 1996. Genetic algorithms in search of groundwater pollution sources: advances in Groundwater Pollution Control and Remediation. NATO ASI Series 2 Environment . Kluwer Academic Publisher: Netherlands; 347– 369.
Aral MM, Guan J, Maslia ML. 2001. Identification of contaminant source location and modflow-gwt release history in aquifers. Journal of Hydrologic Engineering 6(3): 225– 234.
Bear J. 1979. Hydraulics of Groundwater . McGraw-Hill: New York. Cai X, McKinney DC, Lasdon LS. 2001. Solving nonlinear water
management models using a combined genetic algorithm and linear programming approach. Advances in Water Resources 24(6): 667– 676. Cooper VA, Nguyen VTV, Nicell JA. 1997. Evaluation of global optimization methods for conceptual rainfall-runoff model calibration. Water Science Technology 36(5): 53– 60.
Cunha MDC, Sousa J. 1999. Water distribution network design optimization: simulated annealing approach. Journal of Water Resources Planning and Management 125(4): 215– 221.
Deutsch CV, Journel AG. 1998. GSLIB: Geostatistical Software Library and User’s Guide. Second edition. Oxford University Press: Oxford. Doutherty DE, Marryott RA. 1991. Optimal groundwater management: I.
Simulated annealing. Water Resources Research 27(10): 2493– 2508. Guo JQ, Zheng L. 2005. A modified simulated annealing algorithm
for estimating solute transport parameters in streams from tracer experiment data. Environmental Modeling and Software 20(6):
811– 815.
Harbaugh AW, Banta ER, Hill MC, McDonald MG. 2000. MODFLOW-2000, The U.S. Geological survey modular ground-water model-user guide to modularization concepts and the ground-water process. Open-File Report 00– 92.
IMSL. 2003. Fortran Library User’s Guide Stat/Library. Volume 2 of 2. Visual Numerics Inc.: Houston.
Kuo CH, Michel AN, Gray WG. 1992. Design of optimal pump-and-treat strategies for contaminated groundwater remediation system using the simulated annealing algorithm. Advances in Water Resources 15(2): 95– 105.
Li L, Barry DA, Morris J, Stagnitti F. 1999. CXTANNEAL: an improved program for estimating solute transport parameters. Environmental Modelling and Software with Environment Data News 14(6): 607– 611. Lin YC, Yeh HD. 2005. THM species forecast using optimization method: genetic algorithm and simulated annealing. Journal of Computing in Civil Engineering 19(3): 248– 257.
Lo´aiciga HA. 2004. Analytic game-theoretic approach to ground-water extraction. Journal of Hydrology 297: 22– 33.
Mahinthakumar G, Sayeed M. 2005. Hybrid genetic algorithm-local search methods for solving groundwater source identification inverse problem. Journal of Water Resources Planning and Management
131(1): 45– 57.
Mahinthakumar G., Sayeed M. 2006. Reconstructing groundwater source release histories using hybrid optimization approaches. Environmental Forensics 7(1): 45– 54.
Marryott RA. 1996. Optimal ground-water remediation design using multiple control technologies. Ground Water 34(3): 425– 433. Marryott RA, Dougherty DE, Stollar RL. 1993. Optimal groundwater
management: 2. Application of simulated annealing to a field-scale contaminant site. Water Resources Research 29(4): 847– 860. Monem MJ, Namdarian R. 2005. Application of simulated annealing
(SA) techniques for optimal water distribution in irrigation canals. Irrigation and Drainage 54(4): 365– 373.
Pham DT, Karaboga D. 2000. Intelligent optimization techniques: genetic algorithms, tabu search, simulated annealing and neural networks. Springer-Verlag: New York.
Rayward-Smith VJ, Osman IH, Reeves CR, Smith GD. 1996. Modern heuristic search methods. John Wiley & Sons: New York.
Rizzo DM, Dougherty DE. 1996. Design optimization for multiple period groundwater remediation. Water Resources Research 32(8): 2549– 2561.
Romeo F, Sangiovanni-Vincentelli A. 1991. A theoretical framework for simulated annealing. Algorithmica 6(3): 302– 345.
Sayeed M, Mahinthakumar G. 2005. Efficient parallel implementation of hybrid optimization approaches for solving groundwater inverse problems. Journal of Computing in Civil Engineering 19(4): 329– 340. Shieh HJ, Peralta RC. 2005. Optimal in situ bioremediation design by hybrid genetic algorithm-simulated annealing. Journal of Water Resources Planning and Management 131(1): 67– 78.
Todd DK, Mays LW. 2005. Groundwater hydrology. Third edition. Wiley: New Jersey.
Tung CP, Tang CC, Lin YP. 2003. Improving groundwater flow modeling using optimal zoning methods. Environmental Geology 44: 627– 638.
Woodbury AD, Ulrych TJ. 1996. Minimum relative entropy inversion: theory and application to recovering the release history of a groundwater contaminant. Water Resources Research 32(9):
2671– 2681.
Yeh HD, Lin YC, Huang YC. 2007. Parameter identification for leaky aquifers using global optimization methods Hydrological Processes
21(7): 862– 872.
Zheng C, Wang P. 1996. Parameter structure identification using tabu search and simulated annealing. Advances in Water Resources 19(4): 215– 224.
