Aquifer Reclamation Design : Combined Simulation-Optimization Approach

A simulation-management methodology is demonstrated for the rehabilitation of aquifers that have been subjected to seawater intrusion. Finite difference groundwater flow and solute transport simulation are combined with nonlinear optimization. The model is capable of determining well locations plus pumping and injection rates of groundwater quality control. Restrictions was placed on hydraulic heads and concentrations. These restrictions was distributed over space and time. Two design strategies such as optimal pumping and recharge rates are demonstrated for a polluted aquifer. The method is not limited to these cases. It is generally applicable to the optimization of many types of distributed parameter systems.


Introduction
Although groundwater quality management models have been developed for important steady state and transient cases, research is needed for solution of nonlinear groundwater quality control problems such as capturing a contaminant plume that originate from waste disposal sites, saltwater intrusion control problem.Nonlinearity arise from management decisions that create unknown groundwater velocity fields as well as from problems involving chemical interaction.Nonlinearity constraints appear as a result of products of unknown velocity components which occur in advective and dispersive transport terms.Simulation of this nonlinear system has been successful.Research is needed to develop distributed parameter management models of saltwater intrusion that involve simulation of this nonlinear system.One such new approach is groundwater management methodology which optimizes aquifer remediation with the aid of nonlinear programming.In this paper a modeling approach is presented to determine the optimal design of reclamation schemes for saltwater intruded groundwater systems.The planning model combines a nonlinear, distributed parameter, groundwater flow and contaminant transport simulation model with a nonlinear optimization method.The resulting management model enables one to consider a linear or nonlinear water management objective subject to linear or nonlinear restrictions that include a complex, nonlinear aquifer simulation model.This paper deals with two scenarios which are demonstrated to control the seawater intrusion.

Combined Simulation-Optimization Models
Groundwater hydrologists developed simulation models to aid in predicting the migration pattern and fate of contaminants in aquifers.The application of finite difference and finite element methods to groundwater equations has permitted complex, real world systems to be modeled.Numerical simulation models have enabled hydrogeologists to understand the functioning of regional aquifers and to test hypothesis regarding the behaviour of particular facets of groundwater systems.The simulation method has provided a framework for conceptualizing and evaluating aquifer system.Models have become tools to evaluate the long term impacts of sustained water withdrawals, groundwater-surfacewater interaction and the migration of chemical contaminants.Clearly, simulation is a method to explore hydrogeologic problems and a tool to predict impacts upon groundwater systems will continue to be essential to hydrologists and to water managers.However, simulation models are often utilized to explore groundwater management alternatives.In such cases a model is executed repeatedly under various design scenarios which attempt to achieve a particular objective, such as isolating plume of contaminated groundwater, preventing saltwater intrusion, dewatering an excavation area or obtaining a sustainable water supply.Use of such an approach often sidesteps rigorous formulation of groundwater management goals and fails to consider important physical operational restrictions.Determining the proper objective function in a groundwater management model is often difficult but an essential aspect, which should not be avoided.It is unlikely that optimal management alternatives are discovered using simulation techniques alone.The requirement is not a simulation model alone, but a combined simulation and management model.A combined simulation and management model considers the particular behaviour of a given groundwater system and determines the best operating policy under the objectives and restrictions.Remson and Gorelick (1980) considered the optimal locations of minimum pumping rates for the hydraulic containment of contaminated groundwater.His approach was limited utility for problems involving aquifer rehabilitation.The reason is that the model only considers management of the groundwater flow field and does not explicitly contaminant concentrations.The management of pollutant sources under steady state conditions has been addressed by Willis (1976) and Gorelick and Remson (1982).In each of these studies the equations governing groundwater solute transport defined the constraints of a linear program in a management model aimed at maximizing safe waste disposal potential.The dynamic management of pollutant sources was studied by Gorelick et al. (1979), Willis (1979) and Gorelick (1982).Various approaches were adopted in these studies, but in general each relied upon using a simulation model to define constraints in a linear management model.The most severe limitation of the above groundwater quality management models is that the groundwater flow field must be determined before pollutant sources can be managed in order to maintain system linearity.The result is that restrictions could not be imposed simultaneously upon groundwater hydraulics and groundwater contaminant transport.For instance one cannot use the above models to determine a pumping and recharge scheme that would ensure that contaminant concentrations were maintained within purview of specified water quality standards.

