This study presents the first attempt to link the multi-algorithm genetically adaptive search method (AMALGAM) with a groundwater model to define pumping rates within a well distributed set of Pareto solutions. The pumping rates along with three minimization objectives, i.e. minimizing shortage affected by the failure to supply, modified shortage index and minimization of extent of drawdown within prespecified regions, were chosen to define an optimal solution for groundwater drawdown and subsidence. Hydraulic conductivity, specific yield parameters of a modular three-dimensional finite-difference (MODFLOW) groundwater model were first optimized using Cuckoo optimization algorithm (COA) by minimizing the sum of absolute deviation between the observed and simulated water table depths. These parameters were then applied in AMALGAM to optimize the pumping rate variables for an arid groundwater system in Iran. The Pareto parameter sets yielded satisfactory results when maximum and minimum drawdowns of the aquifer were defined in a range of −40 to +40 cm/year. Overall, ‘Modelling – Optimization – Simulation’ procedure was capable to compute a set of optimal solutions displayed on a Pareto front. The proposed optimal solution provides sustainable groundwater management alternatives to decision makers in arid region.

INTRODUCTION

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.

METHODOLOGY

Geographical location of the study area

The Birjand River Basin (BRB) is located in the arid region of South Khorasan province with a drainage area of 3,408 km2. An arid aquifer system namely Birjand aquifer is located in BRB in latitude and longitude of 32 ° 34′ to 33 ° 8′ north and 58 ° 41′ to 59 ° 44′ east, respectively (Figure 1). It underlies an area of approximately 265 km2. According to Domarton continental classification, Birjand plain is classified as arid region with mean annual precipitation and temperature of 170 mm and 17 °C, respectively. Maximum and minimum elevations are, correspondingly, 2,736 and 1,167 m higher than mean sea level. The slope is deep in the east part of the plain whereas it is gentle toward west.
Figure 1

Geographical location and the geology map of the study area. Groundwater depth reduces sharply toward west. The aquifer is mainly composed of young gravel fans and dasht gravel plains.

Figure 1

Geographical location and the geology map of the study area. Groundwater depth reduces sharply toward west. The aquifer is mainly composed of young gravel fans and dasht gravel plains.

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).

Water demands of the study area are divided into four categories: drinking, agricultural, industrial and services; their amounts are listed in Figure 2. It appears agriculture has the highest demand of water during March to July while this demand declines afterward. In addition, drinking demand increases slightly during summer while industrial demand is almost constant during the entire year.
Figure 2

Amounts of drinking, agricultural, industrial and services demands in the study area based on million cubic meters (MCM). Agricultural demand is high during spring and summer while it is low during winter and late fall.

Figure 2

Amounts of drinking, agricultural, industrial and services demands in the study area based on million cubic meters (MCM). Agricultural demand is high during spring and summer while it is low during winter and late fall.

Groundwater model

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.

The governing flow equations in groundwater modeling

The general form of the governing equation in groundwater modeling is: 
formula
1
where, and denote the hydraulic conductivity tensors, h, Ss and R represent pressure head, specific storage and recharge or discharge (positive and negative) components of the aquifer, respectively. In unconfined aquifer, the thickness of saturated layer varies by groundwater-table depth. Some assumptions proposed by Dupuit (1863) are: (a) the flow is horizontal; and (b) the hydraulic gradient is equal to the slope of the free surface.
If there are two-dimensional and transient flows, the equation is then based on Dupuits assumption and the continuity equation is as follows: 
formula
2
where, Sy denotes the specific yield.

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.

A conceptual illustration of aquifer model

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.

To construct the Birjand groundwater model, a one-layer unconfined aquifer with varying thickness from 5 m to 225 m was considered as a conceptual groundwater system. Based on geological analysis, the deepest section of the bedrock is located in the east, ranging from 150 to 225 m below land surface. The groundwater model was developed based on available data, including surface elevation data, well logs, well locations and measurements, geologic map, hydrography and recharge information. Topographic and geologic maps were first employed to define the plain boundary condition. The spatially distributed hydrodynamic properties were then estimated using the single optimization method. In this study, temporal discharge variation was determined by pumping of 190 wells located inside the aquifer boundary. Recharge takes place through nine pathways as inflow paths located in south, east, north and north-west of the aquifer (Figure 3).
Figure 3

