Optimal remediation with well locations and pumping rates selected as continuous decision variables

Journal of Hydrology 221 (1999) 20–42

Optimal remediation with well locations and pumping rates selected as continuous decision variables J. Guan, M.M. Aral* Multimedia Environmental Simulations Laboratory, School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0355, USA Received 6 May 1998; accepted 31 May 1999

Abstract The design of a pump-and-treat groundwater remediation system can be solved as an optimization problem. A common approach in this formulation is to minimize the total cost of the pump-and-treat system, while defining the locations and extraction or injection rates of the candidate pumping wells as continuous decision variables. With this choice, the degree of freedom added to the optimization problem yields significant improvements on the solution. In this approach coupled solution of groundwater simulation models and optimization algorithms are required. The repeated use of the groundwater simulation models throughout the optimization cycle tends to be numerically complex and computationally costly when the governing equations are nonlinear. To overcome this drawback, we propose a new computational procedure, identified as progressive genetic algorithm (PGA), to solve the optimal design problem. PGA is a subdomain method, which combines standard genetic algorithm with ground water simulation models in an iterative solution process and provides a powerful tool for the solution of highly nonlinear optimization problems. Numerical examples are included to demonstrate the feasibility and efficiency of the proposed algorithm. Applications indicate that the proposed approach provides a feasible alternative for the solution of nonlinear optimization problems in groundwater management. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Pump-and-treat; Remediation; Groundwater simulation; Optimization; Linearization; Progressive genetic algorithm; Decisionmaking

1. Introduction The remediation of contaminated groundwater is a complicated, time-consuming and expensive task. For the remediation of dissolved phase contaminant plumes, the pump-and-treat strategy is one of the most commonly used and successful remediation techniques. In the past decade, groundwater flow and contaminant transport modeling and optimization methods have been used in the design of these systems * Corresponding author. Fax: 1 1-404-8945111. E-mail address: [email protected] (M.M. Aral)

in an effort to produce a effective system at low cost and high performance. The pump-and-treat remediation systems involve installation and operation of a well network such that the contaminated groundwater is hydraulically contained and pumped out of subsequent treatment. The purpose of optimal design of such a system is to determine how many wells to install, where these wells should be located and what pumping rate is required from each well, while the objectives of the design are satisfied at minimum cost. This problem can be constructed as an optimization problem and solved by optimization methods embedded in groundwater simulation models.

0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(99)00079-7

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

An extensive review of typical models and algorithms used in the solution of this problem can be found in Gorelick (1983) and Aral et al. (1995). Some of the earlier applications in this field include studies by Gorelick et al. (1984), Atwood and Gorelick (1985), Ahlfeld et al. (1988a,b), Ahlfeld (1990), Ahlfeld and Sawyer (1990), Ratzlaff et al. (1992), and Aral et al. (1995), in which numerous progressive improvements were provided. The common deficiency in most of these earlier studies were either in the use of preselected locations for candidate pumping well sets, or in the use of linearized response matrix approach in the solution of the optimization problem, which restricted the optimal solutions significantly. If non-linear solutions were used, the solution processes were complex and required extensive computational time. The optimal design of a pump-and-treat remediation system tends to yields a highly nonlinear and nonconvex mathematical programming problem, no matter what objective function is selected. The linearization methods, which are based on linear response matrix approach, will result in significant computational error if they are used to solve the groundwater management problems in heterogeneous aquifers, while using complex objective functions. The nonlinear gradient-type algorithms require the derivatives of the objective function with respect to pumping rates, and may not guarantee global optimal solutions. They also require excessive computational effort. Therefore, in recent years, alternative algorithms abstracted from some physical processes, such as simulated annealing and genetic algorithms (GAs), were used in the solution of optimization problems in groundwater management. Dougherty and Marryott (1991) and Kuo et al. (1992) used simulated annealing algorithm to solve this problem. Since the simulated annealing method belongs to a class of probabilistic hill-climbing algorithms, it provided a means for escaping local minima during the solution of nonlinear optimization problem. This approach becomes competitive with gradient-type optimization methods. GAs are relatively new combinatorial search techniques based on the mechanics of natural selection and genetics (Holland, 1975; Goldberg, 1989; Davis, 1991). GAs have no special restriction for linearity, continuity and convexities of the problem studied and provide higher probability of obtaining the global optimal solution. Based on these advantages, GAs have

21

been widely used in the solution of optimization problems in numerous fields including groundwater management problems (Mckinney and Lin, 1994). In their research, authors indicated that the GA model provided solutions, which are as good as or better than, those obtained by linear and nonlinear programming approach. This conclusion was also demonstrated in a dynamic groundwater remediation management problem by Huang and Mayer (1997). In all of the studies reviewed above, the well locations were treated as implicit decision variables and determined by the pumping rates from a preselected candidate well location set, usually coincident with the nodes of the idealized domain. As is well known, different combinations of candidate pumping wells yield different optimal solutions for a given number of candidate pumping wells and the determination of optimal combination of candidate pumping wells is in itself a complicated problem. Approaching this problem, Wang and Ahlfeld (1994) proposed a formulation for the optimal design of aquifer remediation problem, where the well locations were determined explicitly by incorporating the spatial coordinates and extraction rates of wells into the objective function as decision variables. The management model combined the numerical simulation of groundwater flow and contaminant transport with nonlinear optimization model solver MINOS to select the optimal location and pumping rate by moving each well within the problem domain, unrestricted by nodal location. A hypothetical problem was used as an example to examine the numerical behavior of the formulation they proposed. Huang and Mayer (1997) also defined the pumping rates and well locations as decision variables to formulate the dynamic groundwater remediation management in homogeneous and heterogeneous aquifers. In their formulation, pumping rates were defined as continuous variables and well locations were limited within the discrete space containing the specified idealization nodes. GA was used to search for the optimal solution. Their numerical experiments demonstrated that dynamic pumping rates and well locations produce more cost-effective solutions relative to a static model and the fixed well model. As demonstrated in these two recent studies, the well locations play an important role in the optimal remedial design process. However, selecting the well location as an explicit

22

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

decision variable may result in other problems. In Wang’s approach, a barrier function represented by Hermite interpolation function was introduced in order to restrict the well locations within the problem domain. This has the disadvantage of increasing the chance of producing a local minima by the nonlinear optimization solver. In addition, since well locations may move anywhere in the domain, significant computational effort was required to solve this problem. In Huang and Mayer’s model, GA was directly used to search for the optimal solution of the remediation system design. Because groundwater simulation models were executed for each member in the population of each generation, the computational burden was excessive. On average, 250 h of CPU time on Sun Sparcstation 10 workstation was estimated to be required for the solutions of their problems (Huang and Mayer, 1997). In their study, well locations were selected from a discrete space of idealization nodes. This implies that the optimal well locations depend on the mesh size of the idealized domain. In order to obtain the best well locations, a small mesh size would be needed. Consequently, this would further aggravate the computational burden of the solution of this problem. In this study a general nonlinear optimization model is proposed for the design of a pump-andtreat groundwater remediation system, in which both the well locations and the pumping rates are defined as explicit continuous decision variables. In the formulation, three constraints were considered, which include concentration constraints for water quality requirement, velocity constraints for the advective plume containment, and balance constraint to equilibrate the extraction and injection rates. Considering that computational time is dominated by the execution of groundwater simulation models, a new combinatorial solution procedure, identified as progressive genetic algorithm (PGA), is introduced to solve the nonlinear optimization problem. This approach combines standard genetic algorithms with groundwater simulation models in a progressive iteration procedure within subdomains. Throughout the iteration cycle, the location, orientation and size of the subdomain changes continuously. Numerical examples in a randomly heterogenous aquifer were used to investigate the efficiency of the approach proposed and to illustrate the decision-making process of the remediation system design. Furthermore, the effect of well

