Sustainable groundwater modeling using single- and multi-objective optimization algorithms

Modeling groundwater behavior is among the important processes hydrogeologists have been trying to quantify for a long time in order to address emerging groundwater problems. Simulation of groundwater system, especially in arid region, is difficult due to the complex and multi-objective nature of the groundwater system. In an era of increasing demand for groundwater, new paradigms are thus urgently required to devise innovative decision-making tools and optimize the drawdown of aquifer system.

Groundwater simulations have been typically performed using simulation-optimization algorithms (e.g. McKinney & Lin 1994; Wang & Zheng 1998; Giustolisi & Simeone 2006; Wu & Zhu 2006; Zhu et al. 2006; Giustolisi et al. 2008; Hamraz et al. 2015). These models have been used to solve design and operation problems associated with groundwater hydraulic control, water supply, and remediation (e.g. McKinney & Lin 1994). Practical experience with the calibration of groundwater processes suggests that single-objective function, no matter how carefully chosen, is often inadequate to properly measure all aspects and characteristics of groundwater system. Because sustainable management is necessarily a multi-objective problem, no optimal solution in the traditional sense can be found, and decision makers should express their preferences within a set of non-dominated solutions (e.g. McPhee & Yeh 2004). One strategy to explicitly recognize the multi-objective nature of groundwater system is to define several optimization criteria (objective functions) that measure different (complementary) aspects of the system behavior and to use a multi-objective optimization method to identify a set of non-dominated, efficient, or Pareto optimal solutions (e.g. Gupta et al. 1998; Yapo et al. 1998; Boyle et al. 2000; Giustolisi et al. 2008). The multi-objective approach is based on the Pareto dominance criterion and an evolutionary strategy employs to solve the combinatorial optimization problem (Giustolisi & Simeone 2006). The Pareto solutions represent tradeoffs among different incommensurable and often conflicting objectives, having the property that moving from one solution to another results in the improvement of one objective while causing deterioration in one or more others (Vrugt et al. 2003).

A large body of literature related to the application of either deterministic or stochastic optimization methods exists in groundwater studies. Deterministic optimization methods include linear programming, nonlinear programming, and dynamic programming. These methods were adopted by several scholars (e.g. Aquado & Remson 1974; Gorelick et al. 1979, 1984; Willis 1979; Wanakule et al. 1986; Jones et al. 1987; Andricevic 1990; Lee & Kitanidis 1991; McKinney & Lin 1994; Giustolisi & Simeone 2006; Giustolisi et al. 2008). The second type of optimization methods were often referred to as heuristic algorithms includes genetic algorithm (GA), particle swarm optimization (PSO), shuffled complex evolution, simulating annealing, etc. These methods were used in highly nonlinear and multimodal problems, under various complex constraints (e.g. McKinney & Lin 1994; Wang & Zheng 1998; Wu et al. 1999; Giustolisi & Simeone 2006; Wu & Zhu 2006; Zhu et al. 2006; Giustolisi et al. 2008; Ayvas 2009; Gaur et al. 2011a, 2011b; Mategaonkar & Eldho 2012; El-Ghandour & Elsaid 2013) to optimize groundwater characteristics in different regions.

In the literature, many efforts have dealt with the multi-objective optimization of groundwater management problems. For instance, Wang & Zheng (1998) coupled a GA and simulated annealing with the modular three-dimensional finite-difference (MODFLOW) groundwater model, to optimize groundwater remediation design. They applied various simulation periods and included both fixed and operating costs. Park & Aral (2004) presented a multi-objective optimization approach to determine well locations in the coastal regions, which maximized pumping rates and minimized the distance between critical stagnation point and the reference coastline location. Reed et al. (2001) introduced a multi-objective approach to cost effective long-term groundwater monitoring using an Elitist non-dominated sorted GA. Giustolisi et al. (2008) proposed a methodology to determine a dynamic relationship between rainfall and water-table depth for a shallow unconfined aquifer system located in southeast Italy. Siegfried et al. (2009) presented a multi-objective evolutionary algorithm to optimize the placement and the operation of pumping facilities over time. Saafan et al. (2011) coupled a multi-objective GA optimization model with MODFLOW and maximized the pumping rates in El-Farafra oasis, Egypt. They predicted the maximum pumping rates and minimum operation cost as well as the prediction of the future changes in both variables.

