OPTIMUM MANAGEMENT OF GROUNDWATER PUMPING BY USING SIMULATED ANNEALING TECHNIQUE

A two dimensional mathematical model is developed to simulate a groundwater flow regime of the upper part of Quaternary Deposites, in the Teeb area, Missan province. The optimum management of groundwater flow is developed by a linking simulation-optimization model. In this study, MODFLOW, 98 packages are used to simulate the groundwater flow system in unconfined aquifer. This model is integrated with an optimization model (Modular Groundwater optimizer) which is based on the simulated annealying (SA). Fixed well location case was undertaken by running the model with adopted calibrated parameters. The optimum value of pumping rate estimated through this study is (108320 m 3 /day), this value is compared with another previous work which based on genetic algorithm GA for the same area. The results of comparison showed that the optimum value obtained by SA is higher than the one obtained from GA by 18.35% for the same piezometric head draw down that indicate, the SA gives higher value than GA in this study.


INTRODUCTION
At this time, many of the concerns about groundwater resources on or adjacent to public land involve questions about depletion of ground water storage. Increasingly, attention is being placed how to manage ground water resources on public lands in a sustainable manner.
The development of computer simulation models has paralleled the development and increasing availability of faster, more capable, yet less expensive computer systems. One of the most popular and comprehensive deterministic groundwater models available today is the United States Geological Survey's (USGS) MODFLOW modular (Mc Donald & Harbaugh 1988, Delleur, 2007. The model was coded in a modular style to allow and encourage the development of additional packages or modules that can be added on or linked to the original code. It is organized by "Processes", which currently include the groundwater flow, observation, sensitivity, parameter estimation, groundwater management, and groundwater transport processes. Groundwater models often are developed to address specific questions (or problems) that arise in the management of groundwater system. For example: What is the maximum amount of water that can be pumped from a particular aquifer without casing unacceptable groundwaterlevel declines or reductions in stream flow? What is the least-cost pumping strategy to capture or contain a contaminant plume? What is the best strategy to conjunctively use the groundwater and surface-water resources of a basin? (Delleur, 2007).
During the past two decades, much effort has been directed toward developing computer models that can simultaneously account for the physical processes that affect flow and pollutants transport in a groundwater system and the large number of engineering, legal, and economic factors that serve development and management of groundwater resources. Such models couple the capabilities of groundwater quantitative and qualitative simulation with mathematical optimization techniques (such as linear and nonlinear programming) that solve properly formulated management problems. An optimization problem (or formulation) is a model of the management decision-making process, and has three primary components-an objective function, a set of constraints, and a set of decision variables (Ahlfeld & Mulligan, 2000). The decision variables of a management problem are the quantities that are determined by solution of the management problem, for example, the rates, locations, and timing of groundwater withdrawals from a set of available pumping wells. As implied by its name, the objective function defines the objective of the management problem such as "maximize total pumping from the set of pumping wells". The constraints of the problem limit the values that can be taken by the decision variables in the solution of the problem.
The Modular Groundwater Optimizer (MGO) code couples groundwater flow and contaminant transport simulation models MODFLOW and MT3DMS with mathematical optimization techniques to address groundwater management problems (Zheng & Wang, 2003). The development and application of the couple simulation-optimization approach has been an active and fruitful research area in recent years. The simulation-optimization approach is appealing because it can both account for the complex behavior of a groundwater flow system and identify the best management strategy to achieve a given set of prescribed constraints (Wanger, 1995).
The purpose of this research is to link simulation-optimization model for obtaining the optimum management of groundwater flow in Teeb area.
The MGO is used in this research. In the proposed model, MODFLOW packages are used to simulate the flow of the groundwater system. This model is then integrated with an optimization model which is based on the simulated annealing (SA).
SA is a local algorithm (met heuristic) capable to locate global or near-global optima with a high degree of reliability. Many researchers used SA technique in groundwater optimization. Dougherty and Marryott (1991) provide an overview of the simulated annealing algorithm and guidance on its application to groundwater management. They note the simplicity with which simulated annealing is implemented and coupled to flow and transport simulators. Marryott et al. (1993) apply SA to an actual contaminated field site. SA coupled with flow and transport simulator is used to solve two problems. The first is a regulatory constraint formulation in which concentration levels in the aquifer are reduced below regulatory standards by the end of a designated planning period. The second is a velocity constraint formulation in which the goal is plume containment.

STUDY AREA
Study area is located in north and north east of Missan province as shown in Fig. 1. It occurs along the foot of mountains of the Iraqi-Iranian frontier in south of Iraq, between longitudinalline (47°06'-47°36') and latitude-line (32°06'-32°30'). The considered area is about 1860 km 2 . It extended from Teeb area close to the Iraqi-Iranian border to Shikh Fars area. All boundaries for this study area were modeled as head-dependent boundary to allow inflow to the modeled region at a rate proportional to the head difference between the aquifer outside the simulated area and the model boundary. The top of the model was represented as unconfined aquifer. The water table elevation changes as part of the model solution. The bottom of the model was represented as a no flow condition. The vertical location of this boundary was selected to correspond with the base of the aquifer (aquitrad).

