Abstract

This study addresses the issue of optimal management of aquifers using a mathematical simulation- optimization model which relies on the stability of water quality and quantity, considering salinity. In this research first we developed a hydrological model (SWAT) to estimate recharge rates and its spatiotemporal distribution. Then, groundwater simulation of the basin was simulated and calibrated using MODFLOW 2000 and water quality was simulated and calibrated using MT3DMS. Afterwards, a multi-objective optimization model (MOPSO) and embed simulation models as tools to assess the objective function was carried out in order to produce a simulation-optimization model. Finally, a sustainability index to assess Pareto front's answers and three management scenarios (continuing previous operation, 30% increasing and reduction in previous operation) was developed. The results show that the majority of Pareto optimal answers have more sustainability index than a 30% reduction of operation with the best answer of 0.059. Relatively, the sustainability index of 30% reduction of operation is 0.05.

INTRODUCTION

Groundwater is the most significant source of water for a variety of uses, including industrial, irrigation, domestic and drinking purposes and provides approximately one-third of the world's freshwater consumption (Moreaux & Reynaud 2006). Due to increased demand and mechanized pumping technologies, excessive withdrawal of water is taking place from the aquifers. As a result lowering of the groundwater level has been reported in many regions of the world (e.g. Middle East and North Africa, India, China, Japan and Spain) that have faced critical water-resource sustainability issues (Llamas & Custodio 2002). Again, there may be groundwater contamination due to intensive use of fertilizers, herbicides, pesticides, disposal of human and industrial waste. Due to excessive withdrawal of water from the aquifer, more and more contaminants reach the well point due to an increased flow of groundwater towards the well. This leads to an increase in concentration of contaminants at the well locations. Therefore, efficient management strategies should be taken for optimal withdrawal of water and concentration of contaminants at the well location that is below a certain permissible limit. One of the best methods of determining the sustainable management strategy for a groundwater system may be the combined use of simulation-optimization models. While simulation models basically provide solutions of the governing equations of groundwater flow, optimization models identify an optimal management and planning strategy from a set of feasible alternative strategies. After obtaining the feasible management strategy, system behavior is evaluated through the simulation model (Das & Datta 2001).

The groundwater resources management problems have been solved by using simulation (Hunt et al. 1998; Ala Eldin et al. 2000; Pint et al. 2003; Pisinaras et al. 2007; 2013; Rojas & Dassargues 2007; Wen et al. 2007; El Yaouti et al. 2008; Michael & Voss 2008; Budge & Sharp 2009; Wondzell et al. 2009; Zghibi et al. 2011; Singh & Panda 2012; Cao et al. 2013; Narula & Gosain 2013; Yihdego & Becht 2013; El Alfy 2014; Izady et al. 2015; Pulido-Velazquez et al. 2015; Eltarabily et al. 2016) and simulation-optimization (Gorelick et al. 1984; Wang & Zheng 1998; Başaǧaoǧlu & Mariňo 1999; Gordon et al. 2000; Zheng & Wang 2002; McPhee & Yeh 2004; Ayvaz & Karahan 2008; Kourakos & Mantoglou 2008; Singh & Minsker 2008; Ayvaz 2009; Yang et al. 2009; Zou et al. 2009; Gaur et al. 2011; El-Ghandour & Elbeltagi 2013; El-Ghandour & Elsaid 2013; Salcedo-Sánchez et al. 2013; Zekri et al. 2015; Heydari et al. 2016) models. However, previously, researchers have actively sought to combine simulation models with optimization techniques to address the groundwater management problems. Matott et al. (2006) carried out a comparative study by coupling different optimization techniques (genetic algorithm (GA), simulated annealing (SA), conjugate gradient (CG) and particle swarm optimization (PSO)) with a simulation model. From the study they concluded that PSO attains better results with less computational time, compared to the other optimization techniques. Also, many researchers applied PSO for groundwater management problems. They found that PSO converged to the global optimal solution with less computational time and yield better results. Sandoval-Solis et al. (2010) developed Loucks’ (1997) SI in combination with performance criteria proposed with Hashimoto et al. (1982), McMahon et al. (2006) and Moy et al. (1986).

In the present study, the transient groundwater flow and transport processes were first simulated using the MODFLOW and MT3DMS package of the groundwater modeling system (GMS) software by using an integrated modeling framework that consists of a watershed agriculturally based hydrological model (Soil and Water Assessment Tool, SWAT) by a 3-D groundwater flow model developed in MODFLOW, and using a total dissolved solids (TDS) mass-transport model in MT3DMS. The SWAT model outputs (mainly groundwater recharge, considering irrigation needs under evapotranspiration (ET) and precipitation) are used as MODFLOW inputs to simulate changes in groundwater flow and storage. The groundwater velocity field, which was obtained by MODFLOW simulation, are used as MT3DMS inputs for assessing the fate and transport of TDS. For calibration of the quantitative model, a manual method along with the automated method by PEST was used but for calibration of the qualitative model, only the manual method was used. Thereafter a MATLAB function was developed to execute MODFLOW and MT3DMS in a MATLAB environment, by reading the input files generated by GMS and using the executable file of MODFLOW and MT3DMS. Thereafter, the MOPSO based optimization algorithm was coupled with the MODFLOW- and MT3DMS-based simulation models. The coupled MODFLOW-MT3DMS-MOPSO model was used to solve the groundwater management problem of Esfahan-Borkhar unconfined aquifer. The objectives of the groundwater management model are to minimize quantity and quality fluctuation. Then, a sustainability index was developed to select the best point of Pareto front. The sustainability index was calculated using the performance criteria of reliability, resilience and vulnerability for each cell of the simulation models by MATLAB. Innovation of this study is the use of a new developed sustainability index to choose the best answer of Pareto front. The SI developed with Sandoval-Solis et al. (2010) was quantitative and belongs to water users. The SI developed in this study is quantitative and qualitative and belongs to water resources.

