Received 10 November 2015; accepted 10 January 2016; published 13 January 2016
Groundwater contaminant source characterization is the first step of effective groundwater remediation and water resource management. Limited and sparse pollutant concentration measurements data are the main challenge for accurate pollution source characterization precisely in terms of location, magnitude, and duration of activity. The potential pathways for pollution can be determined based on contaminant source characterization. However, adequate measurement data are essential for accurate source identification; budgetary constraints also restrict long-term and widespread spatiotemporal monitoring. In addition, the arbitrary location of well or group of wells may not characterize unknown groundwater pollution sources accurately. Therefore, using an optimal monitoring network to obtain reliable and efficient data is necessary to identify groundwater sources of pollution and describe the aquifer flow and contamination transport precisely. The aim of this paper is to design optimal monitoring network based on initially available sparse contamination measurement data in an illustrative highly complex aquifer such as an abandoned mine site area.
The selected objectives, the aquifer parameters, and specific conditions affect the optimal monitoring network design. Some of the initial contributions towards methodologies for monitoring network design include: designing an optimal monitoring network for early detection of contamination   , locating wells in a monitoring network under conditions of uncertainty  , minimizing the number of monitoring wells  , and reducing the cost of monitoring for groundwater quality monitoring   . Ref.  discussed monitoring network design in terms of the risk, cost, and benefit in a landfill site. Ref.  proposed an optimization based methodology for determining the optimal sampling plan for groundwater quality monitoring. Most of the earlier developed methods were aimed at finding optimal monitoring locations for early detection of pollution along with reducing the cost of monitoring  without considering source characterization.
Some of the more recent design methodologies for monitoring network design include: long-term groundwater monitoring network design for multiple objectives  , cost-effective sampling design  , minimization of temporal redundancy and maximization of spatial accuracy   , sampling strategy in space and time using Kalman filter  , and using monitoring network design for source identification and redundancy reduction with feedback information  - . Ref.  proposed large, long-term groundwater monitoring (LTM) design problems using a new multi-objective evolutionary algorithm (MOEA), and also Ref.  reported a variation of dynamic monitoring network design methodology. The evolutionary optimization algorithms have increased the efficiency and accuracy of optimally detecting contaminant transport processes with objectives of minimizing the number of observation wells. The main optimization algorithms include genetic algorithm (GA)  -  , simulated annealing (SA)  - , adaptive simulated annealing (ASA)  , and genetic programming (GP)   - . Kriging as a geostatistical interpolation technique has been used in groundwater monitoring network design     . Other monitoring network design methodologies incorporated optimal sampling locations for plume detection  , and also reduced the cost in a groundwater quality monitoring network  .
Refs.    used Genetic Programing (GP) based optimal monitoring network design to improve the efficiency and accuracy of pollutant source identification. Potential applicability of GP in groundwater problems was discussed by Ref.  . These methodologies used trained GP models to calculate the impact factor of the sources on the candidate monitoring locations. However, these monitoring networks were static in nature. A number of methodologies have been proposed using different optimization algorithms for improving the monitoring network design results as reported in   .
In this study a linked optimization based approach is developed with three different objective functions: 1) Minimizing the solute mass estimation error which uses a coupled simulated annealing (SA) and kriging routine; 2) Minimizing the concentration gradient based optimization model developed (locate the plume boundary); 3) Selecting the monitoring locations with (potentially) highest concentration. Performance of its efficiency and reliability is evaluated in a contaminated aquifer resembling an abandoned mine site.
Design of the optimal monitoring network for the three different objectives of design was determined by utilizing a developed software package, CARE-GWMND  . The software package provides ease of applicability in complex real life cases and portability of the methodology to multiple scenarios of optimal monitoring network design. Performance of the developed methodology incorporated in the software is demonstrated by solving a problem of optimal observation wells design in a contaminated aquifer resembling an abandoned mine site in Queensland, Australia.
Potential monitoring locations are first specified as the candidate locations for selecting the optimal monitoring well locations. The observed pollutant concentration measurements available from the initial arbitrary located wells are interpolated for the entire study area to estimate the pollutant concentration for all the grid locations. The interpolation is done using Ordinary or Simple Kriging  . A Simulated Annealing (SA) based optimization algorithm is used to solve the optimization problem. Three different optimization models have been formulated to select the well locations which meet a specified objective 1) minimize the solute mass estimation error which uses a coupled SA and kriging routine  , 2) concentration gradient based optimization model developed by  and 3) selecting the monitoring locations with (potentially) highest concentration.
Kriging as a geostatistical interpolation technique provides best linear unbiased prediction of the intermediate unknown data points between known data points  . This method has been popular in estimating hydrogeological parameters, hydraulic head and pollution concentration over large areas by using a relatively small number of observations  .
The first step in kriging application is constructing an experimental semi-variogram and fitting it into a standard model such as spherical or exponential. The semi-variogram describes the relationship between the variance in data values and the distance between data points by plotting the variance against distance as shown in Figure 1. A kriged estimate is a weighted sum of the sampled data values around it.
Figure 1. Typical semivariogram structure.
2.2. Groundwater Flow and Transport Simulation Model
A three-dimensional numerical model MODFLOW  is utilized for simulating the flow process in the selected aquifer study area. The MODFLOW model aims to simulate the steady and transient flow. As several variables affect the groundwater flow, it is generally described as a partial differential equation in space and time. The partial differential Equation (1) for groundwater flow is given by Ref.  :
Kxx, Kyy and Kzz represent the values of hydraulic conductivity along the x, y and z coordinate axes (LT−1);
h is the potentiometric head (L);
W is the volumetric flux per unit volume representing sources and/or sinks (T−1);
Ss is the specific storage of the porous material (L−1);
t is time (T);
x, y and z are the Cartesian co-ordinates (L).
MT3DMS  is a modular three-dimensional transport model in groundwater system which is used to simulate the transport process in terms of advection, dispersion, and chemical reactions of dissolved constituents in groundwater flow systems under general hydrogeological conditions. MT3DMS computes the pollutants transport processes using the MODFLOW flow field results. Equation (2)  described three-dimensional transport of pollutants in groundwater system using the partial differential equation:
C is the concentration of pollutants dissolved in groundwater (ML−3);
t is time (T);
xi is the distance along the respective Cartesian coordinate axis (L);
Dij is the hydrodynamic dispersion coefficient tensor (L2T−1);
vi is the seepage or linear pore water velocity (LT−1); it is related to the specific discharge or Darcy flux through the relationship, vi = qi/θ;
qs is volumetric flux of water per unit volume of aquifer representing fluid sources (positive) and sinks (negative) (T−1);
Cs is the concentration of the sources or sinks (ML−3);
θ is the porosity of the porous medium (dimension less);
Rk is chemical reaction term for each of the N species considered (ML−3T−1).
2.3. Simulated Annealing as Optimization Algorithm
Simulated annealing (SA)  is used as the optimization algorithm for approximating the global optimum of a given function. The basic concept of SA is the physical process of heating a material and then slowly lowering the temperature to decrease defects, thus minimizing the system energy. At each iteration of the SA algorithm, a current solution replaces a random nearby solution based on the difference between the corresponding function values and algorithm control parameters (temperature reduction rate, initial temperature etc.). SIMANN a SA code in FORTRAN has been developed by Ref.  is utilized for the optimization algorithm application in this study.
2.4. Monitoring Network Optimization Problem Formulation
Different criteria may be used for the selection of monitoring well locations, depending on the nature of the problem. Different optimization models can be formulated to achieve different objectives of monitoring network design. Three main objectives of selecting the best monitoring well locations have been formulated and implemented in this study. The details of the optimization models are given below:
1) Minimize the pollutant mass estimation error
This optimization model is designed to select the monitoring locations which provide the best estimate of pollutant mass present in the aquifer, given a set of potential observation well locations and a specified maximum number of wells that can be implemented. To calculate the mass estimation error, the formulation requires the actual mass present in the model domain as an input parameter (M), and the solution is dependent on M. We estimate the actual mass based on the initial observations from arbitrary located monitoring wells and the mass balance estimates from the MT3D model. Total estimated mass at the end of the simulation period = initial mass calculated by interpolating the initial observations plus net increase of mass in the study area obtained from the MT3D simulation results.
An integrated SA optimization and kriging code which dynamically calls a kriging sub-routine is used to solve the non-linear binary programming problem  . Kriging interpolates the concentration to obtain the total mass estimation using the simulated concentration data at observation well locations. SA algorithm minimizes the mass error between actual total mass and total mass estimated by optimizing the location and number of monitoring wells using the objective given by Equation (3):
2) Concentration gradient based optimization model (to locate the plume boundary)
To locate the plume boundary, a mathematical programming model is formulated to choose the monitoring well locations that lie on the boundary of the plume. The underlying logic behind this formulation is that the concentration gradient is steep closer to the source and becomes flat closer to the boundary of the plume  . An optimization problem is formulated using the above logic to find locations with minimum concentration gradient.
3) Selecting the monitoring locations with (potentially) large concentrations
The formulation for selecting the well locations with large concentrations is given below. It should be noted that an optimization model is solvable without the application of a formal optimization algorithm.
n is: potential monitoring well locations;
i is: rows in the model grid;
j is: columns in the model grid;
i,j is: grid cells in the model domain.
is simulated concentration at potential well location n.
M is total contaminant mass present in the study area (actual mass).
NMax is Maximum number of wells to be installed.
Pi,j is porosity of cell i,j.
Vi,j is volume of cell i,j.
Xn is 1 if well location n is selected, 0 otherwise.
CKi,j is interpolated (kriged) concentration in cell i,j. Interpolation is performed based on Xn and, i.e., if Xn is 1, then will be a data (sampled) value for interpolation.
MEST is Total mass estimated based on Xn.
in, jn is cell (row, column) reference for potential monitoring well location n.
is simulated concentration in cell i,j.
Gn is sum of the magnitudes of the four concentration gradients around i,j calculated as Equation (5).
3. Performance Evaluation
Performance of the proposed optimal monitoring network design methodology is evaluated utilizing a polluted aquifer study area, hydrogeologically resembling the Mount Morgan abandoned mine site in Queensland, Australia. In this study the polluted aquifer region is referred to as the “impacted area” and the total aquifer region considered  in this study is called the “aquifer study area” (Figure 2). These data are also utilized to calibrate the numerical simulation model developed using MODFLOW.
3.1. Polluted Aquifer Site Description
The historic mine site utilized as an illustrative contaminated aquifer in this study is located in central Queensland, Australia (Figure 2). This copper, gold and silver mine was operational for nearly 100 years between 1889 and 1990. Dee River forms a natural boundary on the southern side of the study area which inflows into Don and Dawson River and finally into the Fitzroy River. It was one of the largest gold mines in Australia. During its operational life time, approximately 274 tons of gold, 37 tons of silver, and 387,000 tons of Copper were mined from underground and open cut operations  over 100 years mining activities (1889-1990). Acid Rock Drainage (ARD) and Acid Mine Drainage (AMD) generated from the flooded open cut and waste deposits and their tailings, are likely contamination sources at the mine site, and the Dee River downstream of the mine. These distributed pollution sources affect groundwater quality as well as surface water quality. Although several remediation measures were undertaken for this area, identification of contaminated sources and characterization of ARD and AMD chemical sources and pathways and physical behaviors in the groundwater system are required to achieve effective and efficient remediation strategies and optimal groundwater management.
Figure 2. Plan view of the mount Morgan study area and the impacted area.
3.2. Mount Morgan Mine Site: Environmental Issue
During the early years of the mine’s operations, when mining activity was at its height, there was little or no environmental regulation. Mining activities have left a legacy of impacts on the Dee River and the site itself. Indeed, today the environmental pollution is caused by acid and metalliferous drainage also known as acid mine drainage. This pollution is generated when sulfidic rocks, such as pyrite, are exposed to water and oxygen in the surface environment rather than being isolated as crystalline rock underground  . This phenomenon forms sulfuric acid, which in turn dissolves extreme concentrations of salts and metals, including potentially copper, arsenic, nickel, cadmium, zinc, aluminum, iron and many more. Acid mine drainage is a deadly toxic soup to aquatic ecosystems and biodiversity  .
Starting in 1982, tailings were dredged from Sandstone Gully and treated using the Carbon-In-Pulp (CIP) process before being backfilled into the open cut. The CIP process is an extraction technique for recovery of gold which has been liberated into a cyanide solution as part of the gold cyanidation process. This process uses activated carbon to capture the gold tailings. After final closure in 1990, the partially backfilled open cut (and Sandstone Gully) were allowed to flood by natural inflows (surface runoff and groundwater inflow)  .
The main part of waste rock from the open cut is considered as acid-forming. It’s located on the depth of weathering. Furthermore, it contains up to 10% sulphur with sulphide minerals such as pyrite (FeS2, also called fool’s gold). The term “Mundic” describes a copper ore that begins to be melted. The Mundic tailings were placed into the historic drainage channel of Mundic Creek and the other tailings were placed into tailings dams.
Geochemical tests carried out by Ref.  show that the Mundic Red tailings were unreactive whereas the Mundic Grey tailings are highly reactive. These could release sulphate (), iron (Fe), aluminum (Al) and copper (Cu). Thus, if tailings were deposited without proper containment, this would be a real environmental issue   .
Therefore, it is evident that mine activities generate many pollutants which can cause environmental issues. This study serves a very limited scope, i.e., to examine the general applicability of monitoring network design methodologies for designing different optimal networks based on different specified objectives of design in a contaminated abandoned mine site area. This study only serves a demonstration purpose, which includes the implementation of the flow and transport simulation model for a hydrogeologic contaminated site with various sources of contamination. The model calibration process is carried out for the flow modeling only, due to lack of adequate concentration measurement data for calibration. Therefore, the contaminant transport process only models a typical advective dispersive contaminant transport process for a typical conservative contaminant. The spread and detection of such a typical contaminant in the study area hydrogeologically resembling the Mount Morgan mine site in Queensland, Australia. The study site is extensive and complex area (this will be detailed later). Fir this preliminary study, only one type of pollutant which does not undergo reactive process was considered. Therefore, our study will focus on one typical unknown conservative pollutant (such as cobalt, silver etc.).
The climate at the site is seasonal, with average maximum monthly temperatures ranging from 32˚C in January to 23.1˚C in July. The average annual rainfall is approximately 815.5 mm with a large amount of the annual rainfall occurring during the wet season from November until March. This graph (Figure 3) shows the repartition of the rainfall during one year. We can observe that the rainfall is less during the Australian winter (between 22 and 43 mm). While, from October to March, the rainfalls can exceed 100 mm.
The geology of the site has been described in detail in    . Devonian Mount Morgan Tonalite which has almost no permeability is the major ore in this site. Devonian rhyolitic Mine corridor Volcanic have provided Au-Cu mineralization, which is surrounded and intruded by the tonalite widely. Generally, rock formation in Australia are defined without primary permeability. Fractures and faults control any secondary permeability. Moreover, some dykes compartmented the area, and causes groundwater discharge deeper from the mine site  .
Pyrite is the main ore mineral in this site, which is followed by chalcopyrite and trace pyrrhotite, sphalerite, wurtzite, magnetite and Au  . The majority of waste deposits comprise of quartz-pyrite waste rock extracted
Figure 3. Rainfall and temperature during a typical one year at mount Morgan  .
from the mining activities. Sulphide oxidation, pyrite waste rock exposition, and also tailings on the site generate hazardous acidic and metal-rich drainage  . The absence of carbonate rocks to naturalize groundwater pH has exacerbated this situation  .
3.5. Groundwater Flow Modelling and Calibration
Groundwater flow and reactive chemical contaminant transport processes are simulated with specified hydrogeological parameters as given in Table 1. Concentration measurement data at existing concentration monitoring locations are used for calibrating the flow and transport numerical simulation model. In actual application not directed towards performance evaluation of the developed methodology, the actual concentration measurement data available at specified observation locations of the study area aquifer are to be utilized as specified observed concentration inputs to the optimal source characterization decision model  .
The ground topography generally slopes from north east towards the river in south. The Mount Morgan deposit is situated in the Calliope Block, which is an important tectonic zone and occurs along the eastern margin of Australia from Rockhampton to Warwick  . The geology Mount Morgan is mainly constituted of volcanic and metamorphic rocks (quartz, feldspar, basalt etc.). All of the rock formations are considered to have no primary permeability and any secondary permeability is believed to be controlled by structure (fractures). The groundwater flow occurs therefore mainly in the permeable mine waste dumps and in shallow bedrocks that have been fractured by mining and have become more permeable. The area consists of 4 layers  : waste rock dumps and tailing; highly weathered bedrock; partially weathered bedrock; and tight bedrock
For each layer, porosity and hydraulic conductivity are considered as homogeneous. These parameters were estimated by field investigation conducted by Ref.  . The surface area is about 7.7 km2 with an open cut or pit which is over 2.5 km long and over 300 meters deep (Figure 4). The geographical coordinates of this area are 23˚38'28"South and 150˚22'31"East. Figure 4 shows the study area. The Orange curve indicates mountainous boundaries of the study area, and the blue lines indicating the Dee River. At the center is the erstwhile open cut mine, shown in turquoise. So, the model domain is bounded by topographic highs to the west and south and the Dee River and Arnolds Gully to the East and North, respectively. The hydrogeologic parameters of the study area are summarized in Table 1 and are based on the study reported in Refs.   .
The main aim of calibrating the flow model was to assign proper head boundary conditions, and estimate the spatially and temporally averaged recharge. The calibration was partially based on the earlier flow model
Figure 4. Study area resembling the mount Morgan mine site in Queensland, Australia.
Table 1. Hydrogeologic parameters used in flow and transport model of the study area.
calibration for the same study area as reported in  . The evaluation of the calibration was based on limited and sparse measured hydraulic head field data available for the site. Hydraulic head measurement data from 11 observation locations spread across the impacted area were used for simulated model calibration. The hydraulic head data used for model calibration and validation were recorded in 2002 and 2003  . A portion of the average rainfall intensity per year is specified as a recharge for calibration of the simulation model also using measured head data from observation locations. Calibration targets for the developed model were set to be within two meter intervals of the observed hydraulic head value in observation locations with a confidence level of 90 percent. The model boundary conditions were manually adjusted to achieve the calibration targets. The measured and simulated heads were compared at selected location. In Figure 5, the green bars signify that calibration target is achieved for the shown well numbers, whereas yellow bar represents intermediate errors. This calibration process was continued until acceptable calibration results were obtained. More rigorous calibration can be undertaken if a large number of spatiotemporal field measurement data is available.
3.6. Pollutant Transport Simulation in the Impacted Area
The three-dimensional transient transport simulation model, MT3DMS, simulates the pollution transport process
Figure 5. Flow simulation and calibration results. (a) Simulated total head; (b) Map of the total head calibration errors for the study area; (c) Comparison of observed and simulated total heads for the study area at observation wells locations.
originating from a specified source. For the purpose of implementation, a specified pollutant concentration is assumed at the pit. The transport model uses the flow field generated by the flow model to predict the movement of the pollutants in the impacted area of the aquifer over 5 years. The estimated data at the potential well locations are used by linked kriging―optimization to obtain the objective functions. Figure 6, Figure 7 show the pollution distribution in the study area after 5 and 10 years from the initial date, respectively.
3.7. Monitoring Network Design and Performance Evaluation
In the monitoring network design, the potential well locations are selected to be used as input which used in optimization model to design the optimal monitoring networks. In this study thirty potential wells are selected in the study area (Figure 8). The transport simulation model estimates the concentration of the pollutant in the selected wells and the total mass estimation obtained using Kriging based interpolation throughout the aquifer. Majority of the potential wells are around the pit and the river. The total actual mass concentration is estimated using interpolation of the pollution concentration data in available arbitrary observation wells. Three scenarios for designing the monitoring network are considered using three different objective functions to design the optimal monitoring network. The performance of the methodology can be evaluated in terms of how the designed well locations are able to achieve the three main objective functions optimally.
4. Optimal Monitoring Network Design Using CARE-GWMND
CARE-GWMND  is a software package that is used for the optimal monitoring network design. The CARE-GWMND software package was developed in James Cook University, Australia in 2013 under the guidance of the first author, with a research funding provided by CRC-CARE. The software package provides a user friendly interface for optimal monitoring network design. The software essentially comprise of excel user interface for input data, backend VBA routines for generating the input file for simulation and optimization, and FORTRAN executables for solving the source identification problem. SA is used as solution algorithm for solving the optimization problem. The schematic structure of the software is shown in Figure 9. All the relevant details from the calibrated flow, transport and interpolation models are input into the software package as shown in Figure 10. The optimization parameters are chosen and the well location design model is executed from the user interface. The total number of wells and the objective function can be selected in the interface. Once the monitoring network design model has finished execution, the results are exported back to the user interface.
5. Results and Discussion
The performance evaluation results of best (optimal) well locations obtained as solution by using the monitoring network design model are presented in Figure 11(a) and Figure 11(b), Figure 12 for different objective functions. The optimal well locations which minimize the contaminant mass estimation error takes the longest time to run. Nevertheless, we obtained the following result shown in Figure 11(a). The well locations which have the minimum concentration gradient subject to a lower bound are shown in Figure 11(b). These results provide evidence of the capability of the model to locate the plume boundary correctly.
Figure 6. Simulated concentration of the pollutant at mount Morgan area after five years.
Figure 7. Simulated concentration of the pollutant at mount Morgan study area after ten years.
Figure 8. CARE-GWMND, potential monitoring locations in layer 2.
Figure 9. Schematic diagram of CARE-GWMND software. (Datta et al., 2014C  ).
Figure 10. A typical screen shot while running CARE-GWMND.
Figure 11. CARE-GWMND solution, optimal monitoring locations in layer 2 for objective functions 1 and 2 pics (a) and (b) respectively.
Figure 12. CARE-GWMND solution for optimal monitoring locations in layer 2 for objective 3.
Figure 11(b) illustrates the plume boundary with the well locations for which we have the minimum concentration gradient summation. It is noted that most of the points are located near the boundaries with the exception of two points which are the locations 10 and 21. These locations may have been chosen due to low concentration in these locations, and the flat gradient of the concentration values in these zones. These results are consistent with the simulation results of the pollutant plume after five years of simulation. Indeed, in the study area, the gradient is oriented from North West to South East. Therefore, it clearly explains why the majority of optimal monitoring locations are located near the boundaries of the mine site and in the vicinity of the pollutant plume boundary.
The well locations which represent the maximum of pollutant concentration are shown in Figure 12 which located near the open cut. The optimization model selected all the optimal monitoring locations around the pit which is our single main source of pollution. Hence, these are likely locations of large concentrations, and therefore also chosen by the monitoring network design model with objective function 3. Therefore, these preliminary monitoring network designs obtained as the solution of the optimal design methodology, demonstrate that the developed methodology and the developed software can be applied to a very complex contaminated aquifer study area such as an abandoned mine site. The solution results also appear reasonable, given the physical and geochemical characteristics of the study area.
The modeling of the flow and transport processes in a complex contaminated aquifer resembling a real-life abandoned mine site in Queensland, Australia is demonstrated. The developed methodology along with the software incorporating the methodology for optimal design of monitoring networks appears to perform satisfac- torily in designing useful and efficient monitoring networks for a hydrogeologically complex contaminated aquifer site, resembling the contaminated abandoned mine site. The performance evaluation results show that the developed methodology is successful in developing an optimal monitoring plan for the contaminant concen- tration values in a highly complex groundwater pollution scenario. The developed software package CARE- GWMND  is effective in finding the solution results for optimal monitoring network design, and it may be possible to adopt the methodology and the software to various other complex real life scenarios of groundwater pollution. Further studies for incorporating more complex reactive multiple species geochemical processes in the optimal monitoring network design for a mine site are being carried out now.
The first author thanks, CRC-CARE, Australia for providing financial support for this research through Project No. 220.127.116.11.09/10 (2.6.03), CRC-CARE-Bithin Datta, at James Cook University, Australia.