OPTIMIZATION TECHNIQUES
Two sets of variables related in groundwater management problems, decision variables and state variables. Decision variables represent the injection rates or pumping of well and it may be include the "on/off" status of a well and well locations. The best combination of decision variables refer to the optimal management strategy or policy. On the other hand the state variables represent the hydraulic head, which is the dependent variable in the groundwater flow equation. In a connected simulation-optimization model, the simulation model updates the state variables, and the optimization model determines the optimal values for all the decision variables. The objectives of the management must be achieved within a set of constraints. The constraints may be state variables or decision variables and they may take the forms of either equalities or inequalities. The number of candidate wells and the upper and lower bounds for pumping/injection rates, and the candidate well location refers to the constraints on the decision variables. Constraints on the state variables involve the requirements that hydraulic head be maintained above or below a certain level.
A general form of the objective function and a set of commonly used constraints suitable for a wide variety of resources management design problems can be expressed as follows (Zheng & Wang, 2003).

Maximize (or minimize)
Where, in equation (1); J is the management objective in terms of the total costs or in terms of the total amount of pumping or mass removal.
Qi is the pumping/injection rate of a well-represented by parameter i (negative for pumping and positive for injection). Note that the term parameter (Qi) is used to represent the pumping/injection rate associated with a particular well location at a specific management period. For an optimization problem with only a single management period, the flow rate of any well is constant and can be represented by a single parameter. However, for an optimization problem with multiple management periods, the flow rate of any well can vary from one management period to another. Thus, multiple parameters are needed to represent the flow rates of the well at different management periods.
F(q,h,) is any user-supplied cost function which may be dependent on flow rate q and hydraulic head h.
yi is a binary variable equal to either 1 if parameter i is active (i.e., the associated flow rate is not zero) or zero if parameter i is inactive (i.e., the associated flow rate is zero).
di is the depth of well bore associated with parameter i.
Kufa Journal of Engineering, Vol. 7, No. 2, 2016 13 Δti is the duration of pumping or injection associated with parameter i (or the length of the management period for parameter i).
a1 is the fixed capital cost per well in terms of dollars or other currency units; a2 is the included the installation and drilling cost (dollars or other currency units) per unit depth of well bore (e.g., dollars/m); and, a3 is the pumping and/or treatment costs (dollars or other currency units) per unit volume of flow (e.g., dollars/m 3 ).
a4 is the multiplier for an external user-supplied cost function.
Among the constraints equations.
N is the total number of parameters (decision variables) to be optimized.
Equation (2) is a constrain stating that the total number of actual wells at any time period must not exceed a fixed number, NW, out of the total candidate wells.
(Qmin and Qmax) are the minimum and maximum flow rate of a well at any specific management period.
(hmin and hmax) are the lower and upper bounds of the hydraulic head.
(hm) is a hydraulic head at any monitoring location.
(ℎ − ℎ ) is the head difference between an "outside" and an "inside" of considered area, monitoring wells must be greater than a minimum value, Δhmin.
Equation (6) is a constrain stating that the pumping/injection rate of a well at an arbitrary location, Qm, is proportional to the sum of the optimized flow rates represented by parameters I1 through I2 where A and B are proportional constants.

SIMULATED ANNEALING
In Simulated Annealing (SA) (Zheng & Wang, 2003), one can simulate the behavior of a system of particles in thermal equilibrium using a stochastic relaxation technique developed by Metropolis et al. (1953). This stochastic relaxation technique helps the SA procedure escape from the local minima. In a conventional local search method, the search procedure always chooses the best solution among its neighbors. If the best solution is not better than the current solution, the search stops and assumes the optimal solution has been found. The method is efficient for objective functions that are simple and have only one such local minimum point (in a minimization problem). For complex functions, this local minimum point may be quite different from the global minimum, and in this case, the search is said to be trapped in the local minimum point. Through stochastic relaxation, SA has the ability to escape from local minimum points. The essential feature of SA lies in allowing the procedure to move both "uphill" and "downhill". The symbols used in SA technique are shown in Table 1.
The basic steps of SA are outlined in Table 2. There are an inside do-loop and an outside doloop. The outside do-loop reduces the temperature each time by a factor of λ; the inside loop runs the simulation NSimPerIter times under a constant temperature. In step 4, the new solution can be generated randomly or in a deterministic fashion. SA stops either when it has completed all the iterations or when NTOL is reached, i.e., the maximum number of iterations allowed.  The set of control parameters involved in simulated annealing is called an annealing schedule. Important control parameters are the initial temperature T0 , the maximum number of iterations under a constant temperature NSimPerIter, and the temperature reduction factor λ. Determination of the values of the control parameters requires extensive experiment. To guarantee the convergence of the annealing process, the annealing process should start at a high initial temperature T0 and then the temperature T should be dropped slowly. In addition, the number of iterations under each constant temperature NSimPerIter should be sufficiently large. But higher initial temperature and larger number of iterations mean longer simulation time and may severely limit computational efficiency. As a rule of thumb, to balance the computational efficiency and quality of the solution, the initial temperature should be chosen in such a way that the acceptance ratio is somewhere between 70% and 80%. It is recommended that a number of simulation runs should be executed first to estimate the acceptance ratio. Second, the number of iterations under a constant temperature is about ten times the size of the decision variables (Zheng & Wang, 2003).
A major advantage of simulated annealing is reliability in locating global or near-global optima.
Other noted advantages are its flexibility in considering complex, highly nonlinear or discontinuous objectives and/or constraints, and the ease with which it can be interfaced with application-specific simulators.