METHODS

Study area description

Esfahan-Borkhar (Figure 1) is one of the study areas of Iran and is located in Gavkhouni Basin (the terminal basin of the Zayandeh River). The area covers 3,385 km2 (2,688 km2 plain and 1,500 km2 aquifer) and is one of the largest plains in Isfahan Province. The highest altitude of the area is about 2,535 m above mean sea level and the lowest is about 1,540 m above mean sea level.

Figure 1

Location of the Esfahan-Borkhar study area.

Figure 1

Location of the Esfahan-Borkhar study area.

Modeling framework

Traditionally, simulation models were used to obtain the answer of ‘what if’ and optimization models answer the question of ‘what is the best’ under a given set of conditions (Singh 2014). Therefore, the optimal groundwater management alternatives may not be achieved using either simulation or optimization techniques alone. Hence, the combined use of simulation and optimization models is vital (Singh & Panda 2013). During the last three decades, a simulation-optimization modeling approach has been used to solve the real world problems of water resources. Modeling in this research consists of three parts: simulation, simulation-optimization models and development of sustainability index for choosing the best point of Pareto front.

Data for model applications were prepared from physical parameters and hydrogeological properties of the aquifer. Physical parameters included aquifer thickness, hydraulic conductivity, pumping and recharge rates. Hydrogeological properties included geologic formations and a topographic map with wells and locations of boreholes. Then, the study area model domain was identified. The model boundary and boundary conditions were determined using the hydrological and geological data. Table 1 shows the pumping wells tested.

Table 1

Pumping wells tested

No. Wells coordinate (UTM)
 
Transmissibility (m2/day) Depth (m) 
551,380 3,614,183 197 162 
568,500 3,613,899 616 95 
571,823 3,632,840 538 170 
560,100 3,635,380 300 102 
576,185 3,620,462 120 190 
561,546 3,635,298 2,200 135 
555,048 3,636,002 1,415 140 
556,749 3,630,563 650 180 
553,750 3,626,750 350 140 
10 563,469 362,022 300 95 
11 550,600 3,632,500 2,200 106 
12 562,900 3,649,700 1,400 100 
No. Wells coordinate (UTM)
 
Transmissibility (m2/day) Depth (m) 
551,380 3,614,183 197 162 
568,500 3,613,899 616 95 
571,823 3,632,840 538 170 
560,100 3,635,380 300 102 
576,185 3,620,462 120 190 
561,546 3,635,298 2,200 135 
555,048 3,636,002 1,415 140 
556,749 3,630,563 650 180 
553,750 3,626,750 350 140 
10 563,469 362,022 300 95 
11 550,600 3,632,500 2,200 106 
12 562,900 3,649,700 1,400 100 

To simulate groundwater, a link of three models was established. This required the practical integration of operational models that not only represent all of the relevant processes in the hydrologic system in a physically meaningful way, but also are simple enough to allow large-scale basin-wide applications. The characterization of the land phase of the hydrological cycle is essential for assessing the impacts of climate and land use on the temporal and spatial distribution of groundwater recharge. The tool was selected for this purpose was SWAT. While the quasi-distributed SWAT model is capable of properly simulating the spatiotemporal distribution of groundwater recharge rates (at the spatial resolution given by their hydrologic response units, HRUs), its groundwater module is lumped. Therefore, distributed groundwater parameters (such as hydraulic conductivities and storage coefficients) cannot be represented, and the approach is very limited for expressing the spatial distribution of groundwater levels and groundwater flow dynamics. The proposed modeling framework couples the SWAT watershed model with the fully distributed groundwater model MODFLOW (McDonald & Harbaugh 1988), and finally the multi-species transport model MT3DMS (Zheng & Wang 1999; Pulido-Velazquez et al. 2015) for simulating the TDS in the aquifer system. In this approach, SWAT model outputs are used as MODFLOW inputs and MODFLOW outputs are used as MT3DMS inputs.

For optimizing abstraction of groundwater, according to previous studies the MOPSO algorithm usually converges to the global optimal solution with less computational time and better results (Matott et al. 2006). In order to link simulation and optimization model, a link of three models with embedding method was established. After calibration and validation of MODFLOW and MT3DMS in GMS, input files for MODFLOW and MT3DMS based codes were generated using GMS. Then MODFLOW-MT3DMS-MOPSO model were developed. Finally, to choose one point of the Pareto front, a sustainability index was developed. Figure 2 shows the modeling framework adopted.

Figure 2

Modeling framework adopted.

Figure 2

Modeling framework adopted.

Watershed modeling