The aforementioned algorithms have focused on using a single GA approach to compute groundwater characteristics that may not be applicable for large-scale groundwater system. When dealing with large scale and complex system, numerous conflicting objectives can arise so that the number of decision and state variables may increase rapidly with the scale of the problem (e.g. McPhee & Yeh 2004). In such circumstances, single objective techniques may provide unsatisfactory results to decision makers, and therefore multiple optimization solutions should be sought. In many optimization domains the solution of the problem can be multi-dimensional and can be only computed simultaneously by assembling a hierarchy of multiple optimization algorithms.

This study presents the first attempt that has linked a multi-algorithm, genetically adaptive multi-objective (AMALGAM) model, proposed by Vrugt & Robinson (2007), with MODFLOW to optimize pumping rates within a well distributed set of Pareto solutions. The purpose of this research was to find optimal solution to cover maximum demand for arid groundwater management and optimized the pumping rates using three minimization objectives, i.e. minimizing shortage affected by the failure to supply, Modified Shortage Index (MSI) and minimization of extent of drawdown within pre-specified Pareto regions. Focusing on these objectives, AMALGAM simultaneously merges the strengths of the covariance matrix adaptation evolution strategy, GA, differential evolution and particle swarm optimization (PSO) to ensure a reliable and computationally efficient solution is achieved for multi-objective optimization problems.

This research first coupled MODFLOW to an advanced swarm-intelligence-based algorithm (i.e. Cuckoo optimization algorithm (COA); Rajabioun 2011) to compute aquifer hydrodynamic parameters using automatic search method. The ability of the COA approach led to more accurate estimates of error variance and helped us characterize groundwater mechanisms at each subscale (zone). We optimized aquifer properties by formulating an optimization method within the groundwater numerical model and analyzed predictive error in a highly parameterized model. At the second step, this study used a hierarchy of multi-algorithm designed in the AMALGAM model to define optimal Pareto solutions for groundwater drawdown in order to ensure sustainable groundwater management over an arid region in east Iran.

The bedrock of the aquifer is mostly hard rock, such as sandstone, conglomerate and tuff formations, most of the area in the east is young gravel fans and also there is clay flat extends from the east to the center. The area in the west consists of young gravel fans, gravel plains and there is a saline flat as well. Most of the area has low-density rangeland or poor pasture land use land cover (see Hamraz et al. 2015).

Mathematical solution of the groundwater model solves the math form of the mass balance equation in one region and generally produces by continuum an approach to the surroundings. Every component of the mass balance equation indicates a specific value of parameter in unit of surface, volume or time. Briefly, a mathematical model that is applied for the groundwater simulation is a combination of numerical values of different parameters in the balance equation. In other words, the balance equation is written for a limited area of the aquifer while it is generalized for the entire zone. Accordingly, the result of groundwater modeling can be affected by different errors, e.g. the errors arise from the groundwater conceptual model, the approximate solution of groundwater parameters and unknown interaction mechanism between various components and properties.

In practice, the establishment of a groundwater numerical model requires different parameter sets such as the hydraulic conductivity (k), transmissivity, storage coefficient (or specific yield (Sy)) and dispersity. More importantly, because of hydrogeology data scarcity (especially in an arid region), the variability of parameters in space and time, incorrect setting of aquifer properties (location, type, number of layers, distribution, etc.), and the scaling effect of parameters, groundwater modeling can yield different errors and misfit. A conceptual model is the basis of the numerical model, and actual hydrogeological conditions are often simplified incorrectly by groundwater conceptual model (Rojas et al. 2008). In particular, physical consistency in the mass and energy fluxes of groundwater process, which is a continuity mechanism, is usually formulated to make predictions at discrete moments of time. Together, all those variabilities and complexities increase the bias and error in simulation and challenge groundwater modeling particularly in arid regions.

