Application of Simulated Annealing in Water Resources Management: Optimal Solution of Groundwater Contamination Source Characterization Problem and Monitoring Network Design Problems

[Extract] Estimating various characteristics of an unknown groundwater pollutant source can be formulated as an optimization problem using linked simulation-optimization. Meta-heuristics based optimization algorithms such as Simulated Annealing (SA), Genetic Algorithm (GA), Tabu Search etc. are now being accepted as reliable, faster and simpler ways to solve this optimization problem. In this chapter we discuss the suitability of a variant of traditional Simulated Annealing (SA) known as the Adaptive Simulated Annealing (ASA) in solving unknown groundwater pollutant source characterization problem.


Introduction
Estimating various characteristics of an unknown groundwater pollutant source can be formulated as an optimization problem using linked simulation-optimization. Meta-heuristics based optimization algorithms such as Simulated Annealing (SA), Genetic Algorithm (GA), Tabu Search etc. are now being accepted as reliable, faster and simpler ways to solve this optimization problem. In this chapter we discuss the suitability of a variant of traditional Simulated Annealing (SA) known as the Adaptive Simulated Annealing (ASA) in solving unknown groundwater pollutant source characterization problem. Growing anthropogenic activities and improper management of their impacts on groundwater quality has resulted in widespread contamination of groundwater worldwide. Coupled with ever increasing water demand leading to increased reliance on groundwater, it has resulted in a widespread recognition of public health risk posed by contaminated groundwater. This has triggered massive efforts for better management of groundwater quality in general and remediation of contaminated aquifers in particular. The sources of contamination in groundwater are often hidden and inaccessible. Characteristics of these pollutant sources such as their location, periods of activity and contaminant release history are often unknown. Groundwater contaminant source identification problem aims at estimating various characteristics of an unknown groundwater pollutant source using measured contaminant concentrations at a number of monitoring locations over a period of time. It has been widely accepted that for any remediation strategy to work efficiently, it is very important to know the pollutant source characteristic. A detailed account of different categories of source identification problems and various approaches to solve them has been presented in Pinder [26].
Typically, groundwater contamination is first detected by one or more arbitrarily located well or monitoring location. Unknown Groundwater source identification problem specifically attempts to ascertain the following source characteristics: Source type is often obvious. In some cases, information on groundwater contaminant source location may be available from preliminary investigations. If an exhaustive record of pollutant inventory and industrial activities of the area is available, it may be possible to infer start time. Release history of the source, however, is difficult to ascertain as the source is not physically accessible for measurements and hence it is unlikely that any accurate temporal record of contaminant fluxes released from the source exists.
Release history reconstruction problem is one of the most widely studied groundwater source identification problems. Ascertaining release history of the contaminant sources from available contaminant concentration measurements is an inverse problem as it requires solving groundwater flow and transport equations backwards in time and space. The process of solving this inverse problem is essentially the process of finding various unknown characteristics of source using observed information about the transport media and the effects caused by the source. In such circumstances, a solution cannot be guaranteed, especially when observed information is sparse. Even if the solution exists, it may not be unique. This is because different combinations of various source characteristics can produce the same effect at a monitoring location. Moreover, the solution of this problem is highly sensitive to measurement errors either in the observation data or model parameters and hence this problem has been classified as an ill-posed inverse problem [33] . When this inverse problem has to be solved by using inaccurate values of media parameters such as hydraulic conductivity and porosity and contaminant concentration observed at arbitrarily placed monitoring wells, it becomes even more challenging to obtain a reliable solution.
Methods proposed in the past to solve this ill-posed inverse problem can be broadly classified as optimization approaches, analytical solutions, deterministic direct methods and probabilistic and geo-statistical simulation approaches. A detailed review of these methodologies can be found in Atmadja & Bagtzoglou [2]; Michalak & Kitanidis [25]; Bagtzoglou & Atmadja [3] and Sun et al. [31,32]. The most effective of all suggested methods seems to be those based on optimization or probabilistic and geo-statistical simulation. Of the optimization methods, linked simulation-optimization approaches have been established as one of the most efficient methods. In this approach, a numerical groundwater flow and transport simulation model is linked to the optimization model. All the linked simulation-optimization approaches aim at solving a minimization problem with an objective function representing the difference in measured concentration and simulated concentration at various monitoring locations. The optimization model generates candidate solutions for various source characteristics. This is used as an input for the simulation model to generate estimated contaminant concentration observations at designated monitoring locations. The optimization algorithm then calculates the value of objective function by calculating the difference in contaminant concentrations estimated by simulation model and actual observed values at the same monitoring locations over a period of time. Over a number of iterations, the optimization algorithm minimizes the objective function value. Most prominent approaches in this category are linear programming with response matrix approach [9], nonlinear optimization with embedding technique [21][22][23], artificial neural network approach [28][29][30], constrained robust least square approach [31,32], classical optimization based approach [5][6][7], genetic algorithm based approach [1,24,27] etc.
In recent past, heuristic global search approaches such as Genetic Algorithm [12], Harmony Search, Tabu Search, Ant-Colony Optimization, Simulated Annealing [19] have developed rapidly and have been applied to a wide range of optimization problems. One of the major reasons for their popularity is the fact that these optimization methods do not easily get trapped in the local optima, thereby maximizing the probability of achieving a global optimal solution. Genetic algorithm (GA) and its variants, in particular, have been widely applied for solving unknown pollutant source identification. [11,24,27]. Genetic Algorithms are computational optimization algorithms that simulate the laws of natural genetics and natural selection and use it to search for the optimal solution.
Apart from GA or its variants, Simulated Annealing has also been used in solving inverse problems in groundwater management. Simulated annealing is inspired by the physical process of annealing in metallurgy which involves heating and controlled cooling of a material to reduce defects in crystal structure. The atoms are excited by heat and they become agitated while getting into higher energy states. The slow cooling allows a better chance for these atoms to achieve lower energy states than the ones they started with. In simulated annealing, a current solution may be replaced by a random "neighborhood" solution chosen with a probability that depends on the difference between corresponding function values and on a global parameter T (called temperature) that is gradually decreased in the process. Implementations of simulated annealing has been relatively limited because the traditional simulated annealing algorithm is reported to converge slower compared with GA or its variants. However, faster variants of simulated annealing have been developed and one of the most promising variants in terms of convergence speed is Adaptive Simulated Annealing (ASA) [14,15]. The ASA code was first developed in 1987 as Very Fast Simulated Re-annealing (VFSR) [13]. Ingber & Rosen [16] showed that VFSR is at least an order of magnitude superior to Genetic Algorithms in convergence speed and is more likely to find the global optima during a time limited search.
Linked simulation-optimization based approaches are computationally intensive as the simulation model has to be run many thousands of times before an acceptable solution is produced. This has been a deterrent to any desktop based implementation of the simulation-optimization approach. Faster convergence can reduce the computational burden significantly and thereby enhance the possibility of a desktop based implementation of linked simulation-optimization approach. This paper investigates the applicability of ASA to unknown groundwater contaminant source release history reconstruction problem and compares its performance to genetic algorithm based solution. The performance evaluation of competing simulation-optimization approaches are based on a realistic scenario of missing measurement data, where contaminant concentration measurements are available a few years after the sources have ceased to exist. Apart from the convergence speed, the two algorithms are compared for their ability to produce accurate source release histories with moderately erroneous data and with uncertainty in estimation of hydro-geological parameters. One of the most important factors that affects the execution time and accuracy of solutions generated by linked simulation-optimization approaches is the choice of observation locations. Poorly chosen contaminant observation locations often produce misleading results and hence it becomes important that after the initial estimation of the contaminant sources, a monitoring network is designed and implemented. In this study we use a monitoring network designed specifically to enhance the efficiency of source identification. However, a detailed discussion of the methodology used for monitoring network design is beyond the scope of this book.