MANAGEMENT MODEL APPLICATION
The objective of this research is to identify target groundwater supply area having potential elements in the study region. The management case applied for fixed well location in the study area by running the model with adopted calibrated parameters.
There is a thirty five actual number of pumping wells in the study area, having locations as shown in Fig. 2. This problem can be formulated as an optimization problem with the following objective function and constraints.
Subject to And, Where: J :-the objective function (maximize pumped/injected water volume) :-the absolute pumping rates ∆ :-the length of stress period in the flow model (where the number of stress period is equal to (7). ∆ is equal to 30 day).
ℎ :-the initial hydraulic head The location and the hydraulic heads of monitoring wells are presented in Table 3. The number of monitoring wells is nine in the present model. Nine wells are taken in the present model for monitoring the hydraulic heads of groundwater across the study area as shown in Fig. 3. Depending on the observation period (one year) in the study area, the change in the groundwater level is 0.5 m (see Eqs. 10 ,11 ), and it is more realistic to take a value corresponding or greater than this hydraulic head (Al-Aboodi, et al., 2013).
The magnitude of each pumping rates to be optimized ranges from zero as the minimum and 4000 m 3 /day as the maximum, as it shown in equation 9. Generally, the selection of an appropriate value for use as the maximum pumping rate needs several test runs. If it is set too high, the optimization solution may be inefficient.
The number of discretization intervals (NPStep) for each pumping rate parameter is chosen to be 26. Since the minimum and maximum values for each parameter have been specified as 0 and 4000 m 3 /day, respectively, the precision (or resolution) of the identified pumping rates is (4000-0)/ (26-1), or 160 m 3 /day. In other words, the final pumping rates obtained by the SA solution may differ from the actual optimal values by as much as (but not to exceed) 160 m 3 /day. The number of simulation per optimization iteration (NSimPerIter), in SA is set at 100. The initial temperature is set 1E+10. Finally, the reduction factor (λ) used in SA equal to 0.25, this value seems to be a reasonable number for the reduction factor.

RESULTS AND DISCUSSION
After satisfying all the constraints under the case of fixed well location, the distributions of the optimal values for all wells are shown in Table 4. The results are compared with the previous work obtained by Al-Aboodi, et.. al., (2013), which depends on genetic algorithm technique (GA) in their study on the same study area.
It is shown, all wells in SA are active and the optimal pumping rate values range from (1440 m 3 /day to 4000 m 3 /day ) while only thirty three active well in GA and the optimal pumping rate values range from (0 to 3840 m 3 /day). Table 4, the best maximum optimal value found by SA is (108320 m 3 /day), while GA found a better maximum optimal value is (91520 m 3 /day).The difference between the two optimal values is 18.35%. By combining of two different solutions, GA creates anew solution, while SA creates solution by modifying existing solution. Actually this modification makes the algorithm better or worse, is not straightforward, but depends on the problem and the representation. The distribution of hydraulic head based on the optimized pumping rates for first stress period is shown in Table 5.

CONCLUSIONS
Optimum management of groundwater flow using SA is carried out in this research, the optimum value of the objective function is (108320 m 3 /day). This value is then compared with another research based on GA. The results shown that the optimal value of SA is higher than the optimal value of GA by 18.35%. The difference between these values can be attributed to the population size is only one in SA compared to GA, This is a key difference between SA and GA. SA creates a new solution by modifying only one solution with a local move, while GA creates solutions by combining two different solutions. Actually this modification makes the algorithm better or worse, is not straightforward, but depends on the problem and the representation. Thus, it can be decided to recommend using SA approach in this topic. It is recommended for further works to establishing a network of monitoring wells to determine the hydraulic connection between the upper and lower aquifers, and establishing a three dimensional flow model in order to study the influence interaction between upper and lower aquifers.