This study used MODFLOW to simulate steady and unsteady conditions in the Birjand aquifer under different drawdown and management scenarios. MODFLOW is a three-dimensional finite-difference groundwater model which employs the above relationships to compute the level of water-table depth in different points of the aquifer. The model was released in 1984, upgraded and modified to perform simulations in both steady and unsteady conditions. In this research, a distributed hydrogeology model was established using all quantitative information that requires simulating a functional relation between prediction and observation.

In this study, the procedure used to develop the groundwater conceptual model was proposed by Izady et al. (2014). Based on their method, the groundwater model consists of six steps: (1) collection of all available data and information; (2) verification and setup of controlling observations; (3) defining the aquifer geometry; (4) primary estimation of aquifer hydrodynamic properties; (5) identification of aquifer recharge and discharge values; (6) integration of the results from other steps to deliver the overall conceptual model.

In this study, seven data coverages including aquifer boundary conditions, piezometers, pumping wells, surface recharges, drainage information, hydraulic conductivity and specific yield were incorporated in the MODFLOW model to establish groundwater numerical model. Boundary conditions in MODFLOW are established using constant head boundary, head dependent flux (river, drain and general-head boundary packages), and known flux (recharge, evapotranspiration, wells, and stream). This boundary encompasses nine inputs and one output front as a drain. Birjand aquifer acts as a transit aquifer due to its differential hydraulics gradients through various zones. Moreover, this aquifer recharges/discharges by surrounding aquifers at inflow/outflow boundaries which cause conditional and temporary effects. In this study, recharge and discharge locations at inflow and outflow boundaries were identified using a group of cells or grids called water front entrance. For all those boundaries, a specific-head-boundary condition was considered to be constant at each cell number in the model.

Since the study area is categorized as an arid region with small amount of precipitation, only drinking and agricultural backwaters were used as surface recharges for related fronts in the model. In order to construct the aquifer model, MODFLOW was first simulated using steady state and the results were used to define homogeneous zones. Therefore, the aquifer was divided into 17 homogeneous zones to derive the hydraulic conductivity and specific yield values. Next, several setups, including model network design, choice of stress period and time step, determination of initial state, boundary condition, type and number of saturated zones were also defined in order to construct and setup groundwater model. Due to geological heterogeneity, the aquifer was divided into 1,077 cells or grids to solve partial difference equations. The gridding of the plain was designed into 34 rows and 94 columns so that each layer cell is squared to 500 m × 500 m (see Figure 4).

Hereafter, the cells that are positioned outside of the aquifer boundary were assigned to zero code, meaning that they have no influences on the modeling procedure. Though, properties of model cells can be specified individually, grouping the cells with similar properties into homogenous zones significantly improved and eased modeling procedure.

In this study, modeling of the groundwater flow in the Birjand plain was performed for a one year period starting in October 2010 to September 2011 (a hydrologic water year) as the calibration period, whilst October 2011 to September 2012 was used for the validation period. For modeling purposes, 12 month stress periods with a 10 day time step were identified and used for the modeling approach (i.e. three time steps were considered for each stress period). For completing the numerical model, absolute values of bed rock depths, topography and initial water level measurements were interpolated using Kriging method then assigned for the entire cells in the network.

where, Var_hi and Var_low denote the maximum and minimum amounts of variables, respectively, and ‘α’ is an integer digit that controls maximum ELR.

• (4) Letting Cuckoos lay eggs inside their corresponding ELR.

• (5) Eliminating the eggs with lower value of objective function.

• (6) Defining the value of objective function for each adult Cuckoo.

• (7) Limiting the maximum number of Cuckoos in the network.

• (8) Categorizing Cuckoos and define the best habitant.

In this study, K-means clustering method (MacQueen 1967) is used for categorizing the COA model. A group of 3–5 K is considered and the mean profit of each group is determined according to the variation range of the Cuckoo profit function. So each Cuckoo belongs to a group that has the least distance to the mean of the corresponding group.

• (9) Immigration of new Cuckoo towards the best habitant.