locations, constraints, the number of the candidate pumping wells and different starting points on the solution were analyzed. 2. Optimization model The purpose of designing a pump-and-treat groundwater remediation system is to capture and remove the contaminant plume while containing it within a preselected capture zone during the remediation period. Mathematically, this problem can be formulated as an optimization problem. In this formulation, first a capture zone is defined around a contaminant plume observed at a site. Assuming that mw candidate wells will be installed in the aquifer, the solution domain is then divided into mw subdomains, which may be overlapping. Within each subdomain, a candidate well is located. Generally, the subdomains for the extraction wells are placed inside the capture zone and injection wells outside the capture zone. This choice facilitates the satisfaction of the velocity direction constraints, which are defined below. The objective of the optimization problem is to minimize the total cost of the pumpand-treat operation while satisfying the constraints. This problem can be mathematically stated as: minimize f …Q† ˆ

mw X

aj uQj u

…1†

jˆ1

such that Ci;T …Q; X; Y† # Cip ;

uLi # tan mw X

21

Vy i Vx i

i ˆ 1; 2; …; mc

…2†

! # uU i ;

i ˆ 1; 2; …; mv

Qj ˆ 0

…3†

…4†

jˆ1

Qj;min # Qj # Qj;max ;

j ˆ 1; 2; …; mw

…5†

xj;min # xj # xj;max ;

j ˆ 1; 2; …; mw

…6†

yj;min # yj # yj;max ;

j ˆ 1; 2; …; mw

…7†

where mw, mc and mv are the number of the candidate pumping wells, the number of concentration checkpoints and the number of velocity checkpoints

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

respectively; a j is the unit cost of injection and extraction rate in subdomain j; xj and yj are the x and y coordinates in subdomain j; Qj is the pumping rate at location …xj ; yj †, where Qj may be positive, negative or zero. A positive Qj indicates extraction rate and a negative Qj indicates injection rate. If Qj is equal to zero, then a pumping well is not necessary at location …xj ; yj †. If we define Q, X, Y as vectors with components Qj, xj and yj, respectively then Ci;T …Q; X; Y† is the contaminant concentration at concentration checkpoint i at the end of simulation period T, which is a function of Q, X, Y; Cip is the specified maximum concentration at concentration checkpoint i; Qj;min and Qj;max are the bounds of injection or extraction rates; …xj;min ; xj;max † and …yj;min ; yj;max † are the bounds of the subdomain for the well location in x and y directions, respectively. These bounds are defined within the solution domain and define the upper and lower bounds of a subdomain j, inside which the location of pumping well j may continuously vary. uLi and uUi are the angles between lower and upper bounds of the capture zone at the velocity checkpoint i and the x-axis (Ratzlaff et al., 1992; Aral et al., 1995). Vxi and Vyi are velocities in the x and y directions at velocity checkpoint i. In the optimization model described above, Eq. (2) reflects the water quality requirement, Eq. (3) is required for advective plume containment, and Eq. (4) is introduced to maintain an equilibrium between extraction and injection rates used in the groundwater pump-and-treat operation. The optimization model described above is a general model. The optimization model for fixed extraction and injection well locations can be considered to be a special case of this model. When the lower and upper bounds of x and y coordinates of a subdomain are given at the same coordinate location of the well, the above optimization model will reduce to a model in which the well locations are fixed. The optimization model described above can be further simplified, based on the requirements of the application. For example, if the equilibrium between total extraction and total injection rates are not necessary, the equality constraint given in Eq. (4) can be eliminated. Similarly, if one decides that water quality constraints are not necessary, then one may eliminate the constraint given in Eq. (2) and keep the velocity constraints given in Eq. (3).

23

3. Groundwater flow and contaminant transport simulation models In the optimization model given above, the contaminant concentrations in the aquifer, Ci;T …Q; X; Y† and velocities in the elements, Vxi and Vyi, are functions of the characteristics of the groundwater flow field and contaminant fate and transport processes in the aquifer. For steady state flow in a confined aquifer, the two-dimensional behavior of contaminants in the aquifer can be described by the following partial differential equations: 7·…BK·7h† 2

mw X

Qi d…xi ; yi † ˆ 0

…8†

iˆ1

with boundary conditions

a1 7hn ˆ b1 …h 2 hb † 1 g1 h ˆ z1 …x; y†

…9†

Given the pore Darcy velocity field V ˆ2

K 7h f

…10†

the governing equations for the contaminant transport problem can be given as 7·…D·7C† 2 V7·C 2

mw X Qi …C 2 Cp † 2C d…xi ; yi † ˆ R Bf 2t iˆ1

…11† The boundary conditions for this problem may be defined as:

a2 7Cn ˆ b2 …C 2 Cb † 1 g2 C ˆ g3

…12†

The initial condition can be defined as: C ˆ C0 …x; y†

…13†

where a 1, a 2, b 1, b 2, g 1, g 2, g 3 and z 1 are specified functions; h is the vertically averaged piezometric head in the aquifer; 7hn is the piezometric head gradient normal to the boundary; hb is the piezometric head value on the boundary; B is the saturated thickness of the aquifer; K is the aquifer conductivity tensor; mw is the number of pumping wells; Qi is the pumping rate for well i; d…xi ; yi † is the Dirac delta

24

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

function evaluated at …xi ; yi †; C is contaminant concentration; 7Cn is the gradient of the contaminant concentration normal to the boundary; Cb is the specified concentration value on the boundary; R is the retardation factor, defined as …R ˆ 1 1 rb Kd =f† in terms of a linear isotherm; r b is the bulk density of the solid media; Kd is the distribution coefficient; Cp is the concentration of the pumped water; f is the aquifer porosity; C0 is the initial distribution of contaminant concentration in the aquifer; D is the hydrodynamic dispersion tensor with components

Dxy ˆ

aL Vy2 1 aT Vx2 1 D0 ; uVu2

The contaminant concentrations at selected checkpoints are not given in terms of explicit expressions involving variables Q, X, and Y. These values can only be obtained from groundwater flow and contaminant transport simulation models as described above. In order to render Eq. (2) as an explicit function of variables Q, X and Y, a truncated Taylor series expansion will be used for the concentration variable in the neighborhood of …Q0 ; X 0 ; Y 0 † as described below: Ci;T …Q; X; Y† ˆ Ci;T …Q0 1 DQ; X 0 1 DX; Y 0 ; DY†