The SWAT is a basin-scale hydrologic model developed by the United States Department of Agriculture, Agricultural Research Service. SWAT is a conceptual basin-scale continuous-time quasi-distributed watershed simulation model for predicting the impacts of management on hydrology, sediments, agricultural production and chemical yields (Gassman et al. 2007). The SWAT interface tool in ArcGIS (Arc SWAT) is used to develop the model input data sets. The SWAT model requires a wide range of data depending on the model processes. It firstly divides the study area into different sub-basins regarding the river network. The sub-basins are further discretized into HRUs consisting of homogenous soil, slope and land use combinations. At the HRU scale, SWAT simulates the processes specified by the user, being able to perform calculations related to hydrological processes, sediment transport, soil temperature, crop growth and pesticide management (Arnold et al. 1998). Water balance is the driving force behind all the processes in SWAT since it impacts plant growth and the movement of sediments, pesticides and pathogens (Arnold et al. 2012). SWAT calculates the water balance equation at each HRU, which includes daily precipitation, run off, ET, infiltration, and return flow components.

In order to develop the SWAT model simulation for Esfahan-Borkhar the following inputs were used, divided into four categories: the digital elevation model, climate data, soil data and land use data. The HRUs were obtained as combination of land use categories, soil categories and slope category. In a preliminary SWAT HRU definition, an excessive number of HRUs was found. In order to reduce it, a filter was applied by not considering, within one specific sub-basin, land use and soil types whose surface area is below 20% of the total sub-basin area excluding irrigated land. After these operations, a final number of 55 HRUs were obtained. Then, the crop management practices (seed time, irrigation and harvest timetable; the water source, etc.) were introduced in the SWAT model using the Arc SWAT interface.

Groundwater modeling system

The GMS is a comprehensive graphical user environment for performing groundwater simulations. The entire GMS system consists of a graphical user interface (the GMS program) and a number of analysis codes (MODFLOW, MT3DMS, MODPATH, SEEP2D, FEMWATER). The GMS interface was developed by the Engineering Computer Graphics Laboratory of Brigham Young University in partnership with the US Army Engineer Waterways Experiment Station. GMS was designed as a comprehensive modeling environment. Several types of models are supported and facilities are provided to share information between different models and data types. Tools are provided for site characterization, model conceptualization, mesh and grid generation, geostatistics and post-processing. GMS includes a new suite of tools in the Map Module for creating high level representations of groundwater models using geographic information system (GIS) objects. These models can be imported/exported between GMS and Arc/Info or ArcView.

Groundwater flow model

A general form of governing equation used in MODFLOW that describes the three-dimensional movement of groundwater flow of constant density through porous media is shown in Equation (1):  
formula
(1)
where Kxx, Kyy, Kzz are values of hydraulic conductivity along the x, y and z coordinate axes, which are assumed to be parallel to the major axes of hydraulic conductivity (L = T), W is the volumetric flux per unit volume and represents sources (W is negative) and/or sinks (W is positive) of water per unit time (T−1), h is the potentiometric head (L), Ss is the specific storage of the porous material (L−1); and t is time (T). Equation (1), when combined with boundary and initial conditions, describes transient three-dimensional groundwater flows in the anisotropic medium and heterogeneous, provided that the principal axes of hydraulic conductivity are aligned with the coordinate directions.

MODFLOW is a fully distributed model that solves the three-dimensional groundwater flow equation using finite difference method (FDM) approximations. The model calculates the hydraulic head at each cell of the FDM grid (for which the aquifer properties are assumed to be uniform), and from there, the flow between cells, stream–aquifer or lake–groundwater interaction, flows through drains. For that purpose, it requires geological and hydrogeological inputs such as soil layer thicknesses between the surface and bedrock, hydraulic parameters at each grid cell (hydraulic conductivity, storage coefficient), boundary and initial conditions and external stresses. The 3-D Esfahan-Borkhar groundwater model was discretized into 136 columns, 117 rows and one layer, with a cell size of 500 × 500 km2. Based on the water table data availability, time steps have been set such that all hydrological stresses are constant. The total number of pumping wells in the area were 3,970, most of them located in the northern part. The initial heads were interpolated from head observation wells using Kriging algorithm in GIS. The hydraulic conductivity was interpolated from 12 exploration pumping wells by the Kriging interpolation method in GIS, although the latter values were further modified during the calibration process. The hydraulic conductivity varies between 0.0001 and 100 m/day, the horizontal anisotropy, between 10−10 and 100, and the storage coefficient, between 0.006 and 0.6. Recharge values were obtained from the SWAT model outputs, although the latter values were further modified during the calibration process between 10−10 and 0.0008 m/day.

Groundwater TDS transport model

The partial differential equation for three-dimensional transport of contaminants in groundwater flow system is shown in Equation (2):  
formula
(2)
where C is the dissolved concentration of species k (ML−3), θ is porosity of the subsurface porous medium, dimensionless, t is time (T), Dij is hydrodynamic dispersion coefficient tensor (L2 T−1), xi is distance along the respective Cartesian coordinate axis (L), vi is seepage or linear pore water velocity (LT−1), which is related to the specific discharge or Darcy flux through the relationship (vi/qi), qs is volumetric flow rate per unit volume of aquifer representing fluid sources (positive) and sinks (negative) (T−1), Rk is the chemical reaction term (ML−3 T−1), Cs is the concentration of the source or sink flux for species (ML−3).

MT3DMS (Zheng & Wang 1999) is a 3-D groundwater solute transport model that solves the groundwater transport equation using a finite-difference approximation, discretizing the spatial domain into cells with equal characteristics and solute concentrations. TDS transport parameters for the Esfahan-Borkhar aquifer were estimated using the values reported in the literature. The model simulates advection and diffusion in TDS transport in the saturated zone. For the model, initial concentration values were interpolated from concentration observation wells data. The conceptual model of the Esfahan-Borkhar aquifer was developed based on the available topographic map of the study area, and geological, hydrogeological, and physical and chemical parameters associated with the aquifer and salinity data. In this study, three types of boundaries can be considered, specified head boundaries, general head boundaries and no flow boundaries.