Methodology
The linked simulation-optimization approach consists of two parts. An optimization algorithm generates the candidate solutions corresponding to various unknown groundwater source characteristics.
The candidate solutions are used as input in the numerical groundwater transport simulation model to generate the concentration of contaminant in the study area. The generated concentration at designated monitoring locations is matched to the observed values of contaminant concentrations at various time intervals at the same locations. The difference between simulated and observed concentration is used to calculate the objective function value which is utilized by the optimization algorithm to improve the candidate solution. The process continues until an optimal solution is obtained. A detailed schematic representation of this process of using SA as the optimization algorithm in a linked simulation-optimization model is presented in Figure 1. The classical simulated annealing (SA) algorithm has many associated guiding parameters such as the initial parameter temperature, annealing schedule, acceptance probability function, goal function etc. Effective application of the classical simulated annealing to a particular optimization problem normally involves a lot of trials and adjustments to achieve ideal values for all or most of these parameters. ASA, which is a variant of classical SA, helps overcome this difficulty to a certain extent by automating the adjustments of parameters controlling temperature schedule and random step selection thereby making the algorithm less sensitive to user defined parameters compared with classical SA. This additional ability of ASA combined with inherent ability of classical SA to find the global optimal solution even when multiple local optimums exists, makes it a natural choice for solving the groundwater pollutant source identification problem.

