Multi-objective optimization for sustainable groundwater management by developing of coupled quantity-quality simulation-optimization model

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.

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 km² (2,688 km² plain and 1,500 km² 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.

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.

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.

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.

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.

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 km². 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⁻¹⁰ 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⁻¹⁰ and 0.0008 m/day.

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.

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

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.

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.

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

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.

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.

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.