aL Vx2 1 aT Vy2 1 D0 ; Dxx ˆ uVu2 Dyy ˆ

4.1. Linearization of concentration constraints

< Ci;T …Q0 ; X 0 ; Y 0 † 1 …14†

…aL 2 aT †Vx Vy 1 D0 uVu2

V is the pore velocity vector; uVu2 is the Euclidian norm of the velocity vector; Vx and Vy are the components of velocity vector in x and y directions; and aL and aT are the longitudinal and transverse dispersivity coefficients, respectively; D0 is the coefficient of molecular diffusion. In this study, the governing equations described above are solved by standard finite element Galerkin procedures. Triangular elements were used in the idealization of the solution domain. Detailed description of this method can be found in a number of text books and papers, including Pinder and Gray (1977), Aral (1990a,b).

4. Progressive genetic algorithm The optimization model described above is highly nonlinear. Moreover, the concentration and velocity constraints of this optimization problem are not given in terms of explicit functions of the variables involved. Thus, the optimization model will be modified in an effort to develop a computationally efficient algorithm for the solution of this nonlinear problem. These modifications and the resulting optimization model are described below.

1

mw mw X X 2Ci;T 2Ci;T DQj 1 Dxj 2Qj 2xj jˆ1 jˆ1

mw X 2Ci;T Dyj 2y j jˆ1

i ˆ 1; 2; …; mc

…15†

Defining the vector CT …Q; X; Y† ˆ {C1;T ; C2;T ; …; Cmc;T }T will yield CT …Q; X; Y† ˆ CT …Q0 ; X 0 ; Y 0 † 1 ‰MQ ŠDQ 1 ‰MX ŠDX 1 ‰MY ŠDY

…16†

From Eq. (16), Eq. (2) can be simplified as: ‰MQ ŠDQ 1 ‰MX ŠDX 1 ‰MY ŠDY # bc

…17†

where [MQ], [MX] and [MY] are the …mc × mw † Jacobian matrices; bc is a vector with bc ˆ C p 2 C…Q0 ; X 0 ; Y 0 †. Elements of these matrices are evaluated by, mpi;j

Ci;T …p0 1 dpj † 2 Ci;T …p0 † 2Ci;T …Q; X; Y† ˆ < ; 2pj Dp

p ˆ Q; X; Y

…18†

where dpj ˆ …0; …; Dp ; 0; …; 0} are vectors with increments Dp evaluated at the jth location; DQ, DX, DY are the incremental variations for the vectors Q, X, Y, respectively.

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

and combining terms yields

4.2. Transformation of velocity constraints In order to handle the velocity direction constraints given in Eq. (3), coordinate transformations are used to transform the velocity direction constraints into the velocity magnitude constraints, at several preselected checkpoints on the capture zone boundary. The coordinate transformation used for this purpose were defined earlier and will not be repeated here (Ratzlaff et al., 1992; Aral et al., 1995). The resulting equations are: Vx sin uL 2 Vy cos uL # 0

…19†

2Vx sin uU 1 Vy cos uU # 0

…20†

From Eqs. (8) and (10), notice that Vx and Vy can be expressed as functions of Q; X; Y. Given an initial point {Q0 ; X 0 ; Y 0 }, Vx and Vy can be approximated by the truncated Taylor series expansion as

mw X jˆ1

1

! 2V x i 2V y i sin uLi 2 cos uLi DQj 2Qj 2Qj

mw X jˆ1

1

mw X jˆ1

mw X jˆ1

1

mw X

mw X

! 2V x i 2V y i 2 sin uUi 1 cos uUi Dxj 2x j 2x j ! 2V x i 2V y i 2 sin uUi 1 cos uUi Dyj 2y j 2y j

# Vxi …Q0 ; X 0 ; Y 0 †sin uUi 2 Vyi …Q0 ; X 0 ; Y 0 †cos uUi ; i ˆ 1; 2; …; mv

Vy …Q0 1 DQ; X 0 ; DX; Y 0 1 DY† < Vy …Q0 ; X 0 ; Y 0 † 1

mw X jˆ1

1

…26†

These equations can be expressed in matrix notation as

mw X

2V y 2V y DQj 1 Dx 2Qj 2x j j jˆ1

mw X 2V y Dy 2y j j jˆ1

…25†

! 2V x i 2V y i 2 sin uUi 1 cos uUi DQj 2Qj 2Qj

jˆ1

…21†

! 2V x i 2V y i sin uLi 2 cos uLi Dyj 2y j 2y j

i ˆ 1; 2; …; mv

jˆ1 mw mw X X 2V x 2V x DQj 1 Dx < Vx …Q0 ; X 0 ; Y 0 † 1 2 Q 2x j j j jˆ1 jˆ1

! 2V x i 2V y i sin uLi 2 cos uLi Dxj 2x j 2x j

# 2Vxi …Q0 ; X 0 ; Y 0 †sin uLi 1 Vyi …Q0 ; X 0 ; Y 0 †cos uLi ;

1

Vx …Q0 1 DQ; X 0 1 DX; Y 0 1 DY†

mw X 2V x 1 Dy 2y j j jˆ1

25

…22†

Similarly, partial derivatives of the velocity vectors are obtained from Vx …p0 1 dpj † 2 Vxi …p0 † 2V x i < i ; 2pj Dp

p ˆ Q; X; Y …23†

Vy …p0 1 dpj † 2 Vyi …p0 † 2V y i < i ; 2pj Dp

p ˆ Q; X; Y …24†

Substituting Eqs. (21) and (22) into Eqs. (19) and (20)

‰N1Q ŠDQ 1 ‰N1X ŠDX 1 ‰N1Y ŠDY # b1v

…27†

‰N2Q ŠDQ 1 ‰N2X ŠDX 1 ‰N2Y ŠDY # b2v

…28†

where ‰N1Q Š; ‰N1X Š; ‰N1Y Š; ‰N2Q Š; ‰N2X Š; ‰N2Y Š are …mv × mw † matrices corresponding to DQ, DX and DY and the elements of these matrices can be calculated as follows: n1pi;j ˆ

2V x i 2V y i sin uLi 2 cos uLi ; 2pj 2pj

n2pi;j ˆ 2

2V x i 2V y i sin uUi 1 cos uUi ; 2pj 2pj

p ˆ Q; X; Y …29†

p ˆ Q; X; Y …30†

26

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

b1vi ˆ 2Vxi …Q0 ; X 0 ; Y 0 †sin uL 1 Vyi …Q0 ; X 0 ; Y 0 †cos uL …31† b2vi ˆ Vxi …Q0 ; X 0 ; Y 0 †sin uU 2 Vyi …Q0 ; X 0 ; Y 0 †cos uU …32† 4.3. The modified optimization model On the basis of these transformations, the modified optimization model can be restated as follows. Let 2 3 2 3 MQ MX 6 7 6 17 6 7 7 ‰A2 Š ˆ 6 ‰A1 Š ˆ 6 NQ1 7; 4 NX 5; 4 5 NX2 8 9 bc > > > < > = {b} ˆ b1v > > > : 2> ; bv