Governing equations
The three-dimensional transport of contaminants in groundwater can be represented by the following partial differential equation [17] Where C is the concentration of contaminants dissolved in groundwater,ML −3 ; t is time, T; x i , x j is the distance along the respective Cartesian coordinate axis, L; θ is the porosity of the porous medium, dimensionless; N is the number of chemical species considered; ∑ N k=1 R k is the chemical reaction term for each of the N species considered, ML −3 T −1 .
In order to solve this transport equation, linear pore water velocity needs to be known for the study area. Hence, it becomes necessary to first calculate the hydraulic head distribution using a groundwater flow simulation model. The partial differential equation for groundwater flow is given by the following equation: Where K xx , K yy , and K zz are the values of hydraulic conductivity (LT −1 ) along the x, y and z co-ordinate axes respectively; H is the potentiometric head (L); W is the volumetric flux per unit volume representing sources and/or sinks of water (T −1 ); S s is the specific storage of the porous media (L −1 ); and t is time (T).
The flow equation describes transient groundwater flow in three dimensions in a homogeneous anisotropic medium, provided the principal axes of hydraulic conductivity are aligned with the co-ordinate directions. A computer code called MODFLOW is used to solve this groundwater flow equation. MODFLOW was developed by United States Geological Survey (USGS) and is one of the most popular computer programs being used to simulate groundwater flow today. MODFLOW is based on modular finite-difference method which discretizes the study area into a grid of cells. The potentiometric head is calculated at the center of each cell.
To solve the three dimensional ground water transport equation, another computer code called MT3DMS is used. This is also a very popular computer program developed by the USGS and uses modular finite-difference just like MODFLOW. The transport simulation model (MT3DMS) utilizes flow field generated by the flow model (MODFLOW) to compute the velocity field used by the transport simulation model. [34]

Formulation of the optimization problem
It is assumed in this study that information on a set of potential source locations are available. The objective of simulation-optimization method then reduces to regenerating the source release histories at these potential source locations. Spatial and temporal contaminant concentration(C) is known at specific monitoring locations at various point of time. Candidate source fluxes are generated by the optimization algorithm. These values are used for forward transport simulations in MT3DMS. The difference between simulated and observed contaminant concentrations are then used to calculate the objective function. The objective function for this optimization problem is defined as: Where, cest k iob = Concentration estimated by the identification model at observation well location iob and at the end of time period k. nk = Total number of concentration observation time periods; nob= Total number of observation wells; cobs k iob = Observed concentration at well iob and at the end of time period k; w k iob = Weight corresponding to observation location iob, and the time period k.
The weight w k iob can be defined as follows: Where n is a constant, sufficiently large, so that errors at low concentrations do not dominate the solution [18]. It is possible to include other forms of this weight.

Optimization algorithms
Of the various simulated annealing implementations, it is evident in literature that the adaptive simulated annealing algorithm converges faster [16] while maintaining the reliability of results and hence it was preferred over traditional Boltzmann annealing implementation [19]. Its application to the unknown pollutant source identification has been limited but it is potentially a good alternative because its convergence curve is steep, thereby producing better results when execution time is limited.
Currently, the most widely used optimization algorithm for solving groundwater source identification problem using linked simulation-optimization model is Genetic Algorithm and its variants. The effectiveness of ASA in solving this problem is compared against the effectiveness of GA. Genetic algorithms (GAs) are population based search strategies which are popular for many difficult to solve optimization problems including inverse problems. GAs emulate the natural evolutionary process in a population where the fittest survive and reproduce [12]. GA-based search performs well because of its ability to combine aspects of solutions from different parts of the search space. Real coded genetic algorithm was used with a population size of 100, crossover probability of 0.85 and a mutation probability of 0.05. The values were chosen based on a series of numerical experiments.