Parameters that are used in this study for the COA are determined by automatic trial-error approaches (see Rajabioun 2011) that are listed in Table 1.

AMALGAM is initiated using a random initial population size of s, generated by the Latin hypercube sampling method. Then, each parent from the initial population is assigned a rank using the fast non-dominated sorting algorithm (Deb et al. 2002) as used in the NSGA-II algorithm. The population of offspring (with size s) is subsequently created using the multimethod concept: instead of utilizing a single biological or physical model for reproduction, a diverse set of operators is run simultaneously. The individual operators are favored adaptively accordingly to their reproductive success during the search. The offspring and parent population is ranked together and members for the next population are chosen from subsequent non-dominated fronts based on their rank and crowding distance (Deb et al. 2002). The new population is then used to generate offspring and the aforementioned algorithmic steps are repeated until convergence is achieved. Note that the adaptive search properties of AMALGAM reduce the need for tuning of the algorithmic parameters. In fact, the main algorithmic parameter of AMALGAM is its population size. Various synthetic multi-objective benchmark studies have demonstrated that the method works well with a population size of about 50. Readers are referred to Vrugt & Robinson (2007) for further description and information on the AMALGAM algorithm.

In this study 5,000 function evaluations with 50 population sizes were used for optimization purposes. Twelve monthly discharge parameters were then used to ensure a reliable and computationally efficient solution to multi-objective optimization problems. The model reached the convergence criteria after about 3,000 simulation runs. The pumping rates along with three minimization objectives, i.e. minimizing shortage affected by the failure to supply, MSI and minimization of extent of drawdown within pre-specified regions, were then chosen to define optimal solution for groundwater drawdown and subsidence.

A combination of simulator optimizer model was illustrated to simulate and quantify aquifer characteristics in the MATLAB program. Simulator model (i.e. MODFLOW) and optimizer algorithms (COA-AMALGAM) should be coupled in order to calibrate the hydrodynamic parameters of the aquifer and to determine pumping rates so the simulator model and the optimizer algorithms were connected by writing the interface program in MATLAB that leads to the simulator-optimizer model. We used MODFLOW simulator model under MATLAB platform (mfLab program) which was developed by Theo Olsthoorn (2013) at Delft University of Technology, the Netherlands. It is interesting to note that for numerical solving purposes the partial difference equation of flow, MF2005NWT (MODFLOW-2005- Newton-Raphson formulation) program, that contains the Newton-Raphson solver to improve the solution of unconfined groundwater-flow problems along with Upstream Weighting Package (UPW), were also used in this study. UPW was used to specify aquifer properties controlling flow movement between cells in the MF 20005-NWT and MF-OWHM (MODFLOW-2005 – One Water Hydrologic Flow Model) approaches.

Decision variables in the optimization model are based on the sum of the discharge amounts from the aquifer (190 wells) during different months (12 variables). The three objective functions used in this study are defined below.

The resulting hydraulic conductivity field is highest at zone 13 whereas it is low at zone 2 (close to outflow). The specific yield field, on the other hand, maximized at zone 7 whilst minimized at zone 6. The resulting RMSE and NRMSE are, respectively, in a range 0.73–0.89 m and 0.02–0.025 obtained by COA optimization algorithm (see Table 2). High hydraulic conductivity at zone 13 may be temporarily confined to the groundwater regime to change aquifer characteristics. Since both parameters showed a high range of variability through the entire aquifer, it appears the hydraulic response of the aquifer to recharge and pumping processes is highly variable. In other words, groundwater regime may not track a steady-state condition. Variability in the hydraulic conductivity may further destabilize transmissivity (the hydraulic conductivity multiplied by the saturated thickness of the aquifer) at the vertical range of the aquifer. In groundwater studies, what is of the greater interest, however, is the significant positive and negative contributions to groundwater modeling made by hydraulic properties or any other parameters.