Previous Works
Nonlinear programming techniques were used to incorporate fully into groundwater management process (Gorelick et al., 1984, Ahfeld, 1987).The nonlinear techniques were necessary because solute transport modeling is dependent on the velocity of groundwater, which inturn is affected by optimal pumping rates during optimization.Concerns about computational burden in field scale applications, global optimality and difficulty of handling multiple objectives of different cleanup technologies have led researchers to continue improved methods of optimizing groundwater remediation.The original applications (Gorelick et al., (1984) ;Ahfeld, 1987) and most of the studies since then have used MINOS (Murtagh and Saunders, 1980) or NPSOL (Gill et al., 1986).Both are well known optimization algorithms that use the gradient method combined with a projected Lagrangian to manage nonlinear constraints.Ahfeld (1987) found that Central Processing Unit time increases nonlinearly with the number of decision variables and logarithmically as the desired cleanliness increases.MINOS and NPSOL have difficulties to guarantee a global optimum where the problem is nonconvex, which is the case for most field applications.
In light of these concerns, two approaches that have found recent success in other areas of combinatorial optimization, Simulated annealing and Artificial Neural Network, are of interest to examine for use in groundwater research.Simulated annealing was developed by analogy to the statistical mechanics of annealing solids; Artificial Neural Network was developed by analogy to the collective processing behaviour of neurons in the human brain.Both the techniques believed to be robust with techniques to avoid being trapped in the local minimum of a function.Dougherty and Marryott (1991), Marryott et al., (1993) and Marryott (1996) applied the simulated annealing methodology to groundwater management problems and expressed that the relative Central Processing Unit requirements for the simulated annealing approach increased at a slower rate than the gradient-based method of MINOS.They discussed that simulated annealing provides an optimal framework that is flexible enough to incorporate a number of different remedial technologies into the design process.Roger and Dowla (1994) used artificial network to design optimal remediation systems.Ritzel et al. (1994) solved a multi-objective groundwater contamination problem using genetic algorithms.In their study only steady state flow was simulated.Mckinney et al., (1994) applied genetic algorithms to groundwater resource management and design of pump and treat systems in a hypothetical aquifer with a single management period.Johnson and Rogers (1995) combined neural networks and genetic algorithms to select the optimal pumping well locations of groundwater remediation problems.Barlow et al., (1996) used nonlinear programming solver to maximize pumping and found that it is highly efficient when compared with simulation alone.Andricevic (1990) developed a penalty type cost function of two conflicting objectives, aimed to minimize deviation from the desired withdrawal rate and the specified target levels.Cheng et al., (1992) used hyperbolic penalty function to remediate the contaminated aquifer system.Mckinney and Lin (1995) used a polynomial penalty function to minimize the cost of pump and treat aquifer remediation system.Keshari and Datta (1996) used exterior penalty function method along with pattern search Hooke and Jeeves method to solve the unconstrained optimization problem.
In this research a groundwater remediation management model is developed using penalty function method which is used to convert the constrained optimization method to an unconstrained optimization problem.The pattern search algorithm proposed by Hooke and Jeeves is used to solve the resulting unconstrained optimization problem.For the development of this model, three dimensional groundwater flow model MODFLOW (McDonald and Harbaugh, 1988), and solute transport model MT3D (Zheng, 1990) are used as simulation models so that the coupled simulation-optimization model can be applied to real world problems.The methodology is general and the results and suggestions are in particular to Minjur-Mouthambedu Aquifer System, Chennai, India are much useful in decision making to control seawater intrusion and to utilize the resources effectively by the system managers.
Minjur-Mouthambedu Aquifer System consists of alluvial deposit underlined by tertiary rocks, which in turn overlies the massive Gondwana series [United Nations Development Programme Report, (1975)].The aquifer is bounded on the eastern side by the Bay of Bengal, the northern and southern sides by the impervious formation and in the western side by alluvial deposits as shown in Fig. 1.The permeability of the Chennai aquifer varies from 5*10 -4 to 1.21*10 -3 m/s from east to west.Central portion of the aquifer is highly porous having higher ranges of permeability.The aquifer system is conceptualized as a two layer system with an upper and lower formation, separated by an aquitard.The upper formation is an unconfined aquifer and the lower formation is a confined aquifer capped with clay beds and connected to the upper aquifer by vertical leakages.The recharge of the unconfined aquifer is from rainfall and irrigation return flow.The recharge for the confined aquifer is longitudinal and infiltration from exposed parts of lower aquifer.The withdrawal of water for irrigation, industrial and domestic uses are mostly from the confined aquifer.
The main reason for the intrusion was due to over extraction of groundwater to meet the various demands.Heavy pumping in the Minjur -Mouthambedu Aquifer System has created pronounced reversal of hydraulic gradient from sea to land, resulting rapid increase of chloride concentration in the groundwater.Due to continuous concentrated pumping, the interface further moved inland and its movement was back and forth depending on the discharge and recharge.The solute transport modeling involved two parts.The predicted head distribution computed from the flow model is utilized in the solute transport model to estimate the chloride concentration distribution.The model is calibrated for the period 1976 to 1982 and it is tested for the period 1983 to 1996.The head and concentration contours at the end of December 1996 was taken as initial condition for the beginning of the year 1997 and the system was projected for the future period 1997 to 2020.