Inflow and outflow pathways of the Birjand aquifer (Hamraz et al. 2015). Overall, nine fronts were defined in the study aquifer and all fronts were incorporated to the model as inputs.

Figure 3

Inflow and outflow pathways of the Birjand aquifer (Hamraz et al. 2015). Overall, nine fronts were defined in the study aquifer and all fronts were incorporated to the model as inputs.

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.

From a total of 201 wells, 11 are piezometric wells in which monthly observed water-table depth data was available and used for calibration. The rest of the wells are pumping ones including 139 agricultural, 31 drinking water and 20 industrial pumping wells. In addition, six springs and 11 qanats exist in the aquifer that are mostly used for irrigation purposes (see Figure 4) and we incorporated them as several wells in the MODFLOW model. Qanat is a gently sloping underground channel with a series of vertical access shafts, used to transport water from an aquifer under a hill. This technology was developed by the Persian people sometime in the early 1st millennium BC.
Figure 4

Position of the pumping, piezometeric and drainage wells along with springs and qanats over the Birjand aquifer in the MATLAB programming environment.

Figure 4

Position of the pumping, piezometeric and drainage wells along with springs and qanats over the Birjand aquifer in the MATLAB programming environment.

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.

Single objective optimization algorithm

Cuckoo optimization evolutionary algorithm was developed by Rajabioun (2011). This algorithm is inspired from the life of a bird called a cuckoo. The way of spawning and special growth of the cuckoo was the main idea of this algorithm. Some birds release themselves from the trouble of nest-building and their duties of parenting therefore they resort to some kind of perspicacity to raise their chicks. These chicks are known as brood parasites and the cuckoo is one of the most well known. The cuckoo destroys one of the host bird's eggs and puts her egg between them. In this way she puts the responsibility of egg keeping to the host bird. In practice, the cuckoo does this job by mimicking the color and pattern of the available eggs in every nest so that the new egg is similar to the rest. More similar eggs to the host eggs have more opportunity to growth and survival. In general, two types of cuckoo are considered in this model, i.e. adult Cuckoo and eggs. Pseudo code of the optimization algorithm is summarized as follows:
  • (1) Initializing Cuckoo habitats (initial response) with some random points on the profit function.

  • (2) Dedicating some eggs to each Cuckoo.

  • (3) Defining egg-laying radius (ELR) for each Cuckoo.

The amount of ELR is calculated by the number of eggs and their distance to destination for each cuckoo: 
formula
3

where, Varhi and Varlow 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.

Cuckoo moves toward the aim habitant; it traverses only the λ % of the path by the deviation of φ radians (see Figure 5). These two parameters help cuckoos search much more positions in the whole environment (Rajabioun 2011). λ is a random number between 0 and 1 while φ is a random number between -ω and ω (ω is appropriately within π/6).
  • (10) If stop condition is satisfied ‘stop’, if not back to the second step.

Figure 5

Immigration of a sample Cuckoo toward goal function (Rajabioun 2011). A group of 3–5 K is considered in the model and the mean profit of each group is determined according to the variation range of the Cuckoo profit function.

Figure 5

Immigration of a sample Cuckoo toward goal function (Rajabioun 2011). A group of 3–5 K is considered in the model and the mean profit of each group is determined according to the variation range of the Cuckoo profit function.

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.

Table 1

COA parameters

Parameters Selective value 
Initial number of cuckoo 50 
Minimum egg of each Cuckoo 
Maximum egg of each Cuckoo 
Maximum number of alive Cuckoo 50 
Maximum number of iterations 100 
Number of groups or clusters 
Ω π/6 
The convergence criterion 1 × 10−10 
Parameters Selective value 
Initial number of cuckoo 50 
Minimum egg of each Cuckoo 
Maximum egg of each Cuckoo 
Maximum number of alive Cuckoo 50 
Maximum number of iterations 100 
Number of groups or clusters 
Ω π/6 
The convergence criterion 1 × 10−10 

AMALGAM algorithm