The specific yield varies in a range of 0.01 to 0.45. Generally, zones with higher specific yield most likely present lower hydraulics conductivity (see Figures 7 and 8) indicating a reverse relationship between these two parameters. Moreover, variability in specific yields may lead to fluctuation of hydraulic response to transient stress in this aquifer. Together, this variability reflects the highly non-linear characteristics for the groundwater regime under study.

In this study, mfLab program was modified slightly to print the volumetric specific yield after simulating each period. The specific yield values were converted to equivalent water thicknesses divided by the study area (sqkm). A number of MATLAB (Mathworks, Inc.) scripts were developed in this study to facilitate algorithmic runs and to generate MODFLOW input files. All runs were carried out on a PC with an Intel Core i5 3210M 2.5 GHz up to 3.1 GHz CPU. The average running time took about 42 hours, depending on the algorithm configuration.

The error indices for simulated and observed head values were demonstrated and presented in Table 2. The error estimation suggests well calibration results for the model. At the validation period, the estimated error is also low, indicating a successful calibration. The results of the model calibration and validation ultimately confirm the validity of the MODFLOW simulation and the efficiency of the COA to compute groundwater properties. It should be noted that any type of uncertainty, such as error in the measured data, the specified initial and boundary conditions, or the overall conceptual model, have a significant effect on modeling. The focus of this study is on optimization procedures (by reducing the error and misfit) not on uncertainty assessment, though, the parameter uncertainty associated with defining the nine pathways and outflow and ground water parameters were addressed for this aquifer by Hamraz et al. (2015).

Based on error indices, the COA significantly reduced the error and misfit by presenting the lowest error during the calibration and validation periods. COA optimization method described herein can provide modelers with a tool that is efficient to implement and can be used in difficult modeling efforts. Groundwater water-table depths in the beginning and ending of the simulation period were also computed and the results for the entire aquifer are presented in Figures 9 and 10.

Comparing groundwater water-table depths indicate an increasing trend in the head values. This phenomenon may be explained by the influences of local pumping activities (drinking or agricultural wells) at the beginning of the simulation period (begins in the October) and replenishment of flow afterward (mostly by neighboring aquifers). These intensive pumping activities might have smooth influences on the shape of the simulated head curve at the middle of the aquifer where the second and third wells are located. In saying that, the simulated groundwater water-table depth represents the values at the center of the cell or grid and occasionally the boreholes are located at those places.

According to Figure 10, groundwater-table depth declined from east to west and southwest parts of the aquifer so the general direction of the groundwater flow is in the same direction. By comparing Figures 9 and 10, we realized that by continuing the current withdrawals, groundwater level drops down most likely in the west zone. This might be related to over-drawdown, lack of groundwater recharge and natural drainage activities as well.

Optimized parameters obtained using COA optimization algorithm were applied in AMALGAM (5,000 evaluations) to optimize predefined objective functions. The main goal was to visualize the tradeoffs between sustainability of groundwater use and drawdown considerations. For conciseness and clarity, our results include aggregated results as opposed to individual efficient pumping or recharge rates for each well. In order to evaluate the influence of different basic assumptions on the results, three policy-analysis options are tested: optimal solution A (or scenario A) in which maximum pumping rates are assumed with minimizing the sum of shortages while having maximum decline in the aquifer water table; optimal solution B (or scenario B) in which pumping rates are allowed to meet the demand with no decline in the aquifer water-table depth and there should be no variability in the aquifer water table at the beginning and end of a water year (October to September); and optimal solution C (or scenario C) in which pumping rates are assumed to be lowest rates (minimizing the drawdown) thus increasing the aquifer water-table depth during the beginning and end of a water year.