Optimization Methodology
A key component of simulation-optimization modeling is the formulation of the optimization problem that requires definition of management objective, decision variables and management constraints.In this work, the goal is to reduce the salt concentration at the critical zone.This can be achieved through two formulations.One is to find the optimal reduction in the pumping and the other is by optimal recharge for which the concentration is reduced at the end of planning horizon.

Decision variables
Definition of the optimization problem includes decision variables whose values will be determined as part of the optimization solution.In a remediation system, there are two sets of variables.They are decision variables and state variables.Decision variables include the pumping and injection rates for the wells.Other possible decision variables include well locations.The purpose of the design process is to identify the best combination of these decision variables.The state variables are the hydraulic head and the chloride concentration.Any remediation design model shall include two major components.The simulation model updates the state variables and the optimization model selects the optimal values for the decision variables.

Managerial constraints
The formulation of optimization problem includes the definition of management constraints.Usually the constraints are placed on the piezometric heads, total recharge (or pumping) and chloride concentration in the case of quantity and quality management problems.
i.The temporal and spatial distribution of piezometric heads is stated such that it should not drop below the specified minimum values and the same should not rise above the specified maximum values; ii.Total recharge (or pumping) must be within the specified range; iii.The temporal and spatial distribution of pollutant concentrations should not exceed the specified threshold values so that the water quality standards in the aquifer will not fall below the specified value (usually zero) and the concentrations for the intended use; Here h lb , R lb , and c lb are the lower bounds on hydraulic head, artificial recharge and concentration, whereas h ub , R ub , and c ub are the upper bounds on hydraulic head, pumping and concentration, respectively and k denotes the time step.

Objective Function
In most of the pollutant management problems, the goal is to find either the minimum recharge quantity or minimum reduction in pumping that satisfies a set of constraints to maintain the head, recharge and concentration in a specified range.To modify the formulation of the unconstrained nonlinear problem, the constraints are brought into the objective function through penalties.The objective function can be written as where, dh is the deviation of the head from the specified range at locations (i,j); dc is the deviation of the concentration from the specified range at locations (i,j); dr is the deviation of the recharge from the specified range at locations (i,j); (i,j) is a set of specified locations; k is the number of time steps considered within a planning horizon; ph is the penalty for head deviation; pc is the penalty for concentration deviation; pr is the penalty for recharge deviation;