AMALGAM adaptively and simultaneously employs multiple evolutionary multi-objective algorithms to ensure a fast, reliable, and computationally efficient solution to multi-objective optimization problems (Vrugt & Robinson 2007). The algorithm implements a population-based elitism search procedure to find a well distributed set of Pareto solutions within a single optimization run. In other words, the algorithm employs a population of N solutions, whose offspring are created in a genetically adaptive manner by dividing the creation of N offspring among the sub-algorithms proportional to the success achieved by these sub-algorithms during previous generations (Raad et al. 2010). This process is accomplished by means of an offspring partitioning formula, given as 
formula
4
where, is the number of offspring that sub-algorithm j must generate during generation t, and denotes number of offspring that sub-algorithm j contributes to the next generation.

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.

Setting the MODFLOW-COA-AMALGAM framework

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.

Figure 6 illustrates a schematic framework showing the implementation of AMALGAM and COA in the MODFLOW model. In short, aquifer parameters were optimized using COA and used in the AMALGAM model. The master computer runs the algorithmic part of AMALGAM, and generates an offspring population from the parent population using various genetic solutions (e.g. Vrugt & Robinson 2007). This new population is distributed over a predefined number of computational nodes and these nodes execute the simulation model and compute the objective function of the points received. Next, the master computer collects the results, and follows various algorithmic steps to generate the next generation of points. The aforementioned algorithmic processes are repeated until convergence is achieved.
Figure 6

Sequence steps of optimization processes with a linkage among MODFLOW, COA and AMALGAM. We coupled MODFLOW, COA, and AMALGAM using the MATLAB platform.

Figure 6

Sequence steps of optimization processes with a linkage among MODFLOW, COA and AMALGAM. We coupled MODFLOW, COA, and AMALGAM using the MATLAB platform.

Model performance evaluation

Mean error (ME), mean absolute error (MAE), root mean square error (RMSE) and normalized root mean square error (NRMSE; non-dimensional forms of the RMSE) were used as criteria in the calibration period. ME, MAE and RMSE/NRMSE indicate error in the units (or squared units) of the constituent of interest (Moriasi et al. 2007; Etemadi et al. 2014, 2016), which aids in analysis of the results. Error estimation was quantified using the following equations (Equations (5)–(8)). 
formula
5
 
formula
6
 
formula
7
 
formula
8
where, , and denote the observed head, simulated head and mean value for observed head, respectively, n represents the number of piezometers and m is the number of monthly time steps (i ranges from 1 to 11 while t ranges from 1 to 12 months).

Decision variables and objective functions in multi-objective optimization model

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.

Minimizing the sum of shortage was considered as the first objective function. 
formula
9
where, represents the volume of the aquifer discharge in the period of t (MCM), is the amount of demand in the period of t (MCM) and n denotes the total number of time periods (month).
Since this objective function only minimizes the sum of shortages during the period and ignores distribution of shortages over time, a MSI (Equation (10)) was used as the second objective function. This index is important for making economic and societal solutions (see Hsu & Cheng 2002; Chang et al. 2005; Tu et al. 2008). 
formula
10
where, denotes the amount of shortage in the period of t, is the amount of demand in the period of t and n represents the number of the total periods (month). It should be noted that drinking demand is determined for all periods and is taken initially from the discharge values subtracting the shortage amount of other demands (e.g. agriculture, industry and services).
Finally minimizing the drawdown of water-table depth was used as the third objective function. 
formula
11
where, and are the average value for the aquifer water depths (meter) at the beginning and ending of the simulation period, respectively.

RESULTS AND DISCUSSION

Parameter optimization using COA algorithm

COA optimization was performed using 5,000 simulation runs and the model met the convergence criteria after about 2,000 iterations. Two groundwater parameters, i.e. specific yield and hydraulic conductivity, were then optimized for each zone. Optimized values of the hydraulic conductivity (k) and specific yield for the geological units (Sy) are presented in Figures 7 and 8, respectively. Hydraulic conductivity and specific yield vary significantly through the aquifer system indicating a geological heterogeneous and non-homogenous groundwater system.
Figure 7

The variability of hydraulic conductivity through the entire aquifer. Estimated hydraulic conductivity values range from 1 to 150 m/day from west towards the east part of the aquifer. k is very high in zones 4, 13, 12 and 9 while it is lowest in zone 2. Arguably, zones 4, 13, 12 and 9 are mostly saturated compared with other zones (e.g. 2, 11, and 1) with low values for hydraulic conductivity.

Figure 7