According to Equation (4), MSI (caused distribution of shortages in different months) has a direct relationship with shortage. So it is expected that MSI has an inverse relationship with drawdown values (see Figure 12(b)). The sum of optimum discharge amounts (decision variables) in the simulation period is showed for the three samples of optimal solutions from different parts of the Pareto front (solutions A indicated by black dot, solution B indicated by green dot, and solution C indicated by blue dot in Figures 11 and 12). According to Figure 12, solutions A and C (black and blue points), which address zero and 27.24 MCM deficits, respectively, are selected in order to meet the first (minimizing the sum of shortages) and the third (minimizing the drawdown) objective functions. The ideal solution would be a scenario in which all three objectives equally likely contribute to the optimization model. Green point (solution B) represents the optimal solution (ideal scenario) where storage and MSI are almost equal to zero. In practice, when drawdown is close to 0.4 m all demands are met in the aquifer. In solution B, pumping rates at existing wells are allowed to vary over time (however with no variability in the water-table depth at the beginning and ending of a water year), thus creating a more flexible policy that can react to future changes in forcing conditions such as water demands or natural recharges. Each point of the non-dominated sets computed for all solutions represents a specific policy alternative that is to be compared with the others in order to decide management actions. These alternatives have in common the fact that they are different from groundwater management perspectives, that no improvement can be achieved for a particular objective without diminishing the satisfaction associated with the other objectives.

Based on this result decision makers have three optimal solutions for the Birjand aquifer. Specifically, decision makers can select any (drawdown) solution in an optimal range of −40 to +40 cm/year. There are significant differences between water-table depths during maximum (black dot) and minimum (blue dot) drawdowns over the aquifer. If one of the optimal solutions, A, B or C, is used for the groundwater management, it should be ensured that maximum demands of the area are met. Though, the problem of selecting which alternative remains open depending on management policy and may change over time.

It is important to emphasize that optimal solutions B as an intermediate scenario, presented good performances in reducing the objective functions considered in this study. This confirms that for arid case study, intermediate solution was easy for the developed tool to converge. Specifically, optimal solution B was able to fairly satisfy demands, and minimized the drawdown over the aquifer system.

This study first employed a COA algorithm to optimize MODFLOW hydrodynamic parameters then used optimized parameters in the AMALGAM multi-objective algorithm to define Pareto optimal solutions and developed groundwater management scenarios. The proposed multi-objective problem formulation integrated minimizing shortage affected by the failure to supply, MSI and minimization of extent of drawdown within pre-specified regions. Coupling two different algorithms to the MODFLOW model revealed that the developed algorithm was capable of computing a set of optimal solutions displayed on a Pareto front for an arid groundwater system. This study computed optimal solutions for the entire arid aquifer, but examining those solutions at each groundwater zone may provide useful information for arid groundwater management, although it requires a high computational and scripting demand.

Although both COA and AMALGAM algorithms efficiently quantify groundwater characteristics, this proposed methodology has several limitations. For instance, the proposed modeling approach had a long computation time when trying to optimize predefined solutions, which requires running the MODFLOW model at each solution evaluation. For a larger aquifer system than our case study, the computation time may increase dramatically with the number of cells and stress periods. To overcome this limitation, parallel computing can be used by taking advantage of multi-processor machines. Alternatively, artificial neural network (e.g. Giustolisi & Simeone 2006) can be employed to accurately simulate groundwater shortage and reduce the computation time and runs. The concept of fuzzy set theory (Zadeh 1965) can be further used to assess qualitatively each one of the available efficient solutions for the optimized objectives based on additional criteria.

The COA-AMALGAM algorithm developed in this study can be considered as optimal solutions generator and its linkage to the MODFLOW framework makes it possible to simulate optimal groundwater scenarios that can be easily coupled with other environmental models (e.g. hydrology and water quality models). The enrichment of the MODFLOW environment by the optimization capability of COA and AMALGAM formed a looped ‘Modelling – Optimization – Simulation’ procedure useful for decision makers facing groundwater management problems in arid regions.

The authors appreciate those persons and agencies that assisted in accessing research data. Special thanks are owed to Professor Theo Olsthoorn from the Delft University of Technology for his fruitful discussions and technical support on groundwater modeling and processes. Particular acknowledgment is given to Cuckoo developer Dr Ramin Rajabioun from the University of Tehran, for providing Cuckoo MATLAB code. The second author was partially supported by South Carolina Sea Grant Consortium (Grant# 15520-GA11) and the University of South Carolina (Grant# 15520-16-40787). The source code of mfLab-AMALGAM-COA can be obtained from the first author ([email protected]) upon request.