Method of solution
To find the optimal solution, the constrained problem is transformed into an unconstrained problem using penalty application.The unconstrained nonlinear problem is then solved using the Hooke and Jeeves direct search algorithm (Rao, 1991).Hooke and Jeeves method is a sequential technique where each step consists of two kinds of moves, one called the exploratory move and the other called the pattern move.The first kind of move is included to explore the local behaviour of the objective function and the second kind of move is included to take advantage of the pattern direction.This method requires only objective function evaluation and does not require the partial derivatives of the function in finding the optimum point.It starts from an initial trial solution and the step length along the decision parameter axes is kept constant for each cycle of moves and a probe is made first in the positive direction and then in the negative direction of each axis.
Iterative improvement can stuck in a local minimum.Consequently, the solution depends upon the starting configuration.Hence, several different starting points are used to make sure that a better solution is found.A simulation model is used as a sub-model with the Hooke and Jeeves nonlinear programming main model, which passes different configuration of the decision vector to the simulation model and receives an objective function value back.Based on the new value of the objective function, the Hooke and Jeeves model decides to move to the new configuration of the management decision vector.This is repeated until no further improvement in the objective function is possible.

Application
As the aquifer is intruded by seawater, it is essential to control the seawater intrusion either by reducing the pumping or by recharging artificially.The objective of the two strategies is to minimize total pumping and total injection.The transient reclamation design model is demonstrated using complex, heterogeneous aquifer system shown in Fig. 1.

Scenario 1: Optimal reduction in pumping
The aim of this formulation is to find the optimal reduction in the irrigation demand by maintaining the constraints on the head and concentration.This scenario is attempted to know the percentage of pumping in the year 1996 that can reclaim the aquifer optimally.Initially the constraints are kept on both heads and concentrations.After a number of trials, it was found that the reduction in the pumping is directly proportional to the cleanup criteria.The optimal pumping is the zero pumping for which maximum cleanup takes place.The head distribution is not uniform so constraints on head are introduced.
The objective function which constitutes the unconstrained minimization problem for the proposed management model can be expressed as minimize Σph * dh ∀ i,j,k (5)

Results of scenario 1
The optimal reduction in pumping is 28 percent of 1996 irrigation pumping for which the penalty value is 623.All the other neighbouring points show higher penalty values as shown in the Fig. 2. The optimal pumping was found out through bisection method.The constraints are piezometric heads, which should lie between 0.5 and 1.0 m.This range is selected on the idea that maintaining piezometric pressure above mean sea level reduces or arrests seawater intrusion.The deviations from these constraints are penalized.Therefore, the objective is to minimize the penalties.Penalties are calculated for various percentage reductions in pumping.Penalties for zero and 100 percent reduction in pumping are 2781 and 6894 respectively.The bisection trials were carried out in between zero percent and hundred percent reductions in pumping.Penalties for various reductions in pumping were arrived at.Fig. 2 shows the change in penalty with reduction in irrigation demand.
The optimal reduction in pumping gives the minimum deviation from the constraints on head, that is, 28 percent reduction in pumping achieves the objective.After finding the optimal percentage reduction in pumping, the same was used in the simulation model to arrive at various results, such as head and concentration contours, water and salt balance and temporal variations of head and concentration as shown below.
The minimum piezometric head is improved by +8 m in the first nine years and the same was continued.Similarly, the reduction in concentration was steep in first nine years from 8329 mg/l in 1996 and then it reduces slowly and it reaches to 3783 mg/l at 6.5 km in 2020.Fig. 3a and 3b and Table 1 show the improvement in head and concentration; The 1000 mg/l isochlor nearly at 10 km in 1996 retreats back to 8.5 km in 2020 from the coast.Hence, the 1000 mg/l isochlor front moves back by 1500 m as shown in Fig. 4; If the pumping is cut down by 28 percent per year, there is an outflow of 120 mcm of either fresh groundwater or diluted seawater which washes out 0.95 million kilogram of salt in 24 years through sea boundary.When the zero percent reduction in pumping is compared with 28 percent reduction in pumping, the condition was reverse.It provides 131 mcm of saltwater that brings 4.5 million kilogram of salt in 24 years into the aquifer.