The variability of hydraulic conductivity through the entire aquifer. Estimated hydraulic conductivity values range from 1 to 150 m/day from west towards the east part of the aquifer. k is very high in zones 4, 13, 12 and 9 while it is lowest in zone 2. Arguably, zones 4, 13, 12 and 9 are mostly saturated compared with other zones (e.g. 2, 11, and 1) with low values for hydraulic conductivity.

Figure 8

The variability of specific yield through different zones. Maximum Sy was obtained for zones 2, 5, 7 and 10 whereas the lowest value was driven for zones 3 and 6. Specific yield exhibits the same trend as hydraulic conductivity with a range of 0.01 to 0.45.

Figure 8

The variability of specific yield through different zones. Maximum Sy was obtained for zones 2, 5, 7 and 10 whereas the lowest value was driven for zones 3 and 6. Specific yield exhibits the same trend as hydraulic conductivity with a range of 0.01 to 0.45.

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.

Table 2

Performance criteria for calibration and validation periods

Period MAE RMSE ME NRMSE 
Calibration 0.48 0.73 0.17 0.02 
Validation 0.63 0.89 0.32 0.025 
Period MAE RMSE ME NRMSE 
Calibration 0.48 0.73 0.17 0.02 
Validation 0.63 0.89 0.32 0.025 

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.

Calibration and validation results using COA algorithm

According to Figures 7 and 8, aquifer water-table depth was simulated through 17 zones during the calibration period. The model was verified using 12 month stress data from October 2010 to September 2011 (each period contains three 10 day time scales) to predict observed water-table depth. These outcomes allowed us to calculate aquifer properties through model validation. For a particular groundwater parameter, the estimated values in different zones reflect the contribution made by different areas within the model domain to the estimated value of that particular parameter being assessed.
Figure 9

Groundwater level at the beginning of the simulation period in unsteady state where flow characteristics change over time (unsteady state).

Figure 9

Groundwater level at the beginning of the simulation period in unsteady state where flow characteristics change over time (unsteady state).

Figure 10

Groundwater level at the end of simulation period in unsteady state. There is a significant increasing in the groundwater level (compare to the beginning of the simulation) at the middle and the bottom of the aquifer that we assume it is related to recharging activities by neighboring aquifers and drainage activities.

Figure 10

Groundwater level at the end of simulation period in unsteady state. There is a significant increasing in the groundwater level (compare to the beginning of the simulation) at the middle and the bottom of the aquifer that we assume it is related to recharging activities by neighboring aquifers and drainage activities.

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.

Pumping rates optimization

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.

Twelve month discharge values were selected as decision parameters to optimize the AMALGAM algorithm based on predefined objective functions (scenario A, B and C). Specifically, we considered three different optimal solutions located in different parts of the Pareto solution. Hereby, optimal solutions are simulated and communicated with the MODFLOW framework. The Pareto-optimal solution is showed in Figure 11. Arguably, when the shortage amounts decrease, MSI also decreases and the amount of drawdown increases accordingly (negative values of drawdown represents increasing in the aquifer water level). Obviously, if drawdown rises, the shortage amount will decrease so that when drawdown is equal to 0.4 meter (or 40 cm), all demands are met with no shortage.
Figure 11

Non-dominated Pareto front in a three objective function values. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

Figure 11

Non-dominated Pareto front in a three objective function values. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

Next, in order to better illustrate and interpret the optimal solutions, 2D Pareto fronts (two objective functions) were compared mutually (see Figure 12). According to Figure 12(a), drawdown and shortage have an inverse relationship because a declining trend in the aquifer discharge causes reduction in drawdown. The relationship between shortage and drawdown is approximately linear and could result from the same value for aquifer hydrodynamic parameters in various flow depths.
Figure 12

Mutual comparisons of Pareto fronts in a two-dimensional model. Value of 14.4 MCM exhibited deficit with the amount of MSI of 3.95 in order to maintain the algebraic sum of pumping and recharging rates equal to zero: (a) drawdown (m) vs shortage (MCM); (b) drawdown (m) vs MSI; and (c) MSI vs shortage (MCM). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

Figure 12

Mutual comparisons of Pareto fronts in a two-dimensional model. Value of 14.4 MCM exhibited deficit with the amount of MSI of 3.95 in order to maintain the algebraic sum of pumping and recharging rates equal to zero: (a) drawdown (m) vs shortage (MCM); (b) drawdown (m) vs MSI; and (c) MSI vs shortage (MCM). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

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.