NQ2 2

MY

3

6 17 7 ‰A3 Š ˆ 6 4 NY 5; NY2

…33†

where [A1], [A2], [A3] are ‰…mc 1 2mv † × mw Š matrices and {b} is a …mc 1 2mv † vector. Then the optimization model given in Eqs. (1)–(7) in terms of variables Q; X; Y, can now be given in the modified model in terms of the variables {DQ, DX, DY} as minimize f …Q† ˆ

mw X

aj uQ0j 1 DQj u

…34†

jˆ1

DXmax ˆ Xmax 2 X 0

DYmin ˆ Ymin 2 Y 0 ;

DYmax ˆ Ymax 2 Y 0

‰A1 ŠDQ 1 ‰A2 ŠDX 1 ‰A3 ŠDY # {b}

…35†

…40†

4.4. The computational model In order to improve the computational process, the modified optimization model given above will be further modified as follows: (i) The slack variables, si $ 0 for i ˆ 1; 2; …; …mc 1 2mv †areintroducedtoEq.(35)sothat it is expressed as an equality constraint; (ii) since the extraction and injection rates of the pump-and-treat system have only …mw 2 1† degree of freedom, the extraction rate Qk can be selected as the dependent variable whichiscomputedbyEq.(36)and;(iii)the objective function is modified by introducing the penalty terms for the slack variables, si and Qk, and alsothe minimizationof objective function is changed to maximization of the negative objective function. These modifications result in the following computational model: maximize 2 f …Q† ˆ 2

mw X

aj uQ0j 1 DQj u

jˆ1

2bk ‰min…0; Qk 2 Qk;min ; Qk;max 2 Qk †Š2 X 2 bi ‰min…0; si †Š2

…41†

;i

such that ‰A1 ŠDQ 1 ‰A2 ŠDX 1 ‰A3 ŠDY 1 S ˆ {b}

such that,

mw X

DXmin ˆ Xmin 2 X 0 ;

Qk ˆ 2

mw X

…Q0j 1 DQj †

…42† …43†

jˆ1;j±k

…Q0j 1 DQj † ˆ 0

…36†

DQmin # DQ # DQmax

DQmin # DQ # DQmax

…44†

…37†

DXmin # DX # DXmax

…45†

DXmin # DX # DXmax

…38†

DYmin # DY # DYmax

…46†

DYmin # DY # DYmax

…39†

S$0

…47†

jˆ1

where DQmin, DQmax, DXmin, DXmax, DYmin and DYmax are the lower and upper bounds for DQ, DX and DY, which are computed by the following equations: DQmin ˆ Qmin 2 Q0 ;

DQmax ˆ Qmax 2 Q0

where S ˆ {s1 ; s2 ; s3 ; …; s…mc 12mv † }T ; b k and b i are penalty coefficients with large positive constants compared to the value of the objective function. It should be pointed out, that the penalty term introduced above is different then the penalty function

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

methods used in forced constraint optimization. In our application, penalty coefficients are fix constants and do not need to change during iterative solution. 4.5. Computational procedures The optimization model described above looks like a linear constrained programming problem which has been used in earlier applications. However, since it is based on Taylor series expansion, it has some differences with respect to classical linear constrained optimization models. First, all coefficient matrices and vectors in the constraints are a function of the solution step {Q0 ; X 0 ; Y 0 } and will change as {Q0 ; X 0 ; Y 0 } changes; second, the proposed formulation is valid within the neighborhood of the solution {Q0 ; X 0 ; Y 0 }. This implies that if the subdomain size is large, the accuracy of the solution will decrease. Therefore, a new approach, identified as progressive genetic algorithm (PGA), is introduced to solve this model. The fundamental idea in PGA is that, at a given solution step, starting from an initial condition, one builds up the approximate linear constrained model as described above. Then, standard GA is used to search the optimal solution in the neighborhood of {Q0 ; X 0 ; Y 0 } within the subdomain. An iteration algorithm is then applied to obtain the final solution, utilizing the local optimal solution as the new initial solution around which a new subdomain is developed. Based on this idea, the computational steps of this process can be given as follows. Step 1. Definition of starting point {Q0 ; X 0 ; Y 0 } and solution domain regularization: According to the requirements of the design of the remediation system, the solution domain is divided into several subdomains. Within each subdomain a potential well is placed. The initial location and the initial injection and extraction rates of pumping wells can be determined by designer’s experience and knowledge of the hydrogeology of the site under study. As expected, a good initial starting point will reduce the computation time to converge to an optimal solution. The iteration index is defined as k ˆ 0, and objective function value, f0, is selected as a large negative constant. Step 2. Development of a modified optimization model: Groundwater flow and contaminant transport simulation models are run, to generate Jacobian matrices and vectors defined in the equations given

27

above. In general, when the number of candidate pumping wells is mw, groundwater simulation models need to be run for …3 × mw 1 1† times in order to build the modified optimization model using Eqs. (34)–(39). Step 3. Development of the computational model: In the computational model given by Eqs. (41)–(47), slack variables si ;i and Qk are chosen as dependent variables, which are computed from Eqs. (42) and (43). The variables DQ, DX and DY are treated as independent variables which are generated by GA. Step 4. Determine the search subintervals of independent variables: In order for the modified model to be valid, the GA search process is limited to the neighborhood of {Q0 ; X 0 ; Y 0 }. This neighborhood may be determined by the following simple method. First, choose a fixed interval for each independent variable pj;max 2 pj;min dpj ˆ ; p ˆ Q; X; Y …48† k1 The intervals ‰DQlower ; DQupper Š, ‰DXlower ; DXupper Š and ‰DYlower ; DYupper Š are then defined by: Dpj;lower ˆ max{Dpj;min ; 20:5e2gk dpj }; p ˆ Q; X; Y Dpj;upper ˆ max{Dpj;max ; 0:5e2gk dpj }; p ˆ Q; X; Y …49† The intervals are approximately taken as the subdomain of {Q0 ; X 0 ; Y 0 }; where k1 is some integer (e.g. k1 ˆ 5) k is the iteration index; g is a positive coefficient (e.g. g ˆ 0:001) identified as interval reduction coefficient. The subdomain defined is a regular polyhedron with the center located at {Q0 ; X 0 ; Y 0 } and the volume of the polyhedron decreases as the number of iterations increases. When k is zero, the subdomains are bounded by intervals d Qj, d xj and d yj and as the iteration process Approaches the optimal solution, the subdomains will gradually reduce as a factor of e2gk . Based on computational experiments we have conducted, our experience indicates that this approach significantly improves the computational efficiency and convergence rate of the algorithm proposed in this study. Obviously, the subdomains defined by Eq. (49) should satisfy ‰DQlower ; DQupper Š # ‰DQmin ; DQmax Š ‰DXlower ; DXupper Š # ‰DXmin ; DXmax Š ‰DYlower ; DYupper Š # ‰DYmin ; DYmax Š