Scenario 2: Optimal recharge and location
The objective is to minimize the sum of temporal and spatial artificial recharge fluxes for a period of ten years.The purpose of this model is to minimize the total amount of recharge to the aquifer in order to reclaim the seawater intruded aquifer.
i.The objective is to find the optimal (minimum) recharge rate to reclaim the polluted aquifer; ii.Though a 24 year time horizon is considered, only for last ten years, the values of head, concentration and recharge are used in the optimality check.Due to the imposed artificial recharge, the aquifer system may be in a disturbed condition in the initial years and hence the first fourteen years data on head, concentration and recharge are not considered in the optimality check; iii.One year is taken as a cycle.The recharge pattern between cycles is kept constant.Each annual cycle is divided into four stress periods based on rainfall pattern.The four distinct seasons are (i) January to May, (ii) June to September, (iii) October and November and (iv) December; Division of a year into four stress periods leads to four decision variables and the other decision variable is the total recharge.The first four decision variables denote the sharing of the total recharge while the fifth decision variable is the total recharge itself.Thus this becomes a six dimensional problem in which the sixth axis is the penalty axis.The recharge in each stress period is distributed among 20 cells of the discretized aquifer.The simulation model uses all the aquifer parameters and boundary conditions that are finalized after calibration and testing of the aquifer.In this scenario, the optimal recharge is found without adjusting the 1996 pumping pattern.
The objective function that constitutes the unconstrained minimization problem for the proposed management model can be expressed as

Penalty values
The piezometric head h i,j is constrained in the range of 0.5 to 1 m, in order to maintain the piezometric head above the mean sea level.The penalty for the deviation from this range is taken as follows.
Head h (m) Penalty (ph) h < 0.5 (0.5-h) * 10 The recharge quantity is constrained to be less than 20 mcm/year.This has been decided based on the availablity water at the recharge zone.The penalty for the deviation from the specified range of recharge is taken as follows.
Recharge R (mcm/Year) Penalty (pr) R ≤ 20 0 The temporal and spatial distribution of pollutant concentration in the aquifer should not fall below the specified value (usually zero) and these concentrations should not exceed the specified threshold values to meet the water quality standards for the intended use.The penalty for concentration is fixed as follows.
Concentration C (mg/l) Penalty (pc) C < 1000 0 C > 1000 (C -1000) * 1 The nonlinear optimization problem was solved with Hooke and Jeeves method combined with groundwater flow model and the solute transport model.To represent the behaviour of the groundwater flow and chloride concentration in the management model, flow and solute transport models were linked explicitly to the optimization program as a subroutine shown in Fig. 5. Flow and transport models are called repeatedly by Hooke and Jeeves model.Hooke and Jeeves model determines the recharge rates.These recharge rates are passed to flow and solute transport model and progresses towards an optimal point by comparing heads and concentrations with the previous point.This procedure will be continued until optimum point is determined that minimizes the objective.

Location analysis
After finding the optimal total recharge quantity and its distribution, optimal location analysis was carried out.In this aquifer system, two troughs in piezometric surface, one at 6.5 km and the other at 13.5 km from the coast are existing.It is obvious that through control of these troughs, it is possible to control seawater intrusion.The optimal location for recharge was searched between 3.0 and 13.0 km from the coast by using bisectional search method.