Performance evaluation
In order to evaluate the performance of two different optimization algorithms involving comparison of solutions obtained, it is vital to first ensure that only one solution exists. In other words, a unique solution has to be guaranteed. This is possible only under the following idealized assumptions [33]: 1. The numerical models used for simulation of groundwater flow and transport are able to provide exact solution of the governing equations in forward runs.
2. All the model parameters and concentration measurements are known without any associated errors.
3. The unknown parameter is piecewise constant.
The first assumption is valid for cases where grid size and time step used in the numerical solution tends to zero. However since the groundwater simulation models used in this study have been proven to be stable and convergent, this assumption approximately holds. The second assumption however, cannot hold in real life scenarios. Hence it becomes necessary to use synthetically generated observation values initially which can be considered free of measurement errors. The third condition is implemented by assuming that the unknown fluxes are constant in every stress period. In such conditions it approximately resembles a well posed problem. Therefore these evaluations are initially carried out for synthetic data (simulated data) with known parameter values. There is another related issue of unique solutions. Whenever numerical simulation and optimization models are used, the convergence of the solutions may be another issue related to unique solutions. These issues are discussed in [4]. In this study the use of synthetic observation data, with known hydro-geologic parameter values reduces the ill-posed nature of the problem. The uniqueness of the solution cannot be guaranteed. However, sufficient iterations were allowed to ensure convergence to the optimal solution. Performance of the source identification methodology is evaluated using synthetic data from a three dimensional aquifer study area. The synthetic contaminant concentration data is obtained by solving the numerical flow and transport simulation models.

Simulating errors in concentration measurement data
Once the global optimal solution has been obtained for the idealistic assumption, the performance evaluation of developed methodology can take into account the effects of contaminant concentration measurement errors as well as uncertainty associated with the determination of hydro-geological parameters. To test the performance for realistic scenarios, concentration measurement errors are incorporated by introducing varied amounts of synthetically generated statistical noise in the simulated concentration values. The perturbed simulated concentrations represents erroneous measurements and is defined as follows: Where,

C ns = Simulated Concentration
S ud = a uniform random number between -1 and +1 a = a fraction between 0 and 1.0.

Performance evaluation criteria
The execution time of the algorithms is compared based on convergence curves which represent the value of objective function achieved versus time. To compare the ability of competing linked simulation-optimization approaches to produce accurate source histories, the errors in estimating source fluxes accurately is also used as a performance criterion. Normalized absolute error of estimation (NAEE) is used as the measure of errors in estimation of the sources. It can be represented as:

Discussion of solution results
The developed methodology was applied to a hypothetical illustrative study area with synthetically generated concentration measurements over space and time. Advantage of using a hypothetical study area lies in the fact that unknown data errors do not distort the performance evaluation of the methodology. This helps in understanding the drawbacks of developed methodology and improving it further.

Study area
The hypothetical study area is a heterogeneous aquifer measuring 2100m x 1500m x 30m and consisting of three unconfined layers as shown in Figure 2.
The East and west boundaries are constant head boundaries, whereas north and south boundaries are no flow boundaries. There are two sources (S1 and S2) of contamination. S1 is located in the top layer and S2 in middle layer. Five monitoring location (M1 through M5) are located in the first layer as shown in Figure 3. A grid size of 30m x 30m x 10m is used for finite difference based numerical calculation of groundwater flow and transport equations. Transport time step used for MT3DMS is 36.5 days. Other model parameters are listed in Table 1. Only a conservative contaminant is considered. There are two point sources of contaminants. One in the top layer and another one in the middle layer. A time horizon of 16 years is considered. Entire time horizon is divided into 5 different stress periods. The first four stress periods are each 1.5 years long and the final stress period is of 10 years duration. Sources are assumed to be active only in the first four stress periods or in the initial 6 years. Original source fluxes are presented in Table 2. It is assumed that groundwater contamination is detected at five different locations in the study area at the