…50†

28

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

Step 5. GA formulation: In GA, each independent variable is encoded as a k-bit binary, and all independent variables are connected as a binary string to simulate biological chromosome. These independent variables will be searched by GA in the subdomains ‰DQlower ; DQupper Š, ‰DXlower ; DXupper Š and ‰DYlower ; DYupper Š defined by Eq. (49). The binary variables can be mapped into decimal ones by the relationships: Dpi ˆ Dplower;i 1 n

Dpupper;i 2 Dplower;i 2` 2 1

;

…51†

operator randomly changes bits within newly created strings, based on a mutation probability value. Mutation maintains variability in the population, and reduces the chance that the population will prematurely converge on one possible suboptimal solution. The GA computation process continues until convergence is achieved. The resulting local optimal solution is identified as DQ p, DX p and DY p and the new local optimal solution for the original variables is determined as, Q1 ˆ Q0 1 DQp ;

X 1 ˆ X 0 1 DX p ;

…52†

p ˆ Q; X; Y

Y 1 ˆ Y 0 1 DY p

where n is the integer with respect to the binary; ` is the length of the binary. Step 6. Search for local optimal solution: Compute the intermediate local optimal solution Qp ; X p ; Y p by GA in the subdomains defined above (Goldberg, 1989; Holland, 1975). In this process, first the initial population for independent variables, which follows a uniform distribution and contains the initial solution, is generated. Second, the dependent variables slack variables and Qk are computed from Eqs. (42) and (43). The fitness of each member in the current population is estimated by the objective function (41). Third, choosing parents from current population forms a mating pool based on member’s fitness. The larger the fitness is, the larger the probability that the corresponding member is selected as parents. Thus, members with higher fitness values are more likely to be selected for mating. Finally, the new population, whose members have higher quality is reproduced by using the standard GA operators such as direct selection, crossover and mutation. In direct selection operation, the members with highest fitness directly enter the next population at some proportion. These members may also be selected as parents to the mating pool. Thus, using this operator, the best members are protested, and at the same time, this selection assures the monotone increase of the objective function values of the best member during each generation. Crossover is the principal genetic operator. It is the splicing of two parent members from mating pool, at a randomly chosen point. Then keeping the head substrings the same and exchanging the tail substrings creates two children strings, each of whom is made up of “genetic material” from both parents. Mutation

and the corresponding objective function value is f1. Step 7. Evaluation of the convergence of iterations: As is the case for most iterative algorithms, various criteria may be used to test the convergence of the iterative algorithm. For the PGA presented in this study, we chose the magnitude of the relative error of the objective function values between two consecutive iterations as the convergence criteria. If the two consecutive objective function values satisfy the criteria given in Eq. (53), uf1 2 f0 u #1 f0

…53†

then Q1 ; X 1 ; Y 1 is taken as the final solution, otherwise the solution sequence defined above continues with the new starting point defined as: Q 0 ˆ Q1 ; f0 ˆ f1 ;

X0 ˆ X1; k ˆk11

Y 0 ˆ Y 1;

…54†

After this point the computation sequence will start again from step (2) and repeat the computation steps through (7) until a satisfactory convergence level is achieved. In Eq. (53), 1 is a preselected relative allowable error. From the procedure described above, it can be seen that PGA is an iterative algorithm, which uses a series of linear approximations to solve the original nonlinear problem. Matrices [A1], [A2], [A3] as well as {b} in the modified model are actually nonlinear response matrices. Obviously, this approximation is better than solving the problem with linear response matrices as was done in other studies. In solving the

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

29

Fig. 1. Hydraulic conductivity distribution in a randomly heterogeneous aquifer (ft/h).

local optimal solution, infeasible solutions will be eliminated or improved upon using the standard GA operators as a matter of course. Similar to other iterative methods, if the initial subdomain is small, the final solution obtained by PGA may not be a global optimal solution. However, this problem can be avoided easily by selecting a large initial subdomain. PGA has a higher probability in obtaining global optimal solution when compared to other gradient-type nonlinear programming methods. This is because the PGA is an iterative global search algorithm, which combines global/sequential random search approach with standard global evolution approach of GA.

5. Applications A randomly heterogeneous isotropic confined aquifer is selected to demonstrate the performance of the optimization model developed in this study. Triangular finite element discretization and finite element Galerkin method were used to solve the groundwater flow and contaminant transport simulation models.

The confined aquifer with dimensions of 3900 ft by 3300 ft was idealized by 1360 nodes and 2574 elements yielding a mesh size of 100 ft. Impervious boundaries are imposed on the north and south of the aquifer. Groundwater flows from west to east due to the constant piezometric head boundary conditions along the west and east of the aquifer. A spatially correlated random field for hydraulic conductivity was generated, using random simulation approach. The resulting hydraulic conductivity distribution is shown in Fig. 1. Other relevant hydraulic parameters used in the groundwater flow and the contaminant transport simulation models are given in Table 1. Piezometric head distribution and undisturbed groundwater flow field obtained by the finite element solution of the governing equations are shown in Fig. 2. It is assumed that groundwater in the aquifer was contaminated by a constant contaminant concentration of 1000 units at …x ˆ 1000 ft; y ˆ 1600 ft†. The three-year extent of the contaminant plume, obtained from simulation models, is shown in Fig. 3. This plume was selected as the initial plume extent in all applications discussed in this study. In order to

30

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

Table 1 Parameters used in groundwater simulation models Parameter

Value

Time integration step Longitudinal dispersivity Transverse dispersivity Diffusion coefficient Angle of anisotropy Retardation coefficient Effective porosity Aquifer thickness Piezometric head boundary condition at x ˆ 0:0 ft; 0:0 # y # 3300 ft Piezometric head boundary condition at x ˆ 3900:0 ft; 0:0 # y # 3300 ft

2190 h 50.0 ft 30.0 ft 0.0 ft 2/day 08 1.0 0.25 30.0 ft 230.0 ft 130.0 ft

contain the contaminant plume, a capture zone was designed as shown in Fig. 3. Fourteen concentration checkpoints, indicated by solid diamonds, were also placed around this plume for the water quality check points. The maximum allowable concentration at the concentration checkpoints was defined to be two percent of the maximum source concentration, i.e. 20 units of concentration. The contaminant source

was assumed to be removed at the beginning of remedial action. The objective of remediation is to contain and remove most of the contaminant plume by using a pump-and-treat system within a remediation period of three years, throughout which selected constraints are satisfied. In the applications below, the unit cost of a pumping well operation was selected as one and installation costs were not considered. Thus, in this case, the objective function minimizes the total discharge. PGA was applied to search for the optimal pumping rates and locations of the candidate pumping wells. Parameters associated with PGA are given in Table 2. The optimal design of the remediation system is a decision-making process. Sometimes, the solution obtained from a method may be optimal in a mathematical sense, however this may not be a practical solution in application. For example, in some application discussed below, pumping rates for some pumping wells obtained from the optimization model are very small and may be neglected. But the solution, which excludes these small pumping rates may not satisfy the constraints if the problem is not resolved. Therefore, the redesign of the remediation system is