The amount of water-table depth fluctuations also showed for the three selected solutions during simulation period in Figure 13. These three optimal solutions are very close during March while they are far in September. March is an ideal time when drawdown and demands are approximately equal or have minor differences. Although, when a steady decline in the aquifer water level occurs during spring and summer due to various pumping activities, the aquifer recovers substantially from October to March. During March water-table depth in the aquifer is approximately 1,336 to 1,336.25 m while during September (maximum discharge time) this value diminished to 1,335.08 m. Analysis of all optimal solutions indicated that the lower and upper limits of groundwater water-table depths are 1,335.1 m and 1,336.1 m during different alternatives (see Figure 14).
Figure 13

Water-table depth fluctuations for each optimal solution in the simulation period. Solutions A and B have the same fluctuations during October to March. X-axis is based on a water year (October to September). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

Figure 13

Water-table depth fluctuations for each optimal solution in the simulation period. Solutions A and B have the same fluctuations during October to March. X-axis is based on a water year (October to September). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

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.

The performances of optimal solutions among the Pareto front indicate that the contributions of the eleven wells could satisfy the demands and the quality requirements. The optimal contributions ensured maximum (+40 cm) and minimum (−40 cm) drawdowns in the eleven wells at the end of the optimization period. It is obvious from Figures 2 and 14 that agriculture demand is the main reason for water level fluctuation in the Birjand aquifer. Therefore, integrated water resources management such as water saving irrigation, rainfed crops, reducing irrigation area can be proposed as effective alternatives for the Birjand ground water management.
Figure 14

Optimum aquifer discharges in the simulation period for the three examples of solutions. Optimal solution A, B, and C include 0, 14.5, and 27.24 MCM deficits, respectively. X-axis is based on a water year (October to September). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

Figure 14

Optimum aquifer discharges in the simulation period for the three examples of solutions. Optimal solution A, B, and C include 0, 14.5, and 27.24 MCM deficits, respectively. X-axis is based on a water year (October to September). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2016.006.

CONCLUSION AND FUTURE WORK

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.

ACKNOWLEDGEMENTS

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 (sadeghitabas@yahoo.com) upon request.

REFERENCES