Release history estimation with error free data
A set of error free observation data is generated. These observations are then used to evaluate the developed linked simulation-optimization methodology based on both GA and ASA. Input parameters used for GA and ASA are presented in Table 3. Every iteration of ASA based method uses one run of the groundwater transport simulation model (MT3DMS) whereas every generation of GA based method uses 100 (population size) runs of the same simulation model. Irrespective of the method, one run of the groundwater transport simulation model takes 4.281 sec to run on a Dell Optiplex®running an Intel®Core™2 Duo Processor at 2.93GHz. The execution time for one transport simulation run is however dependent on the computing platform. In order to keep the comparison independent of computing platform, both the methods were compared based on number of transport simulation runs which is directly proportional to the execution time. Both the methods were used to estimate source release histories using the error free data. In order to verify the convergence of each optimization method, time of run was made practically unconstrained. It was found that eventually both the optimization algorithms were able to achieve an objective function value very close to zero and identified the release history accurately. The objective function convergence profile as well as estimated fluxes were plotted at the end of 40,000 simulation runs of the groundwater transport model. Minimum value of objective function achieved is plotted against number of runs of the transport simulation model. The estimated flux values for both the sources in each stress period is also plotted against actual source fluxes. Convergence profile and source flux estimates are shown in Figure 4. Convergence profile shows that the objective function value for source identification model converges to a value very close to zero with about 5,000 simulation runs. However, further convergence is accelerated when using ASA algorithm. From these results, it can be concluded that the developed methodology is able to achieve optimal solution for an ideal error free scenario which resembles a well-posed problem.

Release history estimation with erroneous data
Five sets of erroneous observation data are generated with the formulation described in Equation 5. The value of fraction 'a' is specified as 0.1. These erroneous observations are used to reconstruct the release histories of contaminant sources. Linked simulation optimization method using ASA is compared with the one using GA as the optimization algorithm. Parameters used for both the optimization algorithms is presented in Table 3. Unlike the case with error free measurement data, in this case both the methods were used  Table 3. Parameters used in Optimization Algorithms to reconstruct source release histories using the erroneous data with a limit on execution time.
In order to make the comparison consistent by ensuring same number of simulation runs in the ASA and GA based methodologies, the number of simulation runs are restricted to 40,000. This restriction was based on the fact that increasing the number of simulation runs even to 80,000 resulted in very little improvement in the objective function value. Minimum value of objective function achieved is averaged over five solutions and is plotted against number of runs of the transport simulation model. The plot is presented in Figure 5. This Figure 5. Convergence Plot plot clearly shows that ASA based method converges much faster in the beginning and the GA based method is able to achieve comparable objective function values only after a much larger simulation runs. Because of the erroneous measurement data this problem may be ill-posed and the solution may not be unique. Therefore, lower objective function values do not always mean accurate reconstruction of the release histories.
In order to test the effectiveness of the competing methods based on accuracy of solutions produced, reconstructed release histories were compared to the actual release history after every set of 10,000 transport simulation runs. The results are shown in Figure 6. It can be seen that ASA based method is more efficient compared to GA based method after 10,000 and 20,000 simulation runs. However, as the execution time increases further with increase in number of simulation runs, the release histories produced by both methods become similar. This is also confirmed from the calculated values of NAEE presented in Table 4. As the execution time increases, the NAEE of ASA based method appears to increase only slightly. This could be due to statistical variation in the five different solutions and may be attributed to the input data error. Averaging over larger number of solutions may modify this inference. NAEE of GA based method consistently improves. However, the NAEE values obtained using ASA is still better in comparison.