Fig. 2. Piezometric head distribution and the undisturbed groundwater flow field in a randomly heterogeneous aquifer.

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

31

Fig. 3. Initial contaminant plume, capture zone, and concentration checkpoints in a randomly heterogeneous aquifer.

necessary, based on the previous computational results with a revised number of wells in the solution. In general, the final decision is made from this design–revise–redesign process. For example, the pump-and-treat remediation system for the above problem was initially designed with five candidate pumping wells. Then, based on this solution, two candidate pumping wells were used to remove the contaminant plume. Finally, one arrives at the conclusion that only one pumping well may be sufficient to solve this problem. In this process, the effects of well locations, constraints, number of candidate pumping wells and different starting points on the total

discharge were analyzed. In the following figures, solid triangles indicate the injection well locations, solid circles indicate the extraction well locations and solid diamonds indicate the concentration checkpoint locations. In Tables 3 and 4, where the numerical results are summarized, the three number data vector in the “Control Constraints and Initial Conditions” column indicate the values of the initial discharge and the initial x and y coordinates of each well, respectively. The negative discharge in this vector indicates an injection rate and positive discharge indicates an extraction rate.

Table 2 Parameters used in progressive genetic algorithm

5.1. Preliminary design of remediation system

Parameter

Value

Population size String length Probability of crossover Probability of adaptive mutation Probability of direct selection Interval reduction coefficient

30 10 bits 0.95 0.05 0.05 0.01

The optimal design process starts with a five candidate pumping well system, which includes three candidate injection wells and two candidate extraction wells. The overlapping subdomain boundaries of the candidate wells are shown in Fig. 4. Each candidate pumping well can be moved to any location within its subdomain during the iterative process. This problem was investigated in two cases. First, a fixed well

32

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

Table 3 Computational results for five candidate pumping wells Case

(a)

(b)

(c)

Control constraints and initial conditions

Well location

Injection rate (ft 3/h)

X (ft)

Y (ft)

900.00 2200.00 2200.00 1200.00 1900.00

1200.00 1200.00 2000.00 1700.00 1600.00

88.33 6.67 324.17

1. Well locations are decision variables; 2. Concentration, velocity and balance constraints; 3. Starting points: (2299.5, 900.0, 1200.0) (2167.2, 2200.0, 1200.0) (277.6, 2200.0, 2000.0) (201.5, 1200.0, 1700.0) (342.8, 1900.0, 1600.0)

991.10 2298.44 2338.22 1361.32 1964.27

1241.52 1294.16 1820.54 1955.83 1612.81

2.14 32.49 270.35

1. Well locations are decision variables; 2. Concentration, velocity constraints; 3. Starting points: (2299.5, 900.0, 1200.0) (2167.2, 2200.0, 1200.0) (277.6, 2200.0, 2000.0) (201.5, 1200.0, 1700.0) (342.8, 1900.0, 1600.0)

978.00 2255.87 2258.78 1388.36 1991.46

1291.89 1229.97 1814.37 1772.59 1620.54

10.08 0.40 86.03

1. Well locations are fixed; 2. Concentration, velocity and balance constraints.

location and then a moving well location case is solved for comparison purposes. Table 3 lists the computational conditions under which optimization computation was implemented and the numerical results obtained. 5.1.1. Effect of well locations In Case (a) of Table 3, the well locations were fixed, optimization was performed only for pumping rates. The resulting total discharge is equal to 838.45 ft 3/h. In Case (b) of Table 3, the well locations were taken as continuous decision variables along with pumping rates and the optimization computation was started from the fixed well locations. In this case, the total discharge decreases to 609.94 ft 3/h. This is 27.3% less than the case with fixed well locations. In the fixed well location problem, the dimensionality of the decision space is mw, if the candidate pumping well

Extraction rate (ft 3/h)

Total discharge (ft 3/h)

838.45

0.10 419.07 609.94

3.54 301.42

573.60

113.56 363.39

locations are taken as decision variables, the decision space extends to 3 × mw dimensions and a better solution is naturally searched in this space. Numerical results obtained for Cases (a) and (b) are shown in Figs. 5 and 6. From these figures, it can be noted that all the velocity vectors around the boundary of the capture zone point inward towards the capture zone. In the fixed well location case, the only way to control the directions of the velocity vectors is by increasing the discharge. However, in the moving well location case, an additional degree of freedom is available to satisfy this constraint, which yields a better solution. However, since the total discharge in the fixed well solution is larger than in the moving wells case, the former removes more contaminant volume than the latter, and both satisfy all constraints of the problem. Based on the solutions of this example problem, it can be concluded that the candidate well

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

33

Table 4 Computational results for one and two candidate pumping wells Case

(a)

(b)

(c)

Control constraints and initial conditions

Well location

Injection rate (ft 3/h)

X (ft)

Y (ft)

1. Well locations are decision variables; 2. Concentration, velocity constraints; 3. Starting points: (21500.0, 2200.0, 2000.0) (21500.0, 1900.0, 1600.0)

2229.43 1669.45

1486.12 1814.00

28.41

1. Well locations are decision variables; 2. Concentration, velocity constraints; 3. Starting points: (21500.0, 2200.0, 1300.0) (1500.0, 2000.0, 1600.0)

2220.45 1659.97

1471.98 1817.45

19.98

1. Well locations is decision variables; 2. Concentration, velocity constraints; 3. Starting points: (21500.0, 2000.0, 1600.0)

1612.95

1658.74

location is very important in the design of the remediation system and optimal arrangement of the candidate pumping well locations will significantly decrease the total discharge used in pump-and-treat operation. 5.1.2. Effect of constraints The number and type of constraints have a significant effect on optimal solutions. In general, the larger the number of constraints, the smaller the decision space becomes. In the example discussed above three types of constraints were considered, which include concentration constraints, velocity constraints and the injection and extraction balance constraint. In some cases the injection wells may not be necessary. In the application given in Case (c) of Table 3, the problem discussed earlier was solved without the injection and extraction balance constraint. In this case, the resultant total discharge decreases to 573.60 ft 3/h. For this case, although the total discharge was reduced by only 36.34 ft 3/h compared to the results of Case (b), the total injection rate was reduced by 96.51 ft 3/h from 304.98 ft 3/h and the total extraction rate was increased to 476.96 ft 3/h from

Extraction rate (ft 3/h)

Total discharge (ft 3/h)

639.54 611.13

636.57 616.59

–

711.83

711.83

304.98 ft 3/h. This indicates that extraction wells play a more important role in the removal of the contaminant and injection wells may not be necessary. In the design of pump-and-treat system, one should determine whether the injection and extraction rate balance constraint is needed as a design requirement. If such a requirement is not necessary, then this constraint may be avoided to yield better optimal solutions. The importance of the concentration and velocity constraints on the capture and removal of contaminated groundwater has been investigated earlier by Aral and Guan (1997). 5.2. Improvement on remediation system design From the results of five candidate pumping wells, it can be seen that some pumping rates are very small and remedial action was mainly undergone by two pumping wells. Obviously, such a solution is not reasonable in practical applications, although it is optimal mathematically. Thus, two candidate pumping wells may be enough to complete the remedial action for this example. Therefore, the pump-andtreat remediation system needs to be further improved. Based on the previous results, next, we