Multi objective particle swarm optimization

Particle swarm optimization, introduced by Eberhart & Kennedy (1995) is a population-based algorithm in which optimal solutions are searched through a combination of individual learning and social behavior. In PSO, most of the movement of the particles toward the optimum is governed by equations similar to the following:  
formula
(3)
where w is inertia coefficient; c1 and c2 are constants; r1 and r2 are uniform random numbers in the range [0,1]; Pi is the best position vector of particle i in the iteration of nth so far ‘personal’ best; Pg is the best position vector of all particles in the iteration of nth so far global best; xi is the current position vector of particle i in the iteration of (n + 1)th; and vi is the current ‘velocity’ of particle i in the iteration of (n + 1)th (Eberhart & Kennedy 1995). In MOPSO, one of the solutions that is taken from the repository is used instead of Pg since in a multi-objective problem a single solution does not exist (one global best). A number of multi-objective extensions of the particle swarm optimization algorithm have been proposed since the late 1990s. Most of them use Pareto dominance pair-wise comparisons or rankings to drive the search for non-dominated solutions. Most approaches differ basically in two aspects: The way they promote diversity assures that solutions are widespread and cover most of the Pareto front and the way they select the personal and one of the solutions in repository uses to update the position of particles. This paper implements an approach from Coello et al. (2004) where they used an external repository to store non-dominated solutions, and an adaptive grid approach to select one solution from the repository in order to allow this heuristic to handle problems with several objective functions. The multi-objective PSO pseudocode is shown in Figure 3.
Figure 3

MOPSO pseudocode.

Figure 3

MOPSO pseudocode.

Formulation of groundwater optimization model

Decision variable considering in multi-objective optimization problem are pumping flow rates of production wells. The objectives of the groundwater optimization model are to minimize the quantity and quality fluctuation. Specifically, in this aquifer when the water table fluctuates more, the TDS concentration become constant (mainly because of changing of water movement direction and TDS transport influence of abstraction changes) and vice versa:  
formula
(4)
 
formula
(5)
where F1 and F2 are objective functions; Ntp is the number of time steps, which is 125; Nj is the number of cells in model, which is 15,912; Htj is the head of water in the tth time step in the jth cell; H1j is the head of water in the 1th time step in the jth cell; Ctj is the concentration of TDS in the tth time step in the jth cell; C1j is the concentration of TDS in the 1th time step in the jth cell:  
formula
(6)
 
formula
(7)
 
formula
(8)
 
formula
(9)
 
formula
(10)
 
formula
(11)
where is the sum of pumped water from agricultural wells in tp month; tp is the month counter; td is the number of days in tpth month; is the pumping flow rate of kth well in tp month (m3/day); k is the well counter; NW is the number of pumping wells; Ntp is the number of planning months; is the TDS concentration in the kth well in tp month (mg/lit); is the surface water used in tpth month (m3); is the minimum surface water used in tpth month that is zero (m3); is the maximum surface water used in tpth month that is total water demand (m3); is the minimum TDS concentration of used water that is zero; is the maximum TDS concentration of used water that is 2,000 (mg/liter). The demand of water in this region is 600 MCM per year that is supplied by using surface and ground water resources.

After developing the simulation and optimization models mentioned in the previous sections, both models are coupled to solve groundwater management problems. Figure 4 shows the flow chart of the MODFLOW-MT3DMS–MOPSO model. The coupled model is particularly developed to apply the principles of the simulation–optimization approach, where the optimization model repeatedly calls the simulation model to find the optimum solution of the problem. The optimization model calls the simulation model to predict the state variables. The values of those state variables are used to check the constraints and then a penalty will be applied in cases where constraint violations occurred. The whole solution procedure is successively repeated to generate a new solution (discharge of the wells) until the global (or near global) solution is obtained.

Figure 4

Flow chart for coupled MODFLOW-MT3DMS–MOPSO model.

Figure 4

Flow chart for coupled MODFLOW-MT3DMS–MOPSO model.

Performance criteria and sustainability index

Performance criteria are used to evaluate and compare water management policies (Sandoval-Solis et al. 2010). For each time step of calculating the quantity of the sustainability index, deficits are positive when the water is less than the water for the jth point of the source. For each time period t in calculating the quality sustainability index, deficits are positive when the water is more than the water for the jth point of the source. In this study, the water is equal to the water in first time step:  
formula
(12)
 
formula
(13)
Reliability is considered the part of time when the state of source is better than the desirable state of source.  
formula
(14)

Resilience

Resilience is a system's capacity to adapt to changing conditions (World Health Organization (WHO) 2009). According to Hashimoto et al. (1982), resilience is the probability that a system recovers from a period of failure. Moy et al. (1986) used the maximum number of consecutive deficit periods prior to recovery as an alternative definition of resilience. Resilience is the probability that a successful period follows a failure period (the number of times follows ) for all failure periods (the number of times occurred). This statistic assesses the recovery of the source once it has failed:  
formula
(15)

Vulnerability