REFERENCES
Aquado
E.
Remson
I.
1974
Groundwater hydraulics in aquifer management
.
Journal of the Hydraulics Division – ASCE
100
,
103
118
.
Chang
J. F.
Chen
L.
Chang
C. L.
2005
Optimizing the reservoir operating rule curves by generic algorithms
.
Hydrol. Processes
19
,
2277
2289
.
Deb
K.
Pratap
A.
Agarwal
S.
Meyarivan
T. A. M. T.
2002
A fast and elitist multi-objective genetic algorithm: NSGA-II
.
IEEE Transactions on Evolutionary Computation
6
,
182
197
.
Dupuit
J.
1863
Estudes Theoriques et Pratiques sur le Mouvement desEaux
.
Dunod
,
Paris
.
El-Ghandour
Hamdy A.
Elsaid
Ahmed.
2013
Groundwater management using a new coupled model of flow analytical solution and particle swarm optimization
.
International Journal of Water Resources and Environmental Engineering
5
,
1
11
.
Giustolisi
O.
Doglioni
A.
Savic
D. A.
Di Pierro
F.
2008
An evolutionary multi-objective strategy for the effective management of groundwater resources
.
Water Resources Research
44
(
1
),
doi: 10.1029/2006WR005359
.
Gorelick
S. M.
Remson
I.
Cottle
R. W.
1979
Management model of a groundwater system with a transient pollutant source
.
Water Resources Research
15
,
1243
1249
.
Gorelick
S. M.
Voss
C. I.
Gill
P. E.
Murray
W.
Saunders
M. A.
Wright
M. M.
1984
Aquifer reclamation design: the use of contaminant transport simulation combined with nonlinear programming
.
Water Resour. Res.
20
,
415
427
.
Hamraz
B.
Akbarpour
A.
Bilondi
M. P.
SadeghiTabas
S.
2015
On the assessment of ground water parameter uncertainty over an arid aquifer
.
Arabian Journal of Geosciences
8
(
12
),
10759
10773
.
Hsu
N. S.
Cheng
K. W.
2002
Network flow optimization model for basin-scale water supply planning
.
Water Resources Planning and Management
128
,
102
112
.
Izady
A.
Davary
K.
Alizadeh
A.
Ziaei
A. N.
Alipoor
A.
Joodavi
A.
Brusseau
M. L.
2014
A framework toward developing a groundwater conceptual model
.
Arabian Journal of Geosciences
7
,
3611
3631
.
MacQueen
J. B.
1967
Some methods for classification and analysis of multivariate observations
. In:
Proceedings of 5th Berkeley Symposium on Mathematical Statistics and Probability
.
University of California Press
, pp.
281
297
.
Mategaonkar
M.
Eldho
T. I.
2012
Simulation-optimization model for in situ bioremediation of groundwater contamination using mesh-free PCM and PSO
.
Journal of Hazardous, Toxic, and Radioactive Waste
16
,
207
218
.
McKinney
D. C.
Lin
M.
1994
Genetic algorithm solution of groundwater management models
.
Water Resources Research
30
,
doi: 10.1029/94WR00554
.
McPhee
J.
Yeh
W.
2004
Multi-objective optimization for sustainable groundwater management in semiarid regions
.
J. Water Resour. Plann. Manage.
130
(
6
),
490
497
.
Moriasi
D. N.
Arnold
J. G.
Van Liew
M. W.
Bingner
R. L.
Harmel
R. D.
Veith
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Trans. ASABE
50
,
885
900
.
Olsthoorn
T. N.
2013
mfLab: Environment for MODFLOW suite groundwater modeling
.
http://code.google.come/p/mfLab (accessed 19 December 2013)
.
Raad
D.
Sinkse
A.
Vuuren
J.
2010
Multi-objective optimization for water distribution system design using a hyperheuristic
.
Journal of Water Resources Management
136
,
592
596
.
Rajabioun
R.
2011
Cuckoo optimization algorithm
.
Applied Soft Computing
11
,
5508
5518
.
Reed
P.
Minsker
B. S.
Goldberg
D. E.
2001
A multi-objective approach to cost effective long-term groundwater monitoring using an elitist nondominated sorted genetic algorithm with historical data
.
J. Hydroinf.
3
,
71
90
.
Saafan
T. A.
Moharram
S. H.
Gad
M. I.
KhalafAllah
S.
2011
A multi-objective optimization approach to groundwater management using genetic algorithm
.
International Journal of Water Resources and Environmental Engineering
3
,
139
149
.
Siegfried
T.
Bleuler
S.
Laumanns
M.
Zitzler
E.
Kinzelbach
W.
2009
Multi-objective groundwater management using evolutionary algorithms
.
IEEE Transactions on Evolutionary Computation
13
,
229
242
.
Tu
M. Y.
Hsu
N. S.
Tsai
F. T. C.
Yeh
W. W. G.
2008
Optimization of hedging rules for reservoir operations
.
Water Resource Planning and Management
134
,
3
13
.
Vrugt
J. A.
Robinson
B. A.
2007
Improved evolutionary optimization from genetically adaptive multimethod search
.
Proc. Natl. Acad. Sci. U.S.A.
1043
,
708
711
.
Vrugt
J. A.
Gupta
H. V.
Bastidas
L. A.
Bouten
W.
Sorooshian
S.
2003
Effective and efficient algorithm for multi-objective optimization of hydrologic models
.
Water Resour. Res.
39
,
1214
,
doi:10.1029/2002WR001746
.
Wanakule
N.
Mays
L. W.
Lasdon
L. S.
1986
Optimal management of large-scale aquifers: methodology and applications
.
Water Resources Research
22
,
447
465
.
Yapo
P. O.
Gupta
H. V.
Sorooshian
S.
1998
Multi-objective global optimization for hydrologic models
.
J. Hydrol.
204
,
83
97
.
Zadeh
L. A.
1965
Fuzzy sets
.
Inf. Control.
8
,
338
353
.
Zhu
X.
Wu
J.
Wu
J.
2006
Application of SCE-UA to optimize the management model of groundwater resources in deep aquifers of the Yangtze delta
. In:
Proceedings of the First International Multisymposiums on Computer and Computational Sciences (IMSCCS'06) June 20–24 Hangzhou China,
vol.
2
, pp.
303
308
.