Results of scenario 2
The total recharge was initially distributed equally to all the stress periods.About 20 starting points are considered and for each starting point, a best point (local optimal point) is obtained.All the 20 starting points led to the same optimal point.The results are shown in the following Table 2. Optimal recharge quantity of 8.1 mcm/year and optimal location at 6.5 km from the seacoast were arrived through this analysis.
The optimal recharge quantity and its distribution at optimal location was fed into the simulation model to arrive at the various detailed results from which head and concentration contours, isochlor front movement and variations of heads and concentration were drawn.The summarized results of the scenario 2 are shown in the Table 3.
The minimum piezometric head improves steadily from -8.0 m to -4.65 m during the period 1996 to 2020 due to artificial recharge.The reduction in concentration is very steep in the first ten years from 8329 to 2264 mg/l and then it reduces steadily to 502 mg/l as shown in Table 3; The 1000 mg/l isochlor front occupies 9.84 km in 2005 with a slight retrieval of the front as it moves back to 6.94 km in 2020.There is a large retrieval of the front in these five years and then it moved back slowly.On the whole the isochlor moved back by 3390 m in 24 years as shown in Fig. 6; For the 8.1 mcm/year of artificial recharge it was able to wash out half million kilogram of salt through the outflow of diluted seawater quantity of 62 mcm; Fig. 3a indicates that the piezometric head fluctuates steadily.The rate of improvement in the piezometric head for the scenario 1 is better than the scenario 2, because the pumping is cut down throughout the aquifer system; Fig. 3b indicates the variations in the chloride concentration at various locations.In all locations the reduction concentration is more in scenario 1 compared to the scenario 2 except at location 11526.As the recharge is at 6.5 km near the well 11526, the clean up is more in scenario 2 compared to scenario 1;

Discussion
From the scenario 1 and 2, it is evident that the aquifer is not reclaimed completely.It is not possible to reclaim the aquifer within a short time horizon.It will take much longer time.Always the piezometric head reaches steady state quickly whereas the salt concentration will be always in transient condition because of dispersion and diffusion.When the scenarios 1 and 2 are compared, the reduction in pumping takes away more salt from the system than the artificial recharge.The 28 percent reduction in pumping reduces irrigation pumping by 10.4 mcm/year.This reduction in pumping takes away more salt than the 8.1 mcm/year recharge washout.Artificial recharge is expensive and hence reduction in pumping is suggested.This is possible by either changing the agricultural pattern or buying the groundwater rights from the farmers.Then the question of farmers employment arise.It is possible, that they can be employed under special category.Buying the groundwater rights may lead to unemployment problems, which must be solved with alternate arrangements to them.

Conclusion
A general methodology for the management of nonlinear distributed parameter systems has been demonstrated.The procedure enables the solution of a previously intractable class of nonlinear groundwater management problems.The particular problem addressed here regards to reclamation schemes for contaminated groundwater systems.This approach is to unite a contaminant transport simulation model with a nonlinear optimization method.The design model is capable of identifying well locations and determining pumping and recharge rates for optimal aquifer restoration design.Useful extensions will include more complex nonlinear simulations, such as those involving multi-component chemical interactions and heat transport..

Notations
Figure CaptionsFig.1 Study area and its location Fig. 2 Change in penalty with reduction in irrigation demand Fig. 3a Fluctuations of piezometric head for the scenarios 1 and 2 at the well 11526 Fig. 3b Variations of projected chloride ion concentration for scenarios 1 and 2 at the well 11526 Fig. 4 Projected isochlors at an interval of 1000 mg/l for scenario 1 Fig. 5 Flow chart for the combined simulation-optimization model Fig. 6 Projected isochlors at an interval of 1000 mg/l for scenario 2

Table 1 .
Summary of results for scenario 1

Table 2 .
Input parameters and penalty

Table 3 .
Summary of results for scenario 2