Vulnerability is the likely value of deficits, if they occur (Hashimoto et al. 1982). Essentially, vulnerability expresses the severity of failures. Vulnerability can be expressed as (1) the average failure (Loucks et al. 2005), (2) the average of maximum shortfalls over all continuous failure periods (Hashimoto et al. 1982; McMahon et al. 2006), and (3) the probability of exceeding a certain deficit threshold (Mendoza et al. 1997). Equation (16) has been applied in this study:  
formula
(16)
Sandoval-Solis et al. (2010) proposed a variation of Loucks' SI, with the index defined as a geometric average of M performance criteria for the ith water users:  
formula
(17)
In order to compare groups of water sources, sustainability by group (SI) was defined:  
formula
(18)

RESULTS AND DISCUSSION

MODFLOW models' calibration and validation

Figure 5 shows the position of the observation wells in the study area. More bolded points show piezometers position and less bolded points show TDS observation wells. Figure 6 shows boundaries and land use of the study area.

Figure 5

Observation wells position.

Figure 5

Observation wells position.

Figure 6

Boundaries and land use.

Figure 6

Boundaries and land use.

The Esfahan-Borkhar groundwater model was calibrated in two steps: (1) steady state calibration and (2) calibration to transient conditions. MODFLOW calibration was carried in monthly time steps for the 7 October 2002–7 November 2002 period in steady state conditions and in 95 monthly time step for the 7 November 2002–7 October 2010 period in transient conditions using 30 piezometers. MODFLOW validation was carried out for the 7 October 2010–7 March 2013 period (29 monthly time steps).

The model was initially calibrated in steady state mode. The steady state model was calibrated by PEST for hydraulic conductivity and horizontal anisotropy values until the observed and simulated value matched at 30 observation wells distributed over the model domain. The water level in the observation wells was calibrated within a difference of 0.5 m between observed and simulated values. Figure 7 shows the calibrated steady state result for 7 November 2002. Green color shows errors less than 0.5 unit, yellow color shows errors between 0.5 and 1 unit and red color shows errors more than 1 unit. The results of the monitoring well indicates that most of the computed groundwater elevation is within a 0.5 m interval from the observed value. Figure 8 shows the validated result for end of the simulation model, 7 March 2013.

Figure 7

Head of water at end of 1st time step (7 November 2002).

Figure 7

Head of water at end of 1st time step (7 November 2002).

Figure 8

Head of water at 125th time step (7 March 2013).

Figure 8

Head of water at 125th time step (7 March 2013).

The steady state result was used in the transient model simulation as the initial condition. The calibration target in the transient state mode was the water table in 30 observations wells. In the transient simulation, the resultant hydraulic conductivity and horizontal anisotropy from the steady state and recharge from SWAT was modified but the storage coefficient (specific yield) was calibrated. Table 2 depicts the calibration and validation of results obtained in all piezometers located in Esfahan-Borkhar aquifer. The errors of water heads shows that the MODFLOW model was adequately calibrated.

Table 2

Quantitative model errors

Errors Mean error (m) Mean absolute error (m) Root mean squared error (m) 
Transient calibration −0.059 0.315 0.510 
Validation −0.012 0.121 0.337 
Errors Mean error (m) Mean absolute error (m) Root mean squared error (m) 
Transient calibration −0.059 0.315 0.510 
Validation −0.012 0.121 0.337 

MT3DMS models' calibration and validation

MT3DMS calibration was carried out for the 7 October 2002–7 October 2010 period (96 monthly time steps) in transient condition using 15 concentration observation wells. MT3DMS validation was performed for the 7 October 2010–7 March 2013 period (29 monthly time steps).

The parameter adjustment for the MT3DMS model was hindered by the lack of data, especially by the lack of a continuous time series with at least one record per month. So, we fitted regression to the existing data to extract monthly data. The initial concentration values were interpolated from concentration observation wells. The TDS load entering the aquifer was estimated. In calibration porosity, longitudinal dispersivity, ratio of horizontal to longitudinal transverse dispersivity to longitudinal, and ratio of vertical transverse dispersivity to longitudinal dispersivity were obtained. Table 3 shows the calibration and validation of results obtained in all of concentration observation wells. As can be seen, MT3DMS results fit the observed values for the majority of the control points. Figure 9 shows the calibrated errors in each observation well for 7 November 2002 and Figure 10 shows the validated errors in each observation well for 7 March 2013. All of the errors in the observation wells were less than 1 (meter for quantity model and mg/liter for quality model).

Table 3

Errors of qualitative model

Errors Mean error (mg/L) Mean absolute error (mg/L) Root mean squared error (mg/L) 
Transient calibration 0.160 0.957 1.666 
Validation −0.150 0.269 0.919 
Errors Mean error (mg/L) Mean absolute error (mg/L) Root mean squared error (mg/L) 
Transient calibration 0.160 0.957 1.666 
Validation −0.150 0.269 0.919 
Figure 9

TDS concentration of water at end of 1st time step (7 November 2002).

Figure 9

TDS concentration of water at end of 1st time step (7 November 2002).

Figure 10

TDS concentration of water at 125th time step (7 March 2013).

Figure 10

TDS concentration of water at 125th time step (7 March 2013).

Predictive simulation