34

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

Fig. 4. Boundaries and initial locations of five candidate pumping wells.

will select two candidate pumping wells for the remediation system. The boundaries of the candidate wells were defined as …1400; 2100† × …1100; 2100† for the extraction well, and …2000; 2400† × …1000; 2200† for the injection well. In this application, the concentration checkpoints and the capture zone are the same as the five candidate pumping wells, and only the concentration and velocity constraints were considered. Two different starting point sets for moving well locations, given in Cases (a) and (b) of Table 4, were addressed. Table 4 lists the computational conditions and the numerical results. 5.2.1. Effect of number of candidate pumping wells The number of candidate pumping wells has some effect on total discharge. In general, decreasing the number of pumping wells will increase the total discharge, and decrease the costs of well installation and management at the same time. Although in these applications, we did not consider these costs, we can still make decision by the comparison of computational results. Based on the results of the two candidate pumping wells example, the total discharge in

Case (a) of Table 4 has increased about 65.94 ft 3/h when compared to the results of pattern in Case (c) of Table 3. This increase is only 11% and more contaminant will be extracted for this case. Even before considering the operational costs of five wells versus two wells, one may conclude that the solution with two pumping wells would be more reasonable than the one with five pumping wells. In Fig. 7 the piezometric head distribution, groundwater flow field and the three year plume extent is shown for this pattern. Analyzing the pumping rates in Case (a) of Table 4, the injection rate is still very small. This implies that the remediation task may be implemented by one pumping well. Next we will investigate this problem. For this example, the location of the pumping well is allowed to move within the boundary …1200; 2200† × …1300; 2100†; and the well was initially placed at point …2000; 1600† ft. The numerical solution for this example is given in Case (c) of Table 4. In this case, the total discharge is 711.8 ft 3/ h and 72.26 ft 3/h larger than in Case (a) of Table 4, which uses two wells. Using this strategy, one can remove almost the total contaminant plume during a

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

35

Fig. 5. (a) Three-year plume extent; (b) piezometric head distribution and groundwater flow field. Computational conditions: five candidate pumping wells, fixed well locations, velocity, concentration and balance constraints.

36

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

Fig. 6. (a) Three-year plume extent; (b) piezometric head distribution and groundwater flow field. Computational conditions: five candidate pumping wells, well locations as decision variables, velocity, concentration and balance constraints.

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

37

Fig. 7. (a) Three-year plume extent; (b) piezometric head distribution and groundwater flow field. Computational conditions: two candidate pumping wells, well locations as decision variables, velocity and concentration constraints.

38

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

Fig. 8. (a) Three-year plume extent; (b) piezometric head distribution and groundwater flow field. Computational conditions: one candidate pumping well, well location as decision variables, velocity and concentration constraints.

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

39

Fig. 9. Iterative solution routes of well locations from two sets of starting points for two candidate pumping wells.

three-year period, as shown in Fig. 8(a). The flow field for this case is shown in Fig. 8(b). 5.2.2. Effect of starting points The starting points or initial conditions of iterative algorithms may affect the convergence of the algorithm. Generally, for a robust algorithm, the convergence should be independent of the selection of starting points. In order to investigate the sensitivity of the proposed solution algorithm to starting points, two sets of starting points for the extraction and injection wells in the two candidate pumping wells were selected, to design the pump-and-treat system. The computational results obtained for these case are given in Cases (a) and (b) of Table 4. The total discharges obtained for these cases are 639.54 and 636.57 ft 3/h and are not significantly different. The movement route of well

locations during the iteration is shown in Fig. 9. The locations of injection and extraction wells from different starting points converged almost to the same locations. The maximum differences of the locations were 9.48 and 14.14 ft for the extraction and injection well, respectively. From Fig. 9, it can be seen that the extraction well location converges better than the injection well location, and in the injection well location, some oscillations exist as the injection rate decreases during the iterations. The main reason that led to this phenomenon, is the small injection rate, which does not have significant effect on the iterative solution objective function value based on the variation of the injection well location. In Fig. 10 the relationships between the pumping rates and iterations for these two sets of starting points are shown. From these results, it is clearly seen that

40

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

Fig. 10. Relationships between pumping rates and iterations from two sets of starting points for two candidate pumping wells.

the iterative solution converges very fast initially, and as the solution approaches the optimal solution the convergence rate decreases. 5.3. Computational cost For the optimal design of pump-and-treat remediation systems, the computational cost measured by computational time is a very important issue. As is well known, the solution of the groundwater simulation models dominates the computational time, estimated to be about 90% of CPU time on average (Huang and Mayer, 1997). As previously stated, the groundwater simulation models need to be run for …3 × mw 1 1† times during each iteration of the optimization model when PGA is used, where the well locations are defined as decision variables. The examples discussed above were solved on a Silicon Graphics Unix platform. In these applications, each

iteration requires approximately 5 min for the two candidate pumping well case and about 15 min for the five candidate pumping well case. The number of iterations depends on the preselected relative allowable error and the starting point of the algorithm. If the allowable relative error of objective function values between two iterations is given as one percent, the average number of iterations required was about 10 as shown in Fig. 10. This implies that about 50 min was required to solve the two candidate pumping well problem and about 150 min was required to solve the five-candidate well problem. Dougherty and Marryott (1991) compared the computational performance of MINOS, which was reported by Ahlfeld (1990) and simulated annealing algorithms using relationship between relative CPU times and the log peak reduction fraction and found that the relative CPU requirements for simulated annealing increase at a slower rate than the gradient-based method as the level of

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

contaminant removal is increased. In another paper, Kuo et al. (1992) used simulated annealing algorithm to design optimal pump-and-treat strategies for contaminated groundwater remediation for fixed well locations with a solution domain dimensions of 1100 ft × 800 ft. In their algorithm, the simulation models were run for more than 4000 times in order to obtain the optimal solution. Based on the simulation models used in our application, each simulation run requires about 5/7 min. If the simulated annealing algorithm was used by our application, a total of 4000 × 5=7 min (more than 47 h) would be required on a Unix based platform. Huang and Mayer (1997) directly employed GA to solve the remediation system design problem. In their study, the population size was selected as 30, and about 20 generations were needed for the optimal solution. This implies that about 600 runs of the groundwater simulation models were required. Comparing these results, the advantage of PGA is apparent for the solution of groundwater management problems. The inefficiency of these algorithms is caused mainly by the solution time requirement of the repeated execution of the groundwater simulation algorithms. Therefore, the selection of a more efficient contaminant transport simulator was suggested by Huang and Mayer (1997). In our opinion, this choice may partially improve the efficiency, but would not solve the problem completely. The most important way to reduce the CPU time is to reduce the repeated execution of the groundwater simulation models. The linear response matrix approach does use this concept to reduce the computational cost (Gorelick et al., 1984; Aral et al., 1995). This is also the main idea behind the PGA algorithm proposed in this study. In PGA the optimal design of groundwater remediation system was split into two stages: iteration and search stages. In the iteration stage, the simulation models were run to simplify the implicit nonlinear optimization problem into an explicit linear problem, in the neighborhood of the initial solution or any solution step thereafter. In the search stage, GA was directly applied to the linearized problem, in the neighborhood of the initial solution or any solution step thereafter. PGA uses a piecewise linearization method to reduce computational time, and simultaneously improves the computational accuracy compared to the linearization of the problem in the complete solution domain. Another important