Application of simulated annealing for monitoring network design
Monitoring network design in the context of groundwater quality management essentially means specifying the spatial location of monitoring wells and frequency of sampling. Since this is one of the most cost intensive part of most contaminated groundwater remediation problems, an efficient and cost effective design of monitoring network is essential. Monitoring of groundwater quality may be necessitated by a variety of objectives such as: 1. Unknown groundwater source characterization 2. Compliance monitoring for limiting the effects of groundwater contamination 3. Better aquifer characterization 4. Hydro-geological parameter estimation Irrespective of the various objectives, the problem of monitoring network design can be formulated as an optimization problem [8,20]. While designing a monitoring network for estimating unknown groundwater source characteristics, the objective of optimization can be to maximize the reliability of estimated source characteristics or to minimize the total number of monitoring locations in the network or both. Compliance monitoring is aimed at minimizing the area of contamination when the contamination is first detected at monitoring network or maximizing the probability of detection of contaminant in groundwater. Often, only the average values of hydro-geological parameters of the aquifer are known. This results in uncertainty in the modeling results. In order to better characterize an aquifer, spatial distribution of hydro-geological properties should be specified. This objective can be achieved by sampling hydro-geologic parameter at sufficient locations such that the interpolated values can represent actual hydrological parameters accurately. The objective of optimization in this case is to find the minimum number of samples required to accurately represent a population of random hydro-geological parameter values. In all such cases, Adaptive Simulated Annealing can be efficiently used as the tool for optimization. Our attempts to develop classical simulated annealing algorithm for optimal design of a dedicated monitoring network for enhancing the efficiency of source identification was successful to a large extent. However, the mixed integer nature of the decision variables in a monitoring network design problem makes the application of classical simulated annealing algorithm a bit constraining. Adaptive Simulated Annealing is more suitable to solve this monitoring network design problems.

Suitability and sensitivity of adaptive simulated annealing
In the application discussed here, Simulated Annealing is utilized for finding the global minimum of a cost function that characterizes large and complex systems such as transport of pollutants in groundwater.Simulated Annealing, as an algorithm, is very efficient in solving non-convex optimization problems by ensuring that it does not always move downhill on a complex non-convex search space and hence avoids getting trapped in local minimum. Simulated annealing also differs significantly from conventional iterative optimization algorithms in that gross features of the final state of the system are seen at higher temperatures whereas the finer details of the state appear at lower temperatures [10]. The fact that simualted annealing ensures a global optimal solution enhances its suitability for solving ill-posed inverse problems in general and the problem of unknown groundwater pollutant source characterization in particular.
Its ease of use and remarkable efficiency in handling complex objective functions and constraints has made simulated annealing an attractive choice for solving a wide range of complex optimization problems. However, the slow convergence and hence long time of execution of standard Boltzmann-type simulated annealing has been a constraint. Adaptive Simulated Annealing removes that constraint by making the annealing schedules decrease exponentially in annealing-time, thereby making the convergence much faster. A major difference between ASA and traditional Boltzamnn Annealing algorithms is that the ergodic sampling takes place in terms of n parameters and the cost function. In ASA the exponential annealing schedules permit resources to be spent adaptively on re-annealing and on pacing the convergence in all dimensions, ensuring ample global searching in the first phases of search and ample quick convergence in the final phases [15].
Another major advantages of using Adaptive Simulated Annealing is also the fact that the parameters of algorithm are adjusted adaptively and hence the solutions do not vary widely if parameter values are changed within reasonable limits. This is in contrast with Genetic Algorithm where even minor changes to parameters such as mutation probability, cross over probability or population size causes a significant difference in the solutions.

Conclusion
A linked simulation-optimization method for source identification was developed based on adaptive simulated annealing. It was applied to an illustrative study area. The results obtained were compared with those obtained using genetic algorithm, a more commonly used optimization approach. It is evident from the limited numerical experiments that adaptive simulated annealing algorithm based solutions converge to the actual source fluxes faster than genetic algorithm based solutions. This results in substantial saving in computational time. The source fluxes identified by using adaptive simulated annealing are closer to actual fluxes when compared to the results obtained using genetic algorithm, even when the observation data are erroneous and the hydro-geological parameters are uncertain. It can be concluded that adaptive simulated annealing is computationally more efficient for use in simulation-optimization based methods for identification of unknown groundwater pollutant sources, specially in a time constrained environment. Use of ASA has the potential to reduce CPU time required for solution by an order of magnitude. However, with very large number of iterations in the linked simulation-optimization approach, it is possible that the solutions obtained using GA could converge to a marginally better solution compared to that ASA based algorithm. However, it appears that ASA based solutions converge very close to the optimal solution using only a small fraction of iterations required while using GA. The relevance of contaminant monitoring locations is demonstrated. Further studies are required to develop dedicated monitoring networks which can increase the efficiency of source identification.

Author details
Manish Jha and Bithin Datta James Cook University, Townsville and CRC CARE, Adelaide, Australia