The transient model simulation was used to predict the response of the aquifer system to anticipate changes in hydrological stresses for the next 125 monthly time steps. Scenario 1 is set to continue the previous trend in groundwater abstraction. Figure 11 shows the prediction of hydraulic heads and Figure 12 shows the prediction of TDS concentration at the 250th time step (6 August 2023). In the second scenario we assume that the groundwater abstractions increased by 30%. Figure 13 shows hydraulic head predictions and Figure 14 shows the TDS concentration of water at the 250th time step (6 August 2023). In the third scenario we assume a 30% reduction in groundwater abstraction. Figures 15 and 16 show hydraulic head and TDS concentration predictions at the 250th time step (6 August 2023), respectively.

Figure 11

Prediction of head of water at 250th time step (6 August 2023) under the first scenario.

Figure 11

Prediction of head of water at 250th time step (6 August 2023) under the first scenario.

Figure 12

Prediction of TDS concentration of water at the 250th time step (6 August 2023) under the first scenario.

Figure 12

Prediction of TDS concentration of water at the 250th time step (6 August 2023) under the first scenario.

Figure 13

Prediction of head of water at the 250th time step (6 August 2023) under the second scenario.

Figure 13

Prediction of head of water at the 250th time step (6 August 2023) under the second scenario.

Figure 14

Prediction of TDS concentration of water at the 250th time step (6 August 2023) under the second scenario.

Figure 14

Prediction of TDS concentration of water at the 250th time step (6 August 2023) under the second scenario.

Figure 15

Prediction of head of water at the 250th time step (6 August 2023) under the third scenario.

Figure 15

Prediction of head of water at the 250th time step (6 August 2023) under the third scenario.

Figure 16

Prediction of TDS concentration of water at the 250th time step (6 August 2023) under the third scenario.

Figure 16

Prediction of TDS concentration of water at the 250th time step (6 August 2023) under the third scenario.

MODFLOW-MT3DMS-MOPSO models' result

The developed groundwater management model (MODFLOW-MT3DMS-MOPSO) is applied to Esfahan-Borkhar aquifer system. There were 3,970 pumping wells in the region (black points in Figure 17). If we use all of the pumping wells in the simulation-optimization model, the running time increases too much. Because of the reducing time of model running, the main region was divided into 78 regions and all of the wells in each region were accumulated in the center of the region (blue points in Figure 17). The number of runs for the quantitative and qualitative model is 100,500 (500 particles, 200 iterations), and the stop condition is the end of the number of running the model. The innovative aspects of S/O approach is the use of ‘Mutation Operators’ (Yen & Leong 2010) which were adopted to encourage global and local searches and improve performance of MOPSO's codes in water resources management. The results of the model application includes 30 non-dominated or Pareto-optimal solutions, and is presented in Figure 18.

Figure 17

Situation of wells.

Figure 17

Situation of wells.

Figure 18

Pareto front simulated by MOPSO.

Figure 18

Pareto front simulated by MOPSO.

Calculation of sustainability index

MATLAB software was used to calculate the performance criteria and sustainability index. MATLAB receives the hydraulic head and TDS concentration of each cell respectively as the output files of MODFLOW and MT3DMS programs. Then, reliability, resilience, vulnerability and sustainability index is calculated by MATLAB. Finally, the sustainability index calculated is shown for each cell in GMS. The average of all grid cells is considered as the sustainability index for the entire aquifer. Equations (19) and (20) are used to calculate the sustainability index. (Rel × Res × (1–Vul))Quantity are quantity performance criteria and (Rel × Res × (1–Vul))Quality are quality performance criteria:  
formula
(19)
 
formula
(20)
The results of the calculation of the sustainability index for 30 points of Pareto front and three scenarios are presented in Figure 19. The horizontal and vertical axis show Pareto front points and three developed scenarios and SI, respectively. They show that the majority of Pareto optimal answers have more sustainability index than the third scenario (30% reduction of operation). The first point of Pareto front has maximum sustainability index among 30 points of Pareto front. The sustainability index of the first point of the Pareto front is 0.059, relatively, the sustainability index of 30% reduction of operation is 0.05. Thus, optimal management is more sustainable than reduction of operation.
Figure 19

Sustainability index.

Figure 19

Sustainability index.

CONCLUSIONS

The main motivation of this research is addressing the issue of optimal management of aquifer in the form of a mathematical simulation-optimization programming model and relying on the stability of water quality and quantity. Modeling in this research consists of three parts of simulation, simulation-optimization models and development of the sustainability index for choosing the best point of Pareto front.

Part of the simulation model consists of: (1) development of a hydrological model (SWAT) to estimate recharge rates and its spatiotemporal distribution, (2) development, calibration and validation of a quantitative groundwater simulation model (MODFLOW), (3) development, calibration and validation of a qualitative groundwater simulation model (MT3DMS). Groundwater recharge field, which was obtained by SWAT simulation, was used as MODFLOW inputs and groundwater velocity field, which was obtained by MODFLOW simulation, was used as MT3DMS inputs.

Root mean squared error, mean absolute error and mean error for quantitative model calibration were 0.510, 0.315 and −0.059; for quantitative model validation were 0.337, 0.121 and −0.012; for qualitative model calibration were 1.666, 0.957 and 0.160; for qualitative model validation were 0.919, 0.269 and −0.150, respectively. Errors show that the MODFLOW and MT3DMS models were successfully calibrated and validated. Thus, they can be used in simulation-optimization model and under other scenarios.

Part of the simulation-optimization model consists of development of a multi-objective optimization model (MOPSO) and embedded simulation models as tools to assess the objective function and thus produce a simulation-optimization model.

