- Email: [email protected]
Optimal design of pumping networks in coastal aquifers using sharp interface models
Journal of Hydrology (2008) 361, 52– 63
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jhydrol
Optimal design of pumping networks in coastal aquifers using sharp interface models Aristotelis Mantoglou *, Maria Papantoniou Department of Rural and Surveying Engineering, Technical University of Athens, Greece Received 24 June 2008; received in revised form 11 July 2008; accepted 14 July 2008
KEYWORDS Aquifer management; Pumping network design; Coastal aquifers; Salt water intrusion; Genetic algorithms; Pumping optimization
A method for optimal design of pumping networks in coastal aquifers based on genetic algorithms and numerical solution of the governing differential equations of freshwater flow, is developed. The objective is to optimize the well locations subject to constraints that protect the aquifer from saltwater intrusion. Two are developed. The first is based on simultaneous optimization of well locations and pumping rates using genetic algorithms (GA). The second method follows two stages where the first stage concerns optimization of well locations and the second stage concerns optimization of pumping rates. The last method uses GA for searching the well positions and at each generation of GA, a simpler optimization with respect to the pumping rates is performed using nonlinear programming. The formulation of the constraints is based on numerical simulation of freshwater flow equations obtained from the sharp interface and the Ghyben–Herzberg approximations and the single potential formulation of Strack. The optimization methodologies are applied effectively for determining the optimal design of a pumping network and assessment of pumping rates in an island coastal aquifer. ª 2008 Elsevier B.V. All rights reserved.
Summary
Introduction Spatial and time distribution of freshwater often does not agree with the demand. Coastal regions and islands constitute unique cases where the water needs during the summer months are especially high due to tourism while freshwater resources are limited due to droughts. In order to meet the increasing needs for freshwater, coastal aquifers are intensively pumped causing saltwater intrusion. Therefore it is * Corresponding author. E-mail address: [email protected] (A. Mantoglou).
important to develop appropriate management models for optimal design of pumping networks and for assessing the maximum pumping rates in coastal aquifers while protecting them from saltwater intrusion. The problem of pumping optimization using aquifer simulation models has been extensively investigated in the literature (see, e.g. Mayer et al. (2002)). Gorelick (1983) and Ahlfeld and Heidari (1994) present a management approach based on groundwater simulation models and optimization methods. Many different objective functions and sets of constraints have been proposed depending on the problem. Shamir et al. (1984), Hallaji and Yazicigil (1996), Cheng
0022-1694/$ - see front matter ª 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2008.07.022
Optimal design of pumping networks in coastal aquifers using sharp interface models et al. (2000), Mantoglou (2003) and Mantoglou et al. (2004) aim at maximization of the total pumping rate while Das and Datta (1999a) aim at the minimization of the salinity of pumped water. Emch and Yeh (1998) include the pumping cost in the objective function. It is also possible to consider multiple objectives, constituting a multiobjective optimization problem (Shamir et al., 1984; Emch and Yeh, 1998; Das and Datta, 1999b). Often the purpose of optimization is to maximize the total pumping rate from a number of wells while controlling the saltwater intrusion into the aquifer. The constraints restrain the pumping rate between a minimum and a maximum value (e.g. Hallaji and Yazicigil, 1996; Emch and Yeh, 1998; Das and Datta, 1999a,b; Cheng et al., 2000). Constraints imposed by Cheng et al. (2000) and Mantoglou (2003) control the location of the toe. Additional constraints may include maintaining water levels, flow potential or salt concentration of the pumped water at desired levels. The management methods require the use of aquifer simulation models. A common approach for simulation of flow in coastal aquifers is based on the sharp interface approximation and the Ghyben–Herzberg relation (Bear, 1979). Emch and Yeh (1998), Cheng and Ouazar (1999), Cheng et al. (2000), Mantoglou (2003), and Mantoglou et al. (2004) presented coastal aquifer management models based on the sharp interface approximation. Strack’s (1976) flow potential is used in order to trace the toe of saltwater lens (Cheng and Ouazar, 1999; Mantoglou, 2003; Mantoglou et al., 2004). More complex seawater intrusion models consider the transport processes that occur in the mixing zone (Das and Datta, 1999a,b). The differential equations of flow are solved analytically or numerically in order to calculate the aquifer response for various pumping scenarios. Most research articles in the groundwater management literature focus on pumping optimization for known well locations. Gorelick et al. (1984), Wang and Ahlfeld (1994) examine optimal design of pumping networks. The optimization method employed depends on whether the optimization formulation is linear (Ahlfeld and Heidari, 1994; Hallaji and Yazicigil, 1996; Mantoglou, 2003; Mantoglou et al., 2004) or nonlinear (Gorelick et al., 1984; Shamir et al., 1984; Wang and Ahlfeld, 1994; Hallaji and Yazicigil, 1996; Emch and Yeh, 1998; Mantoglou, 2003; Mantoglou et al., 2004). The selected method also depends on the continuity of the objective function, the existence of derivatives of the objective function, the existence of local minimums, etc. Park and Aral (2004) examined the problem of optimization of well locations and pumping rates in coastal aquifers using the analytical solutions of Strack (1976) and Cheng et al. (2000). However, these analytical solutions are applicable only in aquifers of semi-infinite dimensions. Mantoglou (2003) developed analytical solutions based on the method of images that are better suited for orthogonal aquifers of finite dimensions. This paper presents a methodology for optimal design of pumping networks and assessment of pumping rates in coastal aquifers using genetic algorithms and numerical simulation governing equations of flow. The modeling methodology is based on the sharp interface approximation and the Ghyben–Herzberg relation. The flow equations are expressed using the Strack (1976) potential. The governing
53
flow equations are solved numerically using finite differences which overcomes the semi-infinite aquifer geometry limitations required by the analytical solutions of Cheng et al. (2000) and Park and Aral (2004) and the orthogonal geometry limitations of Mantoglou (2003) The management model aims at optimization of the well locations and maximization of the total pumping rate subject to constraints that protect the aquifer from saltwater intrusion. The governing differential equations are solved using finite difference method (MODFLOW). Two optimization methods are developed. The first is based on simultaneous optimization of well locations and pumping rates using genetic algorithms (GA). The second method combines genetic algorithms and nonlinear programming (Sequential Quadratic Programming – SQP) in a two stage optimization. The well position optimization is carried out using GA while at each generation of GA a simpler optimization with respect to the pumping rates is performed using SQP. Genetic algorithms are stochastic search methods that mimic the metaphor of natural biological evolution (see Holland, 1975; Goldberg, 1989; Fogel, 1994). They operate on a population of potential solutions applying the principle of survival of the fittest to produce improved approximations to a solution. A new set of approximations is created at each generation by selecting individuals according to their fitness and breeding them together using operators borrowed from natural genetics. This process leads to the evolution of populations of individuals that are better suited to their environment than their parents, just as in natural adaptation. Genetic algorithms differ substantially from more traditional search and optimization methods. The most significant difference is that GA search a population of possible solutions in parallel rather than a single solution. Additionally, GA do not require calculation of derivatives of the objective function. Therefore genetic algorithms can solve general classes of optimization problems even when the objective function is not continuous and the derivatives of the objective function do not exist. Optimization based on genetic algorithms is very flexible and, if the number of generations is large, it converges to the global rather than local minimum. The proposed optimization methodology is applied to a coastal unconfined aquifer in the Greek island of Kalymnos for optimal design of a pumping network and for assessing the pumping rates while protecting the wells from saltwater intrusion.
Governing equations of flow and saltwater intrusion in coastal aquifers Mantoglou et al. (2004) proposed a numerical aquifer model based on the sharp interface approximations and the flow potential of Strack (1976), assuming that the salt–fresh water interface encounters the aquifer base at a toe (Fig. 1). This requires that the aquifer is of large extend compared to its depth. The model is based on the sharp interface approximation which is appropriate in regional scale problems where the transition zone is narrow relative to the scale of the problem. In the steady state, the position of the interface can be estimated using the
54
A. Mantoglou, M. Papantoniou Well pumping rate Qw
Free surface for Qw = 0
S sea
fresh water interface salt water
τ
Figure 1
Coastal unconfined aquifer pumped by one well.
Ghyben–Herzberg approximation which assumes that horizontally flowing freshwater floats above static saltwater (Essaid, 1999). During transient periods, the behavior of the freshwater–saltwater interface is controlled by both the freshwater and saltwater dynamics since significant amounts of saltwater must be moved into or out of the aquifer when the interface moves. The one fluid sharp interface approach based on Ghyben–Herzberg approximation is more appropriate for modeling long term responses of freshwater lenses or short term responses in aquifers where saltwater can move in and out easily (Essaid, 1999). Fig. 1 shows a vertical cross-section of an unconfined aquifer with a sharp interface separating freshwater from saltwater where the interface intersects the aquifer base at points s (toe). It is assumed practically that the interface has been practically stabilized and is not moving very quickly. Saltwater is considered to be stagnant whereas fresh water flow is assumed practically horizontal and steady state. Variable d represents the aquifer depth from the base to the sea level, b is the total freshwater depth, n is the freshwater depth measured from the sea level and hf is the freshwater head with reference to the impermeable base of the aquifer. Ghyben–Herzberg relation, links hydraulic head hf to depth n as follows: hf d = dn, where d = (qs qf)/qf = (1.025 1)/1 = 0.025, (qs is the density of saltwater and qf is the density of freshwater). As discussed in Strack (1976) and Mantoglou et al. (2004) and outlined in Fig. 1, there exist two zones with different characteristics. In zone 1 the aquifer behaves exactly as an unconfined (phreatic) aquifer while zone 2 is the saltwater intrusion zone where a freshwater lens floats above a denser saltwater layer. Defining flow potential / of Strack (1976) / ¼ ½h2f ð1 þ dÞd 2 =2 for zone 1 and / = [(1 + d)(hf d)2]/ 2d for zone 2, the following governing equation: o o/ o o/ K þ K þNQ ¼0 ð1Þ ox ox oy oy is obtained for both zones 1 and 2 of the aquifer where K is the hydraulic conductivity (L/T), N is the surface recharge (L/T), and Q is the pumping rate (L/T). Eq. (1) is solved with respect to / and the toe of the saltwater lens is determined
by /s = [d(1 + d)d2]/2 as discussed in Strack (1976) and Mantoglou et al. (2004). Strack (1976), Cheng et al. (2000), and Park and Aral (2004) solved these equations analytically and calculated / and the location of the toe for a homogeneous aquifer of semi-infinite dimensions without surface recharge. Mantoglou (2003) derived more general analytical solutions for the case of finite aquifers of rectangular shape and when surface recharge is present. For a non-rectangular geometry of the aquifer and when the hydraulic conductivity and recharge rates are non-uniform, it is not possible to use such analytical solutions. In these cases numerical simulations are required. Notice that differential equation (1) is similar to continuity equation for a confined aquifer, where potential / corresponds to piezometric head and hydraulic conductivity K corresponds to transmissivity T. It is thus possible, with appropriate modifications, to use a standard differential equation numerical code. The well-known code Modflow (McDonald and Harbaugh, 1988), based on finite differences was chosen to solve these equations (Mantoglou, 2003; Mantoglou et al., 2004). The change of variables h ! / and T ! K, is used in order to solve Eq. (1) using MODFLOW and calculate potential /. The aquifer is discretized in a mesh of cells, the locations of which are described in terms of rows and columns. The GUI preprocessor PMWIN (Chiang and Kinzelbach, 2001) was used to prepare the necessary files for MODFLOW. The geometry and dimensions of the aquifer, the hydraulic parameters, the recharge and pumping rates, etc. are imported and saved in ASCII and binary files. These files constitute the necessary input files for MODFLOW. Mantoglou et al. (2004) examined the problem of calculating the optimum pumping rates for specified well coordinates. Here, the more complex problem of calculating the optimal position of the wells in addition to the optimal pumping rates is investigated. Two different methods of optimization are investigated. In the first method, the optimization is carried out in one stage where the well position optimization is simultaneous to pumping rate optimization. The second method follows two stages where at each cycle of the optimization of well positions a simpler optimization with respect to the pumping rates is performed.
Optimal design of pumping networks in coastal aquifers using sharp interface models
Simultaneous optimization of well locations and pumping rates using genetic algorithms
maximize :
Q tot ¼ Q 1 þ Q 2 þ þ Q k
subject to :
/i ðzÞ P 0; i ¼ 1; . . . ; k x si ðzÞ 6 x w i ; i ¼ 1; . . . ; k Q i P 0; i ¼ 1; . . . ; k
z
Q i 6 Q max;i ;
maximize :
Q tot ¼ Q 1 þ Q 2 þ þ Q k
subject to :
/i ðzÞ P 0; i ¼ 1; . . . ; k x ssi ðzÞ 6 x w i ; i ¼ 1; . . . ; k Q i P 0; i ¼ 1; . . . ; k
z
In the first method, the optimization is performed in one stage where the well positions are optimized simultaneously with the pumping rates. Let k be the number of wells pumping the coastal aquifer with pumping rates Qi; i = 1, . . . ,k. Let (xwi, ywi); i = 1, . . . ,k represent the coordinates of the wells, where xwi is the distance from the coast and ywi is the distance from the south boundary. When the well locations are not known, the objective is to determine the optimal location of the wells that yield maximum pumping rates. The optimization problem is expressed as
ð2Þ
i ¼ 1; . . . ; k
where the decision variables z = [Q1, Q2, . . . ,Qk, xw1,xw2, . . . xwk,yw1,yw2, . . . ywk] include the pumping rates and the coordinates of the wells. The first set of constraints (2) states that the value of / must be larger than 0 so that the free surface of the aquifer remains above the seawater level. These constraints can protect from saltwater intrusion only wells located far from the coast. For wells located near the coast these constraints do not provide the required protection since it is possible the toe of the interface to have advanced beyond the wells (Mantoglou et al., 2004). To avoid saltwater intrusion in this case, potential / at the well locations must be larger than the toe potential separating zone 1 from zone 2. This requirement is expressed by the second set of constraints in (2) imposed to avoid seawater intrusion, where xsi expresses the distance of the toe from the coast and xwi expresses the distance of the well from the coast respectively (Mantoglou et al., 2004). Thus the wells are protected against seawater intrusion by not allowing the toe of the interface to reach the wells. Notice that although the objective function is linear with respect to decision variables [Q1, Q2, . . . ,Qk], the optimization problem is nonlinear since the second and third set of constraints is nonlinear with respect to decision variables z. When the number of wells distributed over the aquifer is relatively large, the number of constraints in (2) is often sufficient to effectively protect the wells and the aquifer from seawater intrusion by forming a barrier where the saltwater front cannot advance deeply into the aquifer and contaminate the wells. In some cases however, when the number of wells is small or when the wells are located at a large distances from each other, the constraints in (2) may not effectively protect the aquifer. Mantoglou et al. (2004) proposed making the constraints rather stricter by modifying the second inequality in (2) as follows: xssi(z) 6 xwi, i = 1, . . . , k, where xssi corresponds to a somewhat larger value of potential than the one corresponding to the toe (see numerical examples in later sections). The optimization problem is then expressed as
55
ð3Þ
i ¼ 1; . . . ; k
Q i 6 Q max;i ;
The nonlinear constraint optimization problem (3) is solved using genetic algorithms which are stochastic search methods that mimic the metaphor of natural biological evolution (Holland, 1975; Goldberg, 1989; Fogel, 1994). They operate on a population of potential solutions applying the principle of survival of the fittest to produce improved approximations to a solution. A new set of approximations is created at each generation by selecting individuals according to their fitness and breeding them together using operators borrowed from natural genetics. This process leads to the evolution of populations of individuals that are better suited to their environment than their parents, just as in natural adaptation. Genetic algorithms model natural processes by applying the procedures of selection, recombination, mutation and migration at each iteration of the program. Since genetic algorithms solve unconstraint optimization problems, maximization problem (2) is modified in the following unconstraint form: Q tot ¼ ðQ 1 þ Q 2 þ þ Q k Þ þ pðzÞ
minimize : z
ð4Þ
where p(z) is a penalty term inserted in the optimization problem in order to handle the constraints. Penalty term p(z) is zero when the constraints are satisfied and gets a very high value when constraints are violated. In the following examples penalty term in (4) is defined as follows: pðzÞ ¼ 0;
when /i P 0 10
pðzÞ ¼ 10 ;
and x ssi 6 x w i
when /i < 0
or
xssi > x w i
ð5Þ
Fig. 2 presents the flow chart of the genetic algorithm optimization procedure. At the beginning of the computation, a number of individuals (the population) are randomly initialized and the first initial generation is produced. The objective function is then evaluated for these individuals. If the optimization criteria are not met a new generation is produced using the GA procedures of fitness assignment and parent selection, recombination, mutation and reinsertion. This cycle is performed until the termination criteria are reached. Such criteria are, for example, the number of generations or the maximum computer time. From the above discussion, it can be seen that genetic algorithms differ substantially from more traditional search and optimization methods. The most significant difference is that GA. search a population of possible solutions in parallel rather than a single solution. Additionally, GA do not require calculation of derivatives of the objective function. They require only the objective function and corresponding fitness levels which influence the direction of search. Therefore genetic algorithms can solve general classes of optimization problems even when the objective function is not continuous and the derivatives of the objective function do not exist. They are generally more straightforward to apply, because there are no restrictions for the definition of the objective function. They use probabilistic transition
56
A. Mantoglou, M. Papantoniou the second stage uses nonlinear programming methods (Sequential Quadratic Programming) as follows. The first optimization stage concerning optimization of well locations is expressed as
Initial Population of decision variables z Specification of the range of decision variables
maximize :
Evaluation of the objective function
x
k X
Q i ðxÞ
ð6Þ
i¼1
where the decision variables in (6) include the well coordinates x = [xw1, xw2, . . . ,xwk,yw1, yw2, . . . ywk]. Terms Q i ðxÞ denote the optimal pumping rates from wells i = 1, . . . k for given well locations x. Terms Q i ðxÞ are calculated from a second stage optimization with respect to pumping rates Q = [Q1, Q2, . . ., Qk] for fixed values of x. For known location of the wells x, the following objective function is maximized with respect to decision variables Q = [Q1, Q2, . . ., Qk] subject to constraints that limit saltwater intrusion into the aquifer
Selection Recombination Mutation Reinsertion
New Generation of z
maximize :
Q tot ¼ Q 1 þ Q 2 þ þ Q k
subject to :
/i ðQÞ P 0 i ¼ 1; . . . ; k x ssi ðQÞ 6 x w i ; i ¼ 1; . . . ; k
Q
Evaluation of the objective function
Are termination criteria met?
Q tot ¼
NO
YES Optimal values of well coordinates and pumping rates
Figure 2 Optimization flow chart for well locations using one stage optimization with genetic algorithms.
rules and can produce a number of potential solutions to a given problem. Thus, in cases where the particular problem does not have one individual optimal solution, the genetic algorithms are potentially useful for identifying these alternative solutions simultaneously. Optimization based on genetic algorithms is very flexible and, if the population size and the number of generations are sufficient, it converges to the global rather than a local optimum. However, it is often slow converging and it requires a large number of iterations. Convergence is improved by choosing reasonable initial population values. Increasing the number of individuals in the population and the number of generations improves the results but it requires more computer time. The genetic algorithm optimization program needed for calculation of the optimal well locations and pumping rates, was developed in MATLAB. MODFLOW was used for solving differential equation (1) and evaluating the required potential at each iteration of optimization.
Two stage optimization of well locations and pumping rates The second method of optimization follows two stages where the first stage concerns optimization of the location of wells and the second stage concerns optimization of pumping rates. The first stage uses genetic algorithms and
ð7Þ
Terms Q ðxÞ ¼ ½Q 1 ; Q 2 ; Q k in (6) represent the optimal pumping rates calculated from (7) for given values of x. Maximization of (6) with respect to x, is performed using genetic algorithms while the optimal values Q ðxÞ ¼ ½Q 1 ; Q 2 ; . . . Q k for each value of x are determined from maximization of (7) using nonlinear optimization methods. Since the decision variables in (7) are real and the objective function is continuous and differentiable with respect to the decision variables, nonlinear optimization method based on Sequential Quadratic Programming (SQP) is used for this second stage minimization (see Mantoglou et al. (2004)). Fig. 3 outlines the flow diagram of the two stage maximization based on genetic algorithms and Sequential Quadratic Programming. The range of the well coordinates and pumping rates are specified and an initial population [xini,1, xini,2, . . . xini,N] of N sets of well coordinates is generated. Given the well coordinate sets xini,m, the optimal pumping rates corresponding to xini,m are computed from (7) using SQP. Then a new population of well coordinates is generated based on the fittest members of the population using the selection, recombination, mutation, reinsertion and migration procedures of the genetic algorithm. When the termination criteria specified by the user are met, the GA optimization procedure ends. Such termination criteria are, for example, the number of generations or the computer time. The optimal solution corresponds to a specific set of coordinates which provide the maximum total pumping rate. As stated before, Sequential Quadratic Programming (SQP) method is used for calculation of the optimal pumping rates for specified well locations x. SQP is based on an iterative procedure where at each major iteration a positive definite quasi-Newton approximation of the Hessian of the Lagrangian function is calculated using the BFGS method (Goldfarb, 1970; Shanno, 1970; Fletcher, 1987; Luenberger, 1989). The optimization can be sensitive to numerical approximations such as truncation and roundoff error in the calculation of finite-difference gradients.
Optimal design of pumping networks in coastal aquifers using sharp interface models
57
Initial population of well coordinates x ini Specification of range of well coordinates
x
Initial values of well pumping rates
New Generation of coordinates
Specification of pumping rates
S.Q.P.
Reinsertion
Are termination criteria of S.Q.P. optimization met?
Mutation
NO
Recombination
YES
Selection
Optimal pumping Q* ( x ) for the specific coordinates Q* ( x )
Are termination criteria of GA. optimization met?
NO
YES Optimal well coordinates and pumping rates
Figure 3
Flow diagram for two-stage optimization using genetic algorithms and Sequential Quadratic Programming.
Furthermore, the program may converge to a local minimum rather than the global minimum. Starting the optimization with different initial conditions may help to locate the global minimum. Additionally, by changing the values of the parameters that take part in finite differences calculation or increasing the number of iterations contributes to the computation of the global minimum. Obviously, the above considerations have an important effect on the effectiveness of the two stage method that combines GA and SQP. The genetic algorithm and Sequential Quadratic Programming optimization programs needed for calculation of the optimum well coordinates and pumping rates were developed in MATLAB. MODFLOW was used to solve differential equation (1) and to evaluate the flow potential required at each iteration of optimization.
Application to an aquifer of orthogonal shape An aquifer having a simple rectangular geometry (Fig. 4) is selected in order to test and compare the two optimization methods. The recharge rate at the higher elevations of the aquifer is N = 150 mm/year distributed over an area of A = 9 km2 while the recharge rate at the lower elevations of the aquifer is N = 30 mm/day. The sea boundary constitutes a constant head boundary whereas the north and south boundaries are assumed impermeable. The aquifer thick-
ness from the base to the sea surface is taken as d = 25 m. The potential at the toe of seawater lens is calculated according to equation /s ¼ dð1þdÞ d 2 ¼ 8:0078 m2 . The geom2 etry of this aquifer resembles the aquifer at Vathi in the Greek island of Kalymnos described in detail in Mantoglou et al. (2004). The toe potential in this case is equal to /s ¼ 8:0078 m2 . When the number of wells distributed over the aquifer is relatively large, the constraints in (2) can effectively protect the wells from seawater intrusion by forming a barrier where the saltwater front cannot advance deeply into the aquifer and contaminate the wells. In some cases however, when the number of wells is small or when the wells lie at a large distance from each other, the constraints in (2) might not effectively protect the aquifer. Mantoglou et al. (2004) proposed to modify the second Mantoglou inequality in (2) and make the constraints somewhat stricter as follows: x ssi ðzÞ 6 xw i ; i ¼ 1; . . . ; k where x ssi corresponds to a somewhat larger value of potential than that corresponding to the toe. In the examples presented below, that setting /ss ¼ 8:1 m2 rather than /s ¼ 8:0078 m2 , make constraints x ssi ðzÞ 6 xw i ; i ¼ 1; . . . ; k stricter and effectively protect the wells. As is shown in the figures bellow the suggested sets of constraints allow the wells to yield the maximum pumping rates while protecting them sufficiently from seawater intrusion. More conservative solutions can be obtained by increasing /ss even more.
58
A. Mantoglou, M. Papantoniou
Figure 4
Table 1
GA GA & SQP
A test aquifer of rectangular shape.
Optimization results for three pumping wells for two optimization methods x1 (m)
y1 (m)
Q1 (m3/day)
x2 (m)
y2 (m)
Q2 (m3/day)
x3 (m)
y3 (m)
Q3 (m3/day)
Qtot (m3/day)
6350 5950
1950 2050
1378.3 1513.9
5550 6350
1250 1250
1327.6 1442.1
6350 6050
1450 1150
1604.2 1363.6
4310.1 4319.6
In a first application, the orthogonal aquifer of Fig. 4 with uniform hydraulic conductivity and recharge rate is selected in order to calculate the optimal coordinates of three wells and the optimal pumping rates. The results of optimization obtained from both methods of optimization are presented in Table 1. The first line of Table 1 presents the well coordinates and the optimal pumping rates computed by the first method based on GA while the second line presents the corresponding results obtained from the two stage optimization using GA and SQP. Both optimization methods yield a similar total pumping rate but the individual pumping rates and well
4
3000
16 12
8
2000 4 8 8
12
4
0 0
1000
2000
3000
4000
5000
16 20
8
8
1000
12 16
8
6000
7000
Figure 5 Equipotentials for the optimal well locations and pumping rates according to one stage optimization method.
12
3000 4
16
8 8
2000 12 8
4
8
16
8
0
1000
16
0
8
4
12 8
1000
2000
3000
4000
5000
6000
7000
Figure 6 Equipotentials for the optimal well locations and pumping rates for two stage methodology.
locations are different. Fig. 5 plots the equipotentials corresponding to the calculated optimal values with GA method while Fig. 6 plots the equipotentials corresponding to the optimal well coordinates and pumping rates evaluated using the two stage optimization methodology. While the two stage methodology yields a somewhat better Qtot than the one stage methodology, the computer time required is significantly larger.
Application to coastal aquifer in a Greek island The simulation and optimization methodology developed above is applied to an unconfined aquifer located at Vathi valley in the Greek island of Kalymnos (Fig. 7). As discussed in Mantoglou et al. (2004), the aquifer consists of high permeability limestone which outcrops at the highland areas while the impermeable bottom of the aquifer is composed by schist. The geologic formations found in the valley are highly permeable scree, medium permeability alluvium deposits and almost impermeable tuff (volcanic formations). Although in parts of the aquifer flow is taking place in the limestone fissures, we suppose for the purposes of this analysis that the aquifer can be represented by an equivalent porous medium. Based on existing wells the thickness of the aquifer is taken as: d = 25 m. According to Mantoglou et al. (2004) there are four zones with different surface recharge rates. The large recharge rate N = 150 mm/year in the limestone parts of the aquifer is justified by the large crevices present in the limestone facilitating percolation of surface water. The recharge rate at the scree area is smaller and is estimated at N = 70 mm/ year while the presence of clay in the alluvium deposits reduces the recharge rate to N = 30 mm/year. The tuffs are regarded as impermeable formations with zero recharge. Based on the geology of the valley, the aquifer is divided into four zones having different hydraulic conductivities (Fig. 8). The zonal hydraulic conductivities were evaluated by model calibration so that the model (1) predicts the measured fresh water levels in the existing wells (see Mantoglou et al. (2004)).
Optimal design of pumping networks in coastal aquifers using sharp interface models
Figure 7
Figure 8
59
Kalymnos aquifer.
The four aquifer zones with different hydraulic conductivity.
north boundary. Based on the above, the west, south and north boundaries are considered impermeable in the model. In the optimization case study, the aquifer is pumped by 11 wells and is replenished by surface recharge. Besides the pumping rates, the well locations are calculated using the optimization methodology. Table 2 presents the optimal pumping rates obtained by the one stage methodology while Figs. 9 and 10 plot the cor-
The sea along the east of the aquifer creates a constant hydraulic head boundary. The west and south boundaries of the aquifer correspond to groundwater divides, estimated based on hydrogeologic maps of the region, and are considered as impermeable boundaries in the model. Although the sea is very close to the aquifer in the north boundary, the existing wells are not contaminated by seawater which justifies the assumption of an impermeable boundary at the
Table 2 Optimal pumping rates (m3/day) for the optimal pumping network with 11 wells obtained from the one stage simultaneous optimization method with GA Well
1
2
3
4
5
6
7
8
9
10
11
Sum
Pumping rate
404.1
427.2
486.9
539.6
774.6
617.6
202.8
354.3
213.6
608.5
1101.2
5730.4
60
A. Mantoglou, M. Papantoniou
44
52
2 4
60 52
7 8 32
2000
40 36 2 3 5
28 24
3000
60
96 92
84 76 68
1
84 76 68
4000
92
5000
28
The above results indicate that both methods produce similar results regarding the optimal well location and optimal pumping rates (QGA = 5730.4 m3/day, QGA,SQP = 5639.7 m3/day, respectively). The two stage optimization method (using SQP and GA) requires more computer time than the one stage optimization. The tolerance level of SQP parameters in the two stage optimization was set to relatively coarse values so that the required iterations and the computation time are not excessively large. Decreasing the tolerance level and increasing the number of iterations of SQP, while increasing the population size and the number of generations of GA would im-
responding equipotential and the piezometric head maps. The optimal location of the wells is also shown in these figures. Table 3 plots the optimal pumping rates for the optimal pumping network with 11 wells obtained using the two stage optimization methodology which combines GA and SQP. The population size of the GA is equal to 30 and the number of generations is 50. Figs. 11 and 12 plot the optimal location of wells and the corresponding equipotential and the piezometric head maps, respectively. Fig. 13 plots the total maximum pumping rate obtained as a function of generation number evaluated from the genetic algorithm.
24 20 6 61
3 16
12
9
84
10 8
11
4
8
1000
1000
3000
4000
5000
6000
8000
9000
2. 62 5
1. 62 5 1. 37 5 1. 12 5 0.75 10
11
2000
3000
4000
5000
6000
0.5
1000
1000
0. 875
0.625
9
4
7000
0. 375 0.2 5 0.12 5
2000
6
1.375 1.1 25 0. 87 5
1.875
7
1.625
3.2 5 3
25 5 2. 87 2 1. 3 5 8
1 3.25 3 2.6 25 2. 25
3000
3.625
4000
3.6 25
5000
Figure 10
7000
Equipotentials for the optimal pumping network with 11 wells obtained from one stage method using GA.
3.8 75
Figure 9
2000
8000
9000
Piezometric head map for the optimal pumping network with 11 wells obtained from one stage method using GA.
Table 3 Optimal pumping rates (m3/day) for the optimal pumping network with 11 wells obtained from the two stage method that combines GA and SQP Well
1
2
3
4
5
6
7
8
9
10
11
SUM
Pumping rate
546.6
533.1
547.2
518.5
547.5
547.2
549.4
171.1
552.4
565.2
561.6
5639.7
Optimal design of pumping networks in coastal aquifers using sharp interface models
61
68 60 52
44 0 1 4 36 32
2
2000
24
3
5
28
24 20 16
4 6
28 24 2 1 0 12 6
44 36
3000
68 60 52
76
4000
80
80 76
5000
48 10
9
12
12
8
7
11 4
8
1000
1000
2000
3000
4000
5000
6000
7000
8000
9000
9
10
3000
4000
5000
6000
7000
0.5
2000
11 0.625
1000
1000
5 62 0.
8000
0.25 0.12 5
7
8
1.25 1.12 5 1 0.875 0.75
0.375
4 6
0.6 25
2000
75 1.3
1
3000
5 .87 1 1 25 1.6 2 3 5
0. 75
3
4000
1.25
5000
3 2.625 2 2.25 .62 5 1.87 5 2. 1.625 25
Equipotentials for the optimal pumping network with 11 wells obtained from the two stage method that combines GA
3. 3 3.2 75 5
Figure 11 and SQP.
4
9000
Figure 12 Piezometric head map for the optimal pumping network with 11 wells obtained from the two stage method that combines GA and SQP.
prove the solution obtained by the two stage methodology and it could outperform the one stage method at the cost of more computer time. Mantoglou et al. (2004) assessed the optimum pumping rates in the same aquifer from a set of 11 pumping wells with specified fixed coordinates from a non-optimal pumping network and found a maximum total pumping rate of about Qtot = 5360 m3/day. Notice that the total pumping rate obtained with fixed well locations from a non-optimal pumping network, is quite smaller than the pumping rate obtained here from the optimal pumping network where the well locations are included in the decision variables.
Summary and conclusions A method for optimal design of pumping networks in coastal aquifers was developed and utilized in a coastal aquifer in the Greek island of Kalymnos. The objective is to optimize the well locations and to maximize the total pumping rates
subject to constraints that protect the aquifer from saltwater intrusion. The method is based on genetic algorithms and nonlinear optimization subject to constraints that limit the saltwater intrusion into the aquifer. The simulation model is based on the sharp interface and the Ghyben–Herzberg approximations. The single potential formulation of Stack (Strack (1976) was used and the governing equations were solved numerically using finite differences (Mantoglou et al., 2004). The numerical model is more general than previous analytical solutions (Cheng et al., 2000; Mantoglou, 2003; Park and Aral, 2004) and it can handle aquifers of complex shapes, non-uniform hydraulic conductivity, non-uniform distribution of surface recharge, etc. The governing differential equations are solved using MODFLOW finite difference code. Two optimization methods are developed. The first method is based on simultaneous optimization of well locations and pumping rates using genetic algorithms
62
A. Mantoglou, M. Papantoniou 5650
5600
3
Qtotal (m /day)
5550
5500
5450
5400
5350 0
10
20
30
40
50
generation number Figure 13 to 30.
Maximum total pumping rate as a function of generation number for 50 generations. The population size of GA is equal
(GA). The second method combines genetic algorithms and Sequential Quadratic Programming in a two stage optimization. The well position optimization in the last method is carried out using GA while at each generation of GA a simpler optimization with respect to the pumping rates is performed using nonlinear optimization (Sequential Quadratic Programming). The evaluation of the constraints at each iteration is based on the numerical simulation. The simulation and optimization methodology was applied to an unconfined aquifer located at the valley of Vathi in the Greek island of Kalymnos for optimal design of a pumping network and for assessing the maximum pumping rates while protecting the wells from saltwater intrusion. The aquifer is modeled as heterogeneous with four zones of different hydraulic conductivity and four areas with different recharge rates. Using the two optimization methods developed in the paper, the optimal well locations and the maximum pumping rates were assessed. The methodology for optimal design of pumping networks based on a two stage optimization (first stage concerns well location optimization and second stage concerns pumping optimization) yields similar results with the one stage simultaneous optimization of well location and pumping rates but it requires more computer time. Decreasing the tolerance level and increasing the number of iterations of SQP while increasing the population size and the number of generations of GA would improve the solution of the two stage methodology and is expected to outperform the one stage method at the cost of more computer time. The pumping rates obtained from both methodologies are quite higher than the optimal pumping rate obtained from the pumping optimization method of Mantoglou
et al. (2004) based on pumping optimization with fixed non-optimal well locations. The proposed methodology, based on the sharp interface approximation, is applicable to regional scale problems where the transition zone is narrow relative to the scale of the problem. The proposed numerical solution of the governing differential equation does not have the limitations required by analytical solutions regarding orthogonal shape of the aquifer and homogeneity of hydraulic conductivity (Cheng et al., 2000; Mantoglou, 2003; Park and Aral, 2004) and can solve the problem of optimal well design in cases where the shape of the aquifer is of irregular geometry or when the hydraulic conductivity of the aquifer and the recharge rate is nonuniform.
References Ahlfeld, D.P., Heidari, M., 1994. Applications of optimal hydraulic control to ground-water systems. Journal of Water Resources Planning and Management 120 (3), 350–365. Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York. Cheng, A.H.-D., Ouazar, D., 1999. Analytical solutions. In: Bear, J., Cheng, A.H.-D., Sorek, S., Ouazar, D., Herrera, I. (Eds.), Seawater Intrusion in Coastal Aquifers – Concepts, Methods and Practices. Kluwer Academic Publishers. Cheng, A.H.-D., Halhal, D., Naji, A., Ouazar, D., 2000. Pumping optimization in saltwater-intruded coastal aquifers. Water Resources Research 36 (8), 2155–2165. Chiang, W-H., Kinzelbach, W., 2001. 3D-Groundwater Modeling with PMWIN: A Simulation System for Modeling Groundwater Flow and Pollution. Springer-Verlag, Berlin, Heidelberg, New York. Das, A., Datta, B., 1999a. Development of multiobjective management models for coastal aquifers. Journal of Water Resources Planning and Management 125 (2), 76–87.
Optimal design of pumping networks in coastal aquifers using sharp interface models Das, A., Datta, B., 1999b. Development of management models for sustainable use of coastal aquifers. Journal Irrigation and Drainage Engineering 125 (3), 112–121. Emch, P.G., Yeh, W.W.-G., 1998. Management model for conjunctive use of coastal surface water and ground water. Journal Water Resources Planning and Management 124 (3), 129–139. Essaid, H.I., 1999. USGS SHARP model. In: Bear, J., Cheng, A.H.-D., Sorek, S., Ouazar, D., Herrera, I. (Eds.), Seawater Intrusion in Coastal Aquifers – Concepts, Methods and Practices. Kluwer Academic Publishers. Fletcher, R., 1987. Practical Methods of Optimization. John Wiley & Sons. Fogel, D.B., 1994. An introduction to simulated evolutionary optimization. IEEE Transactions on Neural Networks 5 (1), 3–14. Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization and Machine learning. Addison-Wesley, Reading, MA. Goldfarb, D., 1970. A family of variable – metric methods derived by variational means. Mathematics of Computing 24, 23–26. Gorelick, S.M., 1983. A review of distributed parameter groundwater management modeling methods. Water Resources Research 19 (2), 305–319. Gorelick, S.M., Voss, C.I., Gill, P.E., Murray, W., Saunders, M.A., Wright, M.H., 1984. Aquifer reclamation design: the use of contaminant transport simulation combined with nonlinear programming. Water Resources Research 20 (4), 415–427. Hallaji, K., Yazicigil, H., 1996. Optimal management of a coastal aquifer in Southern Turkey. Journal of Water Resources Planning and Management 122 (4), 233–244. Holland, J.H., 1975. Adaptation in Natural and Artificial Systems. The University of Michigan Press.
63
Luenberger, D.G., 1989. Linear and Nonlinear Programming. Addison-Wesley Publishing Company. Mantoglou, A., 2003. Pumping management of coastal aquifers using analytical models of saltwater intrusion. Water Recourses Research 39 (12), 1–12. Mantoglou, A., Papantoniou, M., Giannoulopoulos, P., 2004. Management of coastal aquifers based on nonlinear optimization and evolutionary algorithms. Journal of Hydrology 297, 209–228. Mayer, A.S., Kelley, C.T., Miller, C.T., 2002. Optimal design for problems involving flow and transport phenomena in saturated subsurface systems. Advances in Water Resources 25, 1233– 1256. McDonald, M.G., Harbaugh, A.W., 1988. A modular three-dimensional finite-difference ground-water flow model. Technical Report, U.S. Geol. Survey, Reston, VA. Park, C.-H., Aral, M.M., 2004. Multi-objective optimization of pumping rates and well placement in coastal aquifers. Journal of Hydrology 290, 80–99. Shamir, U., Bear, J., Gamliel, A., 1984. Optimal annual operation of a coastal aquifer. Water Resources Research 20 (4), 435–444. Shanno, D.F., 1970. Conditioning of quasi-Newton methods for function minimization. Mathematics of Computation 24, 647–656. Strack, O.D.L., 1976. A single-potential solution for regional interface problems in coastal aquifers. Water Resources Research 12 (6), 1165–1174. Wang, W., Ahlfeld, D.P., 1994. Optimal groundwater remediation with well location as a decision variable: model development. Water Resources Research 30 (5), 1605–1618.
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jhydrol
Optimal design of pumping networks in coastal aquifers using sharp interface models Aristotelis Mantoglou *, Maria Papantoniou Department of Rural and Surveying Engineering, Technical University of Athens, Greece Received 24 June 2008; received in revised form 11 July 2008; accepted 14 July 2008
KEYWORDS Aquifer management; Pumping network design; Coastal aquifers; Salt water intrusion; Genetic algorithms; Pumping optimization
A method for optimal design of pumping networks in coastal aquifers based on genetic algorithms and numerical solution of the governing differential equations of freshwater flow, is developed. The objective is to optimize the well locations subject to constraints that protect the aquifer from saltwater intrusion. Two are developed. The first is based on simultaneous optimization of well locations and pumping rates using genetic algorithms (GA). The second method follows two stages where the first stage concerns optimization of well locations and the second stage concerns optimization of pumping rates. The last method uses GA for searching the well positions and at each generation of GA, a simpler optimization with respect to the pumping rates is performed using nonlinear programming. The formulation of the constraints is based on numerical simulation of freshwater flow equations obtained from the sharp interface and the Ghyben–Herzberg approximations and the single potential formulation of Strack. The optimization methodologies are applied effectively for determining the optimal design of a pumping network and assessment of pumping rates in an island coastal aquifer. ª 2008 Elsevier B.V. All rights reserved.
Summary
Introduction Spatial and time distribution of freshwater often does not agree with the demand. Coastal regions and islands constitute unique cases where the water needs during the summer months are especially high due to tourism while freshwater resources are limited due to droughts. In order to meet the increasing needs for freshwater, coastal aquifers are intensively pumped causing saltwater intrusion. Therefore it is * Corresponding author. E-mail address: [email protected] (A. Mantoglou).
important to develop appropriate management models for optimal design of pumping networks and for assessing the maximum pumping rates in coastal aquifers while protecting them from saltwater intrusion. The problem of pumping optimization using aquifer simulation models has been extensively investigated in the literature (see, e.g. Mayer et al. (2002)). Gorelick (1983) and Ahlfeld and Heidari (1994) present a management approach based on groundwater simulation models and optimization methods. Many different objective functions and sets of constraints have been proposed depending on the problem. Shamir et al. (1984), Hallaji and Yazicigil (1996), Cheng
0022-1694/$ - see front matter ª 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2008.07.022
Optimal design of pumping networks in coastal aquifers using sharp interface models et al. (2000), Mantoglou (2003) and Mantoglou et al. (2004) aim at maximization of the total pumping rate while Das and Datta (1999a) aim at the minimization of the salinity of pumped water. Emch and Yeh (1998) include the pumping cost in the objective function. It is also possible to consider multiple objectives, constituting a multiobjective optimization problem (Shamir et al., 1984; Emch and Yeh, 1998; Das and Datta, 1999b). Often the purpose of optimization is to maximize the total pumping rate from a number of wells while controlling the saltwater intrusion into the aquifer. The constraints restrain the pumping rate between a minimum and a maximum value (e.g. Hallaji and Yazicigil, 1996; Emch and Yeh, 1998; Das and Datta, 1999a,b; Cheng et al., 2000). Constraints imposed by Cheng et al. (2000) and Mantoglou (2003) control the location of the toe. Additional constraints may include maintaining water levels, flow potential or salt concentration of the pumped water at desired levels. The management methods require the use of aquifer simulation models. A common approach for simulation of flow in coastal aquifers is based on the sharp interface approximation and the Ghyben–Herzberg relation (Bear, 1979). Emch and Yeh (1998), Cheng and Ouazar (1999), Cheng et al. (2000), Mantoglou (2003), and Mantoglou et al. (2004) presented coastal aquifer management models based on the sharp interface approximation. Strack’s (1976) flow potential is used in order to trace the toe of saltwater lens (Cheng and Ouazar, 1999; Mantoglou, 2003; Mantoglou et al., 2004). More complex seawater intrusion models consider the transport processes that occur in the mixing zone (Das and Datta, 1999a,b). The differential equations of flow are solved analytically or numerically in order to calculate the aquifer response for various pumping scenarios. Most research articles in the groundwater management literature focus on pumping optimization for known well locations. Gorelick et al. (1984), Wang and Ahlfeld (1994) examine optimal design of pumping networks. The optimization method employed depends on whether the optimization formulation is linear (Ahlfeld and Heidari, 1994; Hallaji and Yazicigil, 1996; Mantoglou, 2003; Mantoglou et al., 2004) or nonlinear (Gorelick et al., 1984; Shamir et al., 1984; Wang and Ahlfeld, 1994; Hallaji and Yazicigil, 1996; Emch and Yeh, 1998; Mantoglou, 2003; Mantoglou et al., 2004). The selected method also depends on the continuity of the objective function, the existence of derivatives of the objective function, the existence of local minimums, etc. Park and Aral (2004) examined the problem of optimization of well locations and pumping rates in coastal aquifers using the analytical solutions of Strack (1976) and Cheng et al. (2000). However, these analytical solutions are applicable only in aquifers of semi-infinite dimensions. Mantoglou (2003) developed analytical solutions based on the method of images that are better suited for orthogonal aquifers of finite dimensions. This paper presents a methodology for optimal design of pumping networks and assessment of pumping rates in coastal aquifers using genetic algorithms and numerical simulation governing equations of flow. The modeling methodology is based on the sharp interface approximation and the Ghyben–Herzberg relation. The flow equations are expressed using the Strack (1976) potential. The governing
53
flow equations are solved numerically using finite differences which overcomes the semi-infinite aquifer geometry limitations required by the analytical solutions of Cheng et al. (2000) and Park and Aral (2004) and the orthogonal geometry limitations of Mantoglou (2003) The management model aims at optimization of the well locations and maximization of the total pumping rate subject to constraints that protect the aquifer from saltwater intrusion. The governing differential equations are solved using finite difference method (MODFLOW). Two optimization methods are developed. The first is based on simultaneous optimization of well locations and pumping rates using genetic algorithms (GA). The second method combines genetic algorithms and nonlinear programming (Sequential Quadratic Programming – SQP) in a two stage optimization. The well position optimization is carried out using GA while at each generation of GA a simpler optimization with respect to the pumping rates is performed using SQP. Genetic algorithms are stochastic search methods that mimic the metaphor of natural biological evolution (see Holland, 1975; Goldberg, 1989; Fogel, 1994). They operate on a population of potential solutions applying the principle of survival of the fittest to produce improved approximations to a solution. A new set of approximations is created at each generation by selecting individuals according to their fitness and breeding them together using operators borrowed from natural genetics. This process leads to the evolution of populations of individuals that are better suited to their environment than their parents, just as in natural adaptation. Genetic algorithms differ substantially from more traditional search and optimization methods. The most significant difference is that GA search a population of possible solutions in parallel rather than a single solution. Additionally, GA do not require calculation of derivatives of the objective function. Therefore genetic algorithms can solve general classes of optimization problems even when the objective function is not continuous and the derivatives of the objective function do not exist. Optimization based on genetic algorithms is very flexible and, if the number of generations is large, it converges to the global rather than local minimum. The proposed optimization methodology is applied to a coastal unconfined aquifer in the Greek island of Kalymnos for optimal design of a pumping network and for assessing the pumping rates while protecting the wells from saltwater intrusion.
Governing equations of flow and saltwater intrusion in coastal aquifers Mantoglou et al. (2004) proposed a numerical aquifer model based on the sharp interface approximations and the flow potential of Strack (1976), assuming that the salt–fresh water interface encounters the aquifer base at a toe (Fig. 1). This requires that the aquifer is of large extend compared to its depth. The model is based on the sharp interface approximation which is appropriate in regional scale problems where the transition zone is narrow relative to the scale of the problem. In the steady state, the position of the interface can be estimated using the
54
A. Mantoglou, M. Papantoniou Well pumping rate Qw
Free surface for Qw = 0
S sea
fresh water interface salt water
τ
Figure 1
Coastal unconfined aquifer pumped by one well.
Ghyben–Herzberg approximation which assumes that horizontally flowing freshwater floats above static saltwater (Essaid, 1999). During transient periods, the behavior of the freshwater–saltwater interface is controlled by both the freshwater and saltwater dynamics since significant amounts of saltwater must be moved into or out of the aquifer when the interface moves. The one fluid sharp interface approach based on Ghyben–Herzberg approximation is more appropriate for modeling long term responses of freshwater lenses or short term responses in aquifers where saltwater can move in and out easily (Essaid, 1999). Fig. 1 shows a vertical cross-section of an unconfined aquifer with a sharp interface separating freshwater from saltwater where the interface intersects the aquifer base at points s (toe). It is assumed practically that the interface has been practically stabilized and is not moving very quickly. Saltwater is considered to be stagnant whereas fresh water flow is assumed practically horizontal and steady state. Variable d represents the aquifer depth from the base to the sea level, b is the total freshwater depth, n is the freshwater depth measured from the sea level and hf is the freshwater head with reference to the impermeable base of the aquifer. Ghyben–Herzberg relation, links hydraulic head hf to depth n as follows: hf d = dn, where d = (qs qf)/qf = (1.025 1)/1 = 0.025, (qs is the density of saltwater and qf is the density of freshwater). As discussed in Strack (1976) and Mantoglou et al. (2004) and outlined in Fig. 1, there exist two zones with different characteristics. In zone 1 the aquifer behaves exactly as an unconfined (phreatic) aquifer while zone 2 is the saltwater intrusion zone where a freshwater lens floats above a denser saltwater layer. Defining flow potential / of Strack (1976) / ¼ ½h2f ð1 þ dÞd 2 =2 for zone 1 and / = [(1 + d)(hf d)2]/ 2d for zone 2, the following governing equation: o o/ o o/ K þ K þNQ ¼0 ð1Þ ox ox oy oy is obtained for both zones 1 and 2 of the aquifer where K is the hydraulic conductivity (L/T), N is the surface recharge (L/T), and Q is the pumping rate (L/T). Eq. (1) is solved with respect to / and the toe of the saltwater lens is determined
by /s = [d(1 + d)d2]/2 as discussed in Strack (1976) and Mantoglou et al. (2004). Strack (1976), Cheng et al. (2000), and Park and Aral (2004) solved these equations analytically and calculated / and the location of the toe for a homogeneous aquifer of semi-infinite dimensions without surface recharge. Mantoglou (2003) derived more general analytical solutions for the case of finite aquifers of rectangular shape and when surface recharge is present. For a non-rectangular geometry of the aquifer and when the hydraulic conductivity and recharge rates are non-uniform, it is not possible to use such analytical solutions. In these cases numerical simulations are required. Notice that differential equation (1) is similar to continuity equation for a confined aquifer, where potential / corresponds to piezometric head and hydraulic conductivity K corresponds to transmissivity T. It is thus possible, with appropriate modifications, to use a standard differential equation numerical code. The well-known code Modflow (McDonald and Harbaugh, 1988), based on finite differences was chosen to solve these equations (Mantoglou, 2003; Mantoglou et al., 2004). The change of variables h ! / and T ! K, is used in order to solve Eq. (1) using MODFLOW and calculate potential /. The aquifer is discretized in a mesh of cells, the locations of which are described in terms of rows and columns. The GUI preprocessor PMWIN (Chiang and Kinzelbach, 2001) was used to prepare the necessary files for MODFLOW. The geometry and dimensions of the aquifer, the hydraulic parameters, the recharge and pumping rates, etc. are imported and saved in ASCII and binary files. These files constitute the necessary input files for MODFLOW. Mantoglou et al. (2004) examined the problem of calculating the optimum pumping rates for specified well coordinates. Here, the more complex problem of calculating the optimal position of the wells in addition to the optimal pumping rates is investigated. Two different methods of optimization are investigated. In the first method, the optimization is carried out in one stage where the well position optimization is simultaneous to pumping rate optimization. The second method follows two stages where at each cycle of the optimization of well positions a simpler optimization with respect to the pumping rates is performed.
Optimal design of pumping networks in coastal aquifers using sharp interface models
Simultaneous optimization of well locations and pumping rates using genetic algorithms
maximize :
Q tot ¼ Q 1 þ Q 2 þ þ Q k
subject to :
/i ðzÞ P 0; i ¼ 1; . . . ; k x si ðzÞ 6 x w i ; i ¼ 1; . . . ; k Q i P 0; i ¼ 1; . . . ; k
z
Q i 6 Q max;i ;
maximize :
Q tot ¼ Q 1 þ Q 2 þ þ Q k
subject to :
/i ðzÞ P 0; i ¼ 1; . . . ; k x ssi ðzÞ 6 x w i ; i ¼ 1; . . . ; k Q i P 0; i ¼ 1; . . . ; k
z
In the first method, the optimization is performed in one stage where the well positions are optimized simultaneously with the pumping rates. Let k be the number of wells pumping the coastal aquifer with pumping rates Qi; i = 1, . . . ,k. Let (xwi, ywi); i = 1, . . . ,k represent the coordinates of the wells, where xwi is the distance from the coast and ywi is the distance from the south boundary. When the well locations are not known, the objective is to determine the optimal location of the wells that yield maximum pumping rates. The optimization problem is expressed as
ð2Þ
i ¼ 1; . . . ; k
where the decision variables z = [Q1, Q2, . . . ,Qk, xw1,xw2, . . . xwk,yw1,yw2, . . . ywk] include the pumping rates and the coordinates of the wells. The first set of constraints (2) states that the value of / must be larger than 0 so that the free surface of the aquifer remains above the seawater level. These constraints can protect from saltwater intrusion only wells located far from the coast. For wells located near the coast these constraints do not provide the required protection since it is possible the toe of the interface to have advanced beyond the wells (Mantoglou et al., 2004). To avoid saltwater intrusion in this case, potential / at the well locations must be larger than the toe potential separating zone 1 from zone 2. This requirement is expressed by the second set of constraints in (2) imposed to avoid seawater intrusion, where xsi expresses the distance of the toe from the coast and xwi expresses the distance of the well from the coast respectively (Mantoglou et al., 2004). Thus the wells are protected against seawater intrusion by not allowing the toe of the interface to reach the wells. Notice that although the objective function is linear with respect to decision variables [Q1, Q2, . . . ,Qk], the optimization problem is nonlinear since the second and third set of constraints is nonlinear with respect to decision variables z. When the number of wells distributed over the aquifer is relatively large, the number of constraints in (2) is often sufficient to effectively protect the wells and the aquifer from seawater intrusion by forming a barrier where the saltwater front cannot advance deeply into the aquifer and contaminate the wells. In some cases however, when the number of wells is small or when the wells are located at a large distances from each other, the constraints in (2) may not effectively protect the aquifer. Mantoglou et al. (2004) proposed making the constraints rather stricter by modifying the second inequality in (2) as follows: xssi(z) 6 xwi, i = 1, . . . , k, where xssi corresponds to a somewhat larger value of potential than the one corresponding to the toe (see numerical examples in later sections). The optimization problem is then expressed as
55
ð3Þ
i ¼ 1; . . . ; k
Q i 6 Q max;i ;
The nonlinear constraint optimization problem (3) is solved using genetic algorithms which are stochastic search methods that mimic the metaphor of natural biological evolution (Holland, 1975; Goldberg, 1989; Fogel, 1994). They operate on a population of potential solutions applying the principle of survival of the fittest to produce improved approximations to a solution. A new set of approximations is created at each generation by selecting individuals according to their fitness and breeding them together using operators borrowed from natural genetics. This process leads to the evolution of populations of individuals that are better suited to their environment than their parents, just as in natural adaptation. Genetic algorithms model natural processes by applying the procedures of selection, recombination, mutation and migration at each iteration of the program. Since genetic algorithms solve unconstraint optimization problems, maximization problem (2) is modified in the following unconstraint form: Q tot ¼ ðQ 1 þ Q 2 þ þ Q k Þ þ pðzÞ
minimize : z
ð4Þ
where p(z) is a penalty term inserted in the optimization problem in order to handle the constraints. Penalty term p(z) is zero when the constraints are satisfied and gets a very high value when constraints are violated. In the following examples penalty term in (4) is defined as follows: pðzÞ ¼ 0;
when /i P 0 10
pðzÞ ¼ 10 ;
and x ssi 6 x w i
when /i < 0
or
xssi > x w i
ð5Þ
Fig. 2 presents the flow chart of the genetic algorithm optimization procedure. At the beginning of the computation, a number of individuals (the population) are randomly initialized and the first initial generation is produced. The objective function is then evaluated for these individuals. If the optimization criteria are not met a new generation is produced using the GA procedures of fitness assignment and parent selection, recombination, mutation and reinsertion. This cycle is performed until the termination criteria are reached. Such criteria are, for example, the number of generations or the maximum computer time. From the above discussion, it can be seen that genetic algorithms differ substantially from more traditional search and optimization methods. The most significant difference is that GA. search a population of possible solutions in parallel rather than a single solution. Additionally, GA do not require calculation of derivatives of the objective function. They require only the objective function and corresponding fitness levels which influence the direction of search. Therefore genetic algorithms can solve general classes of optimization problems even when the objective function is not continuous and the derivatives of the objective function do not exist. They are generally more straightforward to apply, because there are no restrictions for the definition of the objective function. They use probabilistic transition
56
A. Mantoglou, M. Papantoniou the second stage uses nonlinear programming methods (Sequential Quadratic Programming) as follows. The first optimization stage concerning optimization of well locations is expressed as
Initial Population of decision variables z Specification of the range of decision variables
maximize :
Evaluation of the objective function
x
k X
Q i ðxÞ
ð6Þ
i¼1
where the decision variables in (6) include the well coordinates x = [xw1, xw2, . . . ,xwk,yw1, yw2, . . . ywk]. Terms Q i ðxÞ denote the optimal pumping rates from wells i = 1, . . . k for given well locations x. Terms Q i ðxÞ are calculated from a second stage optimization with respect to pumping rates Q = [Q1, Q2, . . ., Qk] for fixed values of x. For known location of the wells x, the following objective function is maximized with respect to decision variables Q = [Q1, Q2, . . ., Qk] subject to constraints that limit saltwater intrusion into the aquifer
Selection Recombination Mutation Reinsertion
New Generation of z
maximize :
Q tot ¼ Q 1 þ Q 2 þ þ Q k
subject to :
/i ðQÞ P 0 i ¼ 1; . . . ; k x ssi ðQÞ 6 x w i ; i ¼ 1; . . . ; k
Q
Evaluation of the objective function
Are termination criteria met?
Q tot ¼
NO
YES Optimal values of well coordinates and pumping rates
Figure 2 Optimization flow chart for well locations using one stage optimization with genetic algorithms.
rules and can produce a number of potential solutions to a given problem. Thus, in cases where the particular problem does not have one individual optimal solution, the genetic algorithms are potentially useful for identifying these alternative solutions simultaneously. Optimization based on genetic algorithms is very flexible and, if the population size and the number of generations are sufficient, it converges to the global rather than a local optimum. However, it is often slow converging and it requires a large number of iterations. Convergence is improved by choosing reasonable initial population values. Increasing the number of individuals in the population and the number of generations improves the results but it requires more computer time. The genetic algorithm optimization program needed for calculation of the optimal well locations and pumping rates, was developed in MATLAB. MODFLOW was used for solving differential equation (1) and evaluating the required potential at each iteration of optimization.
Two stage optimization of well locations and pumping rates The second method of optimization follows two stages where the first stage concerns optimization of the location of wells and the second stage concerns optimization of pumping rates. The first stage uses genetic algorithms and
ð7Þ
Terms Q ðxÞ ¼ ½Q 1 ; Q 2 ; Q k in (6) represent the optimal pumping rates calculated from (7) for given values of x. Maximization of (6) with respect to x, is performed using genetic algorithms while the optimal values Q ðxÞ ¼ ½Q 1 ; Q 2 ; . . . Q k for each value of x are determined from maximization of (7) using nonlinear optimization methods. Since the decision variables in (7) are real and the objective function is continuous and differentiable with respect to the decision variables, nonlinear optimization method based on Sequential Quadratic Programming (SQP) is used for this second stage minimization (see Mantoglou et al. (2004)). Fig. 3 outlines the flow diagram of the two stage maximization based on genetic algorithms and Sequential Quadratic Programming. The range of the well coordinates and pumping rates are specified and an initial population [xini,1, xini,2, . . . xini,N] of N sets of well coordinates is generated. Given the well coordinate sets xini,m, the optimal pumping rates corresponding to xini,m are computed from (7) using SQP. Then a new population of well coordinates is generated based on the fittest members of the population using the selection, recombination, mutation, reinsertion and migration procedures of the genetic algorithm. When the termination criteria specified by the user are met, the GA optimization procedure ends. Such termination criteria are, for example, the number of generations or the computer time. The optimal solution corresponds to a specific set of coordinates which provide the maximum total pumping rate. As stated before, Sequential Quadratic Programming (SQP) method is used for calculation of the optimal pumping rates for specified well locations x. SQP is based on an iterative procedure where at each major iteration a positive definite quasi-Newton approximation of the Hessian of the Lagrangian function is calculated using the BFGS method (Goldfarb, 1970; Shanno, 1970; Fletcher, 1987; Luenberger, 1989). The optimization can be sensitive to numerical approximations such as truncation and roundoff error in the calculation of finite-difference gradients.
Optimal design of pumping networks in coastal aquifers using sharp interface models
57
Initial population of well coordinates x ini Specification of range of well coordinates
x
Initial values of well pumping rates
New Generation of coordinates
Specification of pumping rates
S.Q.P.
Reinsertion
Are termination criteria of S.Q.P. optimization met?
Mutation
NO
Recombination
YES
Selection
Optimal pumping Q* ( x ) for the specific coordinates Q* ( x )
Are termination criteria of GA. optimization met?
NO
YES Optimal well coordinates and pumping rates
Figure 3
Flow diagram for two-stage optimization using genetic algorithms and Sequential Quadratic Programming.
Furthermore, the program may converge to a local minimum rather than the global minimum. Starting the optimization with different initial conditions may help to locate the global minimum. Additionally, by changing the values of the parameters that take part in finite differences calculation or increasing the number of iterations contributes to the computation of the global minimum. Obviously, the above considerations have an important effect on the effectiveness of the two stage method that combines GA and SQP. The genetic algorithm and Sequential Quadratic Programming optimization programs needed for calculation of the optimum well coordinates and pumping rates were developed in MATLAB. MODFLOW was used to solve differential equation (1) and to evaluate the flow potential required at each iteration of optimization.
Application to an aquifer of orthogonal shape An aquifer having a simple rectangular geometry (Fig. 4) is selected in order to test and compare the two optimization methods. The recharge rate at the higher elevations of the aquifer is N = 150 mm/year distributed over an area of A = 9 km2 while the recharge rate at the lower elevations of the aquifer is N = 30 mm/day. The sea boundary constitutes a constant head boundary whereas the north and south boundaries are assumed impermeable. The aquifer thick-
ness from the base to the sea surface is taken as d = 25 m. The potential at the toe of seawater lens is calculated according to equation /s ¼ dð1þdÞ d 2 ¼ 8:0078 m2 . The geom2 etry of this aquifer resembles the aquifer at Vathi in the Greek island of Kalymnos described in detail in Mantoglou et al. (2004). The toe potential in this case is equal to /s ¼ 8:0078 m2 . When the number of wells distributed over the aquifer is relatively large, the constraints in (2) can effectively protect the wells from seawater intrusion by forming a barrier where the saltwater front cannot advance deeply into the aquifer and contaminate the wells. In some cases however, when the number of wells is small or when the wells lie at a large distance from each other, the constraints in (2) might not effectively protect the aquifer. Mantoglou et al. (2004) proposed to modify the second Mantoglou inequality in (2) and make the constraints somewhat stricter as follows: x ssi ðzÞ 6 xw i ; i ¼ 1; . . . ; k where x ssi corresponds to a somewhat larger value of potential than that corresponding to the toe. In the examples presented below, that setting /ss ¼ 8:1 m2 rather than /s ¼ 8:0078 m2 , make constraints x ssi ðzÞ 6 xw i ; i ¼ 1; . . . ; k stricter and effectively protect the wells. As is shown in the figures bellow the suggested sets of constraints allow the wells to yield the maximum pumping rates while protecting them sufficiently from seawater intrusion. More conservative solutions can be obtained by increasing /ss even more.
58
A. Mantoglou, M. Papantoniou
Figure 4
Table 1
GA GA & SQP
A test aquifer of rectangular shape.
Optimization results for three pumping wells for two optimization methods x1 (m)
y1 (m)
Q1 (m3/day)
x2 (m)
y2 (m)
Q2 (m3/day)
x3 (m)
y3 (m)
Q3 (m3/day)
Qtot (m3/day)
6350 5950
1950 2050
1378.3 1513.9
5550 6350
1250 1250
1327.6 1442.1
6350 6050
1450 1150
1604.2 1363.6
4310.1 4319.6
In a first application, the orthogonal aquifer of Fig. 4 with uniform hydraulic conductivity and recharge rate is selected in order to calculate the optimal coordinates of three wells and the optimal pumping rates. The results of optimization obtained from both methods of optimization are presented in Table 1. The first line of Table 1 presents the well coordinates and the optimal pumping rates computed by the first method based on GA while the second line presents the corresponding results obtained from the two stage optimization using GA and SQP. Both optimization methods yield a similar total pumping rate but the individual pumping rates and well
4
3000
16 12
8
2000 4 8 8
12
4
0 0
1000
2000
3000
4000
5000
16 20
8
8
1000
12 16
8
6000
7000
Figure 5 Equipotentials for the optimal well locations and pumping rates according to one stage optimization method.
12
3000 4
16
8 8
2000 12 8
4
8
16
8
0
1000
16
0
8
4
12 8
1000
2000
3000
4000
5000
6000
7000
Figure 6 Equipotentials for the optimal well locations and pumping rates for two stage methodology.
locations are different. Fig. 5 plots the equipotentials corresponding to the calculated optimal values with GA method while Fig. 6 plots the equipotentials corresponding to the optimal well coordinates and pumping rates evaluated using the two stage optimization methodology. While the two stage methodology yields a somewhat better Qtot than the one stage methodology, the computer time required is significantly larger.
Application to coastal aquifer in a Greek island The simulation and optimization methodology developed above is applied to an unconfined aquifer located at Vathi valley in the Greek island of Kalymnos (Fig. 7). As discussed in Mantoglou et al. (2004), the aquifer consists of high permeability limestone which outcrops at the highland areas while the impermeable bottom of the aquifer is composed by schist. The geologic formations found in the valley are highly permeable scree, medium permeability alluvium deposits and almost impermeable tuff (volcanic formations). Although in parts of the aquifer flow is taking place in the limestone fissures, we suppose for the purposes of this analysis that the aquifer can be represented by an equivalent porous medium. Based on existing wells the thickness of the aquifer is taken as: d = 25 m. According to Mantoglou et al. (2004) there are four zones with different surface recharge rates. The large recharge rate N = 150 mm/year in the limestone parts of the aquifer is justified by the large crevices present in the limestone facilitating percolation of surface water. The recharge rate at the scree area is smaller and is estimated at N = 70 mm/ year while the presence of clay in the alluvium deposits reduces the recharge rate to N = 30 mm/year. The tuffs are regarded as impermeable formations with zero recharge. Based on the geology of the valley, the aquifer is divided into four zones having different hydraulic conductivities (Fig. 8). The zonal hydraulic conductivities were evaluated by model calibration so that the model (1) predicts the measured fresh water levels in the existing wells (see Mantoglou et al. (2004)).
Optimal design of pumping networks in coastal aquifers using sharp interface models
Figure 7
Figure 8
59
Kalymnos aquifer.
The four aquifer zones with different hydraulic conductivity.
north boundary. Based on the above, the west, south and north boundaries are considered impermeable in the model. In the optimization case study, the aquifer is pumped by 11 wells and is replenished by surface recharge. Besides the pumping rates, the well locations are calculated using the optimization methodology. Table 2 presents the optimal pumping rates obtained by the one stage methodology while Figs. 9 and 10 plot the cor-
The sea along the east of the aquifer creates a constant hydraulic head boundary. The west and south boundaries of the aquifer correspond to groundwater divides, estimated based on hydrogeologic maps of the region, and are considered as impermeable boundaries in the model. Although the sea is very close to the aquifer in the north boundary, the existing wells are not contaminated by seawater which justifies the assumption of an impermeable boundary at the
Table 2 Optimal pumping rates (m3/day) for the optimal pumping network with 11 wells obtained from the one stage simultaneous optimization method with GA Well
1
2
3
4
5
6
7
8
9
10
11
Sum
Pumping rate
404.1
427.2
486.9
539.6
774.6
617.6
202.8
354.3
213.6
608.5
1101.2
5730.4
60
A. Mantoglou, M. Papantoniou
44
52
2 4
60 52
7 8 32
2000
40 36 2 3 5
28 24
3000
60
96 92
84 76 68
1
84 76 68
4000
92
5000
28
The above results indicate that both methods produce similar results regarding the optimal well location and optimal pumping rates (QGA = 5730.4 m3/day, QGA,SQP = 5639.7 m3/day, respectively). The two stage optimization method (using SQP and GA) requires more computer time than the one stage optimization. The tolerance level of SQP parameters in the two stage optimization was set to relatively coarse values so that the required iterations and the computation time are not excessively large. Decreasing the tolerance level and increasing the number of iterations of SQP, while increasing the population size and the number of generations of GA would im-
responding equipotential and the piezometric head maps. The optimal location of the wells is also shown in these figures. Table 3 plots the optimal pumping rates for the optimal pumping network with 11 wells obtained using the two stage optimization methodology which combines GA and SQP. The population size of the GA is equal to 30 and the number of generations is 50. Figs. 11 and 12 plot the optimal location of wells and the corresponding equipotential and the piezometric head maps, respectively. Fig. 13 plots the total maximum pumping rate obtained as a function of generation number evaluated from the genetic algorithm.
24 20 6 61
3 16
12
9
84
10 8
11
4
8
1000
1000
3000
4000
5000
6000
8000
9000
2. 62 5
1. 62 5 1. 37 5 1. 12 5 0.75 10
11
2000
3000
4000
5000
6000
0.5
1000
1000
0. 875
0.625
9
4
7000
0. 375 0.2 5 0.12 5
2000
6
1.375 1.1 25 0. 87 5
1.875
7
1.625
3.2 5 3
25 5 2. 87 2 1. 3 5 8
1 3.25 3 2.6 25 2. 25
3000
3.625
4000
3.6 25
5000
Figure 10
7000
Equipotentials for the optimal pumping network with 11 wells obtained from one stage method using GA.
3.8 75
Figure 9
2000
8000
9000
Piezometric head map for the optimal pumping network with 11 wells obtained from one stage method using GA.
Table 3 Optimal pumping rates (m3/day) for the optimal pumping network with 11 wells obtained from the two stage method that combines GA and SQP Well
1
2
3
4
5
6
7
8
9
10
11
SUM
Pumping rate
546.6
533.1
547.2
518.5
547.5
547.2
549.4
171.1
552.4
565.2
561.6
5639.7
Optimal design of pumping networks in coastal aquifers using sharp interface models
61
68 60 52
44 0 1 4 36 32
2
2000
24
3
5
28
24 20 16
4 6
28 24 2 1 0 12 6
44 36
3000
68 60 52
76
4000
80
80 76
5000
48 10
9
12
12
8
7
11 4
8
1000
1000
2000
3000
4000
5000
6000
7000
8000
9000
9
10
3000
4000
5000
6000
7000
0.5
2000
11 0.625
1000
1000
5 62 0.
8000
0.25 0.12 5
7
8
1.25 1.12 5 1 0.875 0.75
0.375
4 6
0.6 25
2000
75 1.3
1
3000
5 .87 1 1 25 1.6 2 3 5
0. 75
3
4000
1.25
5000
3 2.625 2 2.25 .62 5 1.87 5 2. 1.625 25
Equipotentials for the optimal pumping network with 11 wells obtained from the two stage method that combines GA
3. 3 3.2 75 5
Figure 11 and SQP.
4
9000
Figure 12 Piezometric head map for the optimal pumping network with 11 wells obtained from the two stage method that combines GA and SQP.
prove the solution obtained by the two stage methodology and it could outperform the one stage method at the cost of more computer time. Mantoglou et al. (2004) assessed the optimum pumping rates in the same aquifer from a set of 11 pumping wells with specified fixed coordinates from a non-optimal pumping network and found a maximum total pumping rate of about Qtot = 5360 m3/day. Notice that the total pumping rate obtained with fixed well locations from a non-optimal pumping network, is quite smaller than the pumping rate obtained here from the optimal pumping network where the well locations are included in the decision variables.
Summary and conclusions A method for optimal design of pumping networks in coastal aquifers was developed and utilized in a coastal aquifer in the Greek island of Kalymnos. The objective is to optimize the well locations and to maximize the total pumping rates
subject to constraints that protect the aquifer from saltwater intrusion. The method is based on genetic algorithms and nonlinear optimization subject to constraints that limit the saltwater intrusion into the aquifer. The simulation model is based on the sharp interface and the Ghyben–Herzberg approximations. The single potential formulation of Stack (Strack (1976) was used and the governing equations were solved numerically using finite differences (Mantoglou et al., 2004). The numerical model is more general than previous analytical solutions (Cheng et al., 2000; Mantoglou, 2003; Park and Aral, 2004) and it can handle aquifers of complex shapes, non-uniform hydraulic conductivity, non-uniform distribution of surface recharge, etc. The governing differential equations are solved using MODFLOW finite difference code. Two optimization methods are developed. The first method is based on simultaneous optimization of well locations and pumping rates using genetic algorithms
62
A. Mantoglou, M. Papantoniou 5650
5600
3
Qtotal (m /day)
5550
5500
5450
5400
5350 0
10
20
30
40
50
generation number Figure 13 to 30.
Maximum total pumping rate as a function of generation number for 50 generations. The population size of GA is equal
(GA). The second method combines genetic algorithms and Sequential Quadratic Programming in a two stage optimization. The well position optimization in the last method is carried out using GA while at each generation of GA a simpler optimization with respect to the pumping rates is performed using nonlinear optimization (Sequential Quadratic Programming). The evaluation of the constraints at each iteration is based on the numerical simulation. The simulation and optimization methodology was applied to an unconfined aquifer located at the valley of Vathi in the Greek island of Kalymnos for optimal design of a pumping network and for assessing the maximum pumping rates while protecting the wells from saltwater intrusion. The aquifer is modeled as heterogeneous with four zones of different hydraulic conductivity and four areas with different recharge rates. Using the two optimization methods developed in the paper, the optimal well locations and the maximum pumping rates were assessed. The methodology for optimal design of pumping networks based on a two stage optimization (first stage concerns well location optimization and second stage concerns pumping optimization) yields similar results with the one stage simultaneous optimization of well location and pumping rates but it requires more computer time. Decreasing the tolerance level and increasing the number of iterations of SQP while increasing the population size and the number of generations of GA would improve the solution of the two stage methodology and is expected to outperform the one stage method at the cost of more computer time. The pumping rates obtained from both methodologies are quite higher than the optimal pumping rate obtained from the pumping optimization method of Mantoglou
et al. (2004) based on pumping optimization with fixed non-optimal well locations. The proposed methodology, based on the sharp interface approximation, is applicable to regional scale problems where the transition zone is narrow relative to the scale of the problem. The proposed numerical solution of the governing differential equation does not have the limitations required by analytical solutions regarding orthogonal shape of the aquifer and homogeneity of hydraulic conductivity (Cheng et al., 2000; Mantoglou, 2003; Park and Aral, 2004) and can solve the problem of optimal well design in cases where the shape of the aquifer is of irregular geometry or when the hydraulic conductivity of the aquifer and the recharge rate is nonuniform.
References Ahlfeld, D.P., Heidari, M., 1994. Applications of optimal hydraulic control to ground-water systems. Journal of Water Resources Planning and Management 120 (3), 350–365. Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York. Cheng, A.H.-D., Ouazar, D., 1999. Analytical solutions. In: Bear, J., Cheng, A.H.-D., Sorek, S., Ouazar, D., Herrera, I. (Eds.), Seawater Intrusion in Coastal Aquifers – Concepts, Methods and Practices. Kluwer Academic Publishers. Cheng, A.H.-D., Halhal, D., Naji, A., Ouazar, D., 2000. Pumping optimization in saltwater-intruded coastal aquifers. Water Resources Research 36 (8), 2155–2165. Chiang, W-H., Kinzelbach, W., 2001. 3D-Groundwater Modeling with PMWIN: A Simulation System for Modeling Groundwater Flow and Pollution. Springer-Verlag, Berlin, Heidelberg, New York. Das, A., Datta, B., 1999a. Development of multiobjective management models for coastal aquifers. Journal of Water Resources Planning and Management 125 (2), 76–87.
Optimal design of pumping networks in coastal aquifers using sharp interface models Das, A., Datta, B., 1999b. Development of management models for sustainable use of coastal aquifers. Journal Irrigation and Drainage Engineering 125 (3), 112–121. Emch, P.G., Yeh, W.W.-G., 1998. Management model for conjunctive use of coastal surface water and ground water. Journal Water Resources Planning and Management 124 (3), 129–139. Essaid, H.I., 1999. USGS SHARP model. In: Bear, J., Cheng, A.H.-D., Sorek, S., Ouazar, D., Herrera, I. (Eds.), Seawater Intrusion in Coastal Aquifers – Concepts, Methods and Practices. Kluwer Academic Publishers. Fletcher, R., 1987. Practical Methods of Optimization. John Wiley & Sons. Fogel, D.B., 1994. An introduction to simulated evolutionary optimization. IEEE Transactions on Neural Networks 5 (1), 3–14. Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization and Machine learning. Addison-Wesley, Reading, MA. Goldfarb, D., 1970. A family of variable – metric methods derived by variational means. Mathematics of Computing 24, 23–26. Gorelick, S.M., 1983. A review of distributed parameter groundwater management modeling methods. Water Resources Research 19 (2), 305–319. Gorelick, S.M., Voss, C.I., Gill, P.E., Murray, W., Saunders, M.A., Wright, M.H., 1984. Aquifer reclamation design: the use of contaminant transport simulation combined with nonlinear programming. Water Resources Research 20 (4), 415–427. Hallaji, K., Yazicigil, H., 1996. Optimal management of a coastal aquifer in Southern Turkey. Journal of Water Resources Planning and Management 122 (4), 233–244. Holland, J.H., 1975. Adaptation in Natural and Artificial Systems. The University of Michigan Press.
63
Luenberger, D.G., 1989. Linear and Nonlinear Programming. Addison-Wesley Publishing Company. Mantoglou, A., 2003. Pumping management of coastal aquifers using analytical models of saltwater intrusion. Water Recourses Research 39 (12), 1–12. Mantoglou, A., Papantoniou, M., Giannoulopoulos, P., 2004. Management of coastal aquifers based on nonlinear optimization and evolutionary algorithms. Journal of Hydrology 297, 209–228. Mayer, A.S., Kelley, C.T., Miller, C.T., 2002. Optimal design for problems involving flow and transport phenomena in saturated subsurface systems. Advances in Water Resources 25, 1233– 1256. McDonald, M.G., Harbaugh, A.W., 1988. A modular three-dimensional finite-difference ground-water flow model. Technical Report, U.S. Geol. Survey, Reston, VA. Park, C.-H., Aral, M.M., 2004. Multi-objective optimization of pumping rates and well placement in coastal aquifers. Journal of Hydrology 290, 80–99. Shamir, U., Bear, J., Gamliel, A., 1984. Optimal annual operation of a coastal aquifer. Water Resources Research 20 (4), 435–444. Shanno, D.F., 1970. Conditioning of quasi-Newton methods for function minimization. Mathematics of Computation 24, 647–656. Strack, O.D.L., 1976. A single-potential solution for regional interface problems in coastal aquifers. Water Resources Research 12 (6), 1165–1174. Wang, W., Ahlfeld, D.P., 1994. Optimal groundwater remediation with well location as a decision variable: model development. Water Resources Research 30 (5), 1605–1618.