41

advantage of PGA is, that the solution time is reduced as the iterations increase, since the solution domain shrinks as the iterations progress.

6. Conclusions An optimization model was proposed for the optimal design of pump-and-treat systems used to cleanup contaminated sites. In order to contain and extract the contaminant plume, three types of constraints were considered, which include concentration, velocity, injection and extraction balance constraints. These constraints provide standard criteria for the design of a groundwater remediation system. In the proposed model, the locations and pumping rates of extraction and injection wells were defined as continuous decision variables to have the decision space extend to 3 × mw from mw. A new combinatorial algorithm, identified as progressive genetic algorithm (PGA), was presented. PGA combines the genetic algorithm with groundwater flow and contaminant transport simulation models in an iterative manner and provides an alternative procedure for solving nonlinear optimization problems. In PGA, the solution to nonlinear optimization problem is approached by a series of linear subdomain problems. This approach overcomes the weakness of linearization of the problem in the complete solution domain. The well locations are allowed to move within selected subdomains, of which location and size changes from iteration to iteration. Subdomain optimization problems are solved by standard genetic algorithm which yields global optimal solutions within the subdomain. Since genetic algorithm is an efficient global random search method and since the selection of subdomains are not restricted in any manner and the subdomain selected, progressively moves to an optimal location, PGA yields an efficient method for obtaining optimal solution when compared to other gradient-based search methods. It must be pointed out, that PGA does not impose restrictions on the objective function, however, the constraint functions has to be continuously differentiable in order to compute the derivatives of the constraints with respect to the independent variables. In the design of pump-and-treat groundwater remediation system, this condition satisfied because the contaminant transport

42

J. Guan, M.M. Aral / Journal of Hydrology 221 (1999) 20–42

process is described by the parabolic partial differential equations. Through a design process of a pump-and-treat remediation system, the effects of pumping well locations, constraints, the number of candidate pumping wells and starting points on the optimal solution were analyzed. The locations of candidate pumping wells are very important in the optimal design of the remediation system. Defining them as continuous decision variables can significantly reduce the total discharge. The velocity constraints may be used to contain the contaminant plume within the capture zone and concentration constraints may be used to satisfy the groundwater quality standards. It is suggested that the concentration and velocity constraints should be considered in all practical engineering design, and the balance between injection and extraction rates should only be considered, if it is a design requirement. Increasing the number of candidate pumping wells may decrease the total discharge, but may increase the cost of installation and management of pumping wells. The final decision should be made by a trade-off between these factors. The starting points have no significant effect on the final solution. Numerical experiments performed in this study demonstrate that the formulation and algorithms presented in this study are reliable and robust. The computational experience we have gained in this study is encouraging and suggests that the PGA concept may be effectively utilized in the solution of this or other nonlinear optimization problems. Acknowledgements We thank Dr Vivek Kapoor of the School of Civil and Environmental Engineering, Georgia Institute of Technology for providing us with the random data used in this study and we also acknowledge the support provided by ATSDR, USDHHS in funding this project as a part of Exposure-Dose Reconstruction study.

References Ahlfeld, D.P., Sawyer, C.S., 1990. Well location in capture zone design using simulation and optimization techniques. Ground Water 28 (4), 507–512. Ahlfeld, D.P., Mulvey, J.M., Pinder, G.F., Wood, E.F., 1988.

Contaminated groundwater remediation design using simulation, optimization and sensitivity theory—1. Model development. Water Resour. Res. 24 (3), 431–441. Ahlfeld, D.P., Mulvey, J.M., Pinder, G.F., Wood, E.F., 1988. Contaminated groundwater remediation design using simulation, optimization, and sensitivity theory—2. Analysis of a field site. Water Resour. Res. 24 (3), 443–452. Ahlfeld, D.P., 1990. Two-stage groundwater remediation design. J. Water Resour. Plan. Mgmt. 116 (4), 517–529. Aral, M.M., Guan, J., 1997. Optimal groundwater remediation system design with well locations selected as decision variables, MESL-01-97, March 1997. Aral, M.M., Shea, C., Al-Kkayyal, F., 1995. Optimal Design of Pump-and-Treat Well Networks Applications of Management Science, 8, JAI Press, pp. 213–246. Aral, M.M., 1990. Ground Water Modeling in Multilayer Aquifers—Steady Flow, Lewis, London, 110pp. Aral, M.M., 1990. Ground Water Modeling in Multilayer Aquifers—Unsteady Flow, Lewis, London, 114pp. Atwood, D.F., Gorelick, S.M., 1985. Hydraulic gradient control for groundwater contaminant removal. J. Hydrol. 76, 85–106. Davis, L., 1991. Handbook of Genetic Algorithms, Van Nostrand Reinhold, New York, NY, 385pp. Dougherty, D.E., Marryott, R.A., 1991. Optimal groundwater management—1. Simulated annealing. Water Resour. Res. 27 (10), 2493–2508. Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, Reading, MA, 412pp. Gorelick, S.M., Voss, C.I., Gill, P.E., Murry, W., Saunders, M.A., Wright, M.H., 1984. Aquifer reclamation design: the use of contaminant transport simulation combined with nonlinear programming. Water Resour. Res. 20 (4), 415–427. Gorelick, S.M., 1983. A review of distributed parameter groundwater management modeling methods. Water Resour. Res. 19 (2), 305–319. Holland, J.H., 1975. Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor, MI, 183pp. Huang, C., Mayer, A.S., 1997. Pump-and-treat optimization using well locations and pumping rates as decision variables. Water Resour. Res. 33 (5), 1001–1012. Kuo, C.-W., Michel, A.N., Gray, W.G., 1992. Design of optimal pump-and-treat strategies for contaminated groundwater remediation using the simulated annealing algorithm. Adv. Water Resour. 15, 95–105. McKinney, D.C., Lin, M.D., 1994. Genetic algorithm solution of groundwater management models. Water Resour. Res. 30 (6), 1897–1906. Pinder, G.F., Gray, W.G., 1977. Finite Element Simulation in Surface and Subsurface Hydrology, Academic Press, New York, NY, 295pp. Ratzlaff, S.A., Aral, M.M., Al-Khayyal, F., 1992. Optimal design of groundwater capture systems using segmental velocity-direction constraints. Ground Water 30 (4), 607–612. Wang, W., Ahlfeld, D.P., 1994. Optimal groundwater remediation with well location as a decision variable: model development. Water Resour. Res. 30 (5), 1605–1618.