Part of the development of the sustainability index includes the development of a sustainability index to assess and compare Pareto front's answers and three management scenarios (continuing previous operation, 30% increase and reduction in previous operation). The purpose of developing the SI is the selection of the best Pareto front answer.

The results indicate that the SI in the simulation period is 0.052 and the SI under the first, second and third scenarios are 0.04, 0.033 and 0.05, respectively. While the majority of Pareto optimal answers have a higher sustainability index than the third scenario (30% reduction of operation). The sustainability index with the best answer of Pareto front is 0.059, relatively, and the sustainability index of the third scenario is 0.05.

REFERENCES

REFERENCES
Ala Eldin
,
M. E. H.
,
Sami Ahmed
,
M.
,
Gurunadha Rao
,
V. V. S.
&
Dhar
,
R. L.
2000
Aquifer modelling of the Ganga–Mahawa sub-basin, a part of the Central Ganga Plain, Uttar Pradesh, India
.
Hydrol. Process.
14
(
2
),
297
315
.
Arnold
,
J. G.
,
Srinivasan
,
R.
,
Muttiah
,
R. S.
&
Williams
,
J. R.
1998
Large area hydrologic modeling and assessment part I: Model development 1
.
JAWRA
34
(
1
),
73
89
.
Arnold
,
J. G.
,
Moriasi
,
D. N.
,
Gassman
,
P. W.
,
Abbaspour
,
K. C.
,
White
,
M. J.
,
Srinivasan
,
R.
&
Kannan
,
N.
2012
SWAT: Model use, calibration, and validation
.
Trans. ASABE
55
(
4
),
1491
1508
.
Başaǧaoǧlu
,
H.
&
Mariňo
,
M. A.
1999
Joint management of surface and ground water supplies
.
Ground Water
37
(
2
),
214
222
.
Cao
,
G.
,
Zheng
,
C.
,
Scanlon
,
B. R.
,
Liu
,
J.
&
Li
,
W.
2013
Use of flow modeling to assess sustainability of groundwater resources in the North China Plain
.
Water Resour. Res.
49
(
1
),
159
175
.
Coello
,
C. A. C.
,
Pulido
,
G. T.
&
Lechuga
,
M. S.
2004
Handling multiple objectives with particle swarm optimization
.
IEEE Trans. Evol. Comput.
8
(
3
),
256
279
.
Eberhart
,
R. C.
&
Kennedy
,
J.
1995
A new optimizer using particle swarm theory
. In:
Proceedings of the Sixth International Symposium on Micro Machine and Human Science
,
Vol. 1
,
Nagoya, Japan
, pp.
39
43
.
El Yaouti
,
F.
,
El Mandour
,
A.
,
Khattach
,
D.
&
Kaufmann
,
O.
2008
Modelling groundwater flow and advective contaminant transport in the Bou-Areg unconfined aquifer (NE Morocco)
.
J. Hydro Environ. Res.
2
(
3
),
192
209
.
Eltarabily
,
M. G.
,
Negm
,
A. M.
,
Yoshimura
,
C.
&
Saavedra
,
O. C.
2016
Modeling the impact of nitrate fertilizers on groundwater quality in the southern part of the Nile Delta, Egypt
.
Water Sci. Technol. Water Supply
17
(
2
),
561
570
.
Gassman
,
P. W.
,
Reyes
,
M. R.
,
Green
,
C. H.
&
Arnold
,
J. G.
2007
The soil and water assessment tool: historical development, applications, and future research directions
.
Trans. ASABE
50
(
4
),
1211
1250
.
Gordon
,
E.
,
Shamir
,
U.
&
Bensabat
,
J.
2000
Optimal management of a regional aquifer under salinization conditions
.
Water Resour. Res.
36
(
11
),
3193
3203
.
Gorelick
,
S. M.
,
Voss
,
C. I.
,
Gill
,
P. E.
,
Murray
,
W.
,
Saunders
,
M. A.
&
Wright
,
M. H.
1984
Aquifer reclamation design: the use of contaminant transport simulation combined with nonlinear programing
.
Water Resour. Res.
20
(
4
),
415
427
.
Hashimoto
,
T.
,
Stedinger
,
J. R.
&
Loucks
,
D. P.
1982
Reliability, resiliency, and vulnerability criteria for water resource system performance evaluation
.
Water Resour. Res.
18
(
1
),
14
20
.
Heydari
,
F.
,
Saghafian
,
B.
&
Delavar
,
M.
2016
Coupled quantity-quality simulation-optimization model for conjunctive surface-groundwater use
.
Water Resour. Manage.
30
(
12
),
4381
4397
.
Izady
,
A.
,
Davary
,
K.
,
Alizadeh
,
A.
,
Ziaei
,
A. N.
,
Akhavan
,
S.
,
Alipoor
,
A.
&
Brusseau
,
M. L.
2015
Groundwater conceptualization and modeling using distributed SWAT-based recharge for the semi-arid agricultural Neishaboor plain, Iran
.
Hydrogeol. J.
23
(
1
),
47
68
.
Llamas
,
M. R.
&
Custodio
,
E.
, (eds).
2002
Intensive Use of Groundwater: Challenges and Opportunities
.
CRC Press
,
The Netherlands
.
Loucks
,
D. P.
1997
Quantifying trends in system sustainability
.
Hydrol. Sci. J.
42
(
4
),
513
530
.
Loucks
,
D. P.
,
Van Beek
,
E.
,
Stedinger
,
J. R.
,
Dijkman
,
J. P.
&
Villars
,
M. T.
2005
Water Resources Systems Planning And Management: An Introduction To Methods, Models And Applications
.
UNESCO
,
Paris
.
Matott
,
L. S.
,
Rabideau
,
A. J.
&
Craig
,
J. R.
2006
Pump-and-treat optimization using analytic element method flow models
.
Adv. Water Resour.
29
(
5
),
760
775
.
McDonald
,
M. G.
&
Harbaugh
,
A. W.
1988
A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model
.
Techniques of water-resources investigations of the US Geological Survey
,
USA
.
McMahon
,
T. A.
,
Adeloye
,
A. J.
&
Zhou
,
S. L.
2006
Understanding performance measures of reservoirs
.
J. Hydrol.
324
(
1
),
359
382
.
McPhee
,
J.
&
Yeh
,
W. W. G.
2004
Multiobjective optimization for sustainable groundwater management in semiarid regions
.
J. Water Resour. Plan. Manage.
130
(
6
),
490
497
.
Mendoza
,
V. M.
,
Villanueva
,
E. E.
&
Adem
,
J.
1997
Vulnerability of basins and watersheds in Mexico to global climate change
.
Clim. Res.
9
(
1–2
),
139
145
.
Pint
,
C. D.
,
Hunt
,
R. J.
&
Anderson
,
M. P.
2003
Flowpath delineation and ground water age, Allequash Basin, Wisconsin
.
Ground Water
41
(
7
),
895
902
.
Pisinaras
,
V.
,
Petalas
,
C.
,
Tsihrintzis
,
V. A.
&
Zagana
,
E.
2007
A groundwater flow model for water resources management in the Ismarida plain, North Greece
.
Environ. Model. Assess.
12
(
2
),
75
89
.
Pisinaras
,
V.
,
Petalas
,
C.
,
Tsihrintzis
,
V. A.
&
Karatzas
,
G. P.
2013
Integrated modeling as a decision-aiding tool for groundwater management in a Mediterranean agricultural watershed
.
Hydrol. Process.
27
(
14
),
1973
1987
.
Pulido-Velazquez
,
M.
,
Peña-Haro
,
S.
,
García-Prats
,
A.
,
Mocholi-Almudever
,
A. F.
,
Henriquez-Dole
,
L.
,
Macian-Sorribes
,
H.
&
Lopez-Nicolas
,
A.
2015
Integrated assessment of the impact of climate and land use changes on groundwater quantity and quality in the Mancha Oriental system (Spain)
.
Hydrol. Earth Syst. Sci.
19
(
4
),
1677
1693
.
Salcedo-Sánchez
,
E. R.
,
Esteller
,
M. V.
,
Hoyos
,
S. E. G.
&
Martínez-Morales
,
M.
2013
Groundwater optimization model for sustainable management of the Valley of Puebla aquifer, Mexico
.
Environ. Earth Sci.
70
(
1
),
337
351
.
Sandoval-Solis
,
S.
,
McKinney
,
D. C.
&
Loucks
,
D. P.
2010
Sustainability index for water resources planning and management
.
J. Water Resour. Plan. Manage.
137
(
5
),
381
390
.
Singh
,
A.
&
Minsker
,
B. S.
2008
Uncertainty-based multiobjective optimization of groundwater remediation design
.
Water Resour. Res.
44
(
2
),
404
423
.
Singh
,
A.
&
Panda
,
S. N.
2013
Optimization and simulation modelling for managing the problems of water resources
.
Water Resour. Manage.
27
(
9
),
3421
3431
.
Wang
,
M.
&
Zheng
,
C.
1998
Ground water management optimization using genetic algorithms and simulated annealing: formulation and comparison
.
JAWRA
34
(
3
),
519
530
.
Wen
,
X. H.
,
Wu
,
Y. Q.
,
Lee
,
L. J. E.
,
Su
,
J. P.
&
Wu
,
J.
2007
Groundwater flow modeling in the Zhangye basin, Northwestern China
.
Environ. Geol.
53
(
1
),
77
84
.
World Health Organization
2009
Summary and Policy Implications Vision 2030: The Resilience of Water Supply and Sanitation in the Face of Climate Change
.
WHO Press
,
Switzerland
.
Yen
,
G. G.
&
Leong
,
W. F.
2010
Constraint handling procedure for multiobjective particle swarm optimization
. In:
Evolutionary Computation (CEC), 2010 IEEE Congress
.
IEEE
,
Barcelona
,
Spain
.
Zekri
,
S.
,
Triki
,
C.
,
Al-Maktoumi
,
A.
&
Bazargan-Lari
,
M. R.
2015
An optimization-simulation approach for groundwater abstraction under recharge uncertainty
.
Water Resour. Manage.
29
(
10
),
3681
3695
.
Zghibi
,
A.
,
Zouhri
,
L.
&
Tarhouni
,
J.
2011
Groundwater modelling and marine intrusion in the semi-arid systems (Cap-Bon, Tunisia)
.
Hydrol. Process.
25
(
11
),
1822
1836
.
Zheng
,
C.
&
Wang
,
P. P.
1999
MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User's Guide
.
Alabama University
,
AL
,
USA
.
Zou
,
Y.
,
Huang
,
G. H.
,
He
,
L.
&
Li
,
H.
2009
Multi-stage optimal design for groundwater remediation: a hybrid bi-level programming approach
.
J. Contam. Hydrol.
108
(
1
),
64
76
.