Multi-objective optimal model for sustainable management of groundwater resources in an arid and semiarid area using a coupled optimization-simulation modeling

Excessive exploitation of groundwater resources can increase the concentration of pollutants in addition to the progressive drawdown of groundwater table. In this research, to achieve aquifer quantitative and qualitative (QQ) sustainable development, an optimal scenario for withdrawing from operation wells is proposed. At the first step, the aquifer QQ simulation was carried out with the GMS model. The developed code in MATLAB2018b in the second step provides the link between the simulation and the NSGA-II optimization tools. In the third step, a multi-objective coupled optimization-simulation model based on GMS and NSGA-II developed. Finally, optimal scenario was chosen based on applying the multiple criteria decision-making (MCDM) and Berda Aggregation Method (BAM). The results show that reducing the current withdrawal rate to 51.55% can establish the QQ stability of the aquifer. This decrease in groundwater abstraction has led to a 4.6 m increase in groundwater level (GWL) over 3 years (average 19 cm per month). The spatial and temporal distribution of nitrate concentration after applying the optimal discharge of wells shows the nitrate concentration in central and eastern parts of the aquifer has greatly reduced. Developed sustainable management model can be used to provide a real operation planning of wells to improvement of the QQ status of groundwater in each unconfined aquifer.


Introduction
Groundwater resources are one of the main and most important components of water supply in the miscellaneous sectors. Lack of precipitation, high potential of evapotranspiration, and increasing consumption along with their development during the past decades have caused an increasing pressure on water resources, especially groundwater resources in arid and semiarid regions (Safavi et al. 2010). In Iran, more than 400 plains out of 609 plains due to intensive use have been declared as "forbidden plains." A total of 500 BCM (billion cubic meters) of strategic underground reservoirs have been identified in Iran. Of this amount, 200 BCM is brackish water and only 300 BCM in terms of quality can be used by various consumers. So far, more than 110 BCM have been extracted from this valuable resource. (IWRMC 2017). Excessive withdrawals, drilling of illegal wells, and insufficient monitoring of the number of withdrawals beyond the exploitation permission have led to a significant decrease in well yields, decrease in groundwater level (GWL), increase in pumping costs, decrease in groundwater quality, etc. (Rejani et al. 2009;Zhang et al. 2014;Xiang et al. 2020).
It is not possible to create sustainability in QQ conditions of groundwater resources in many arid and semiarid Communicated by Marcus Schulz. regions, but it is possible to minimize the adverse effects of over-abstraction by using appropriate operation policies. For this purpose, many studies have been proposed and successfully applied to real-world groundwater systems since the 1970s to manage groundwater resources and create sustainable conditions using simulation and optimization tools. In these models, due to the lack of access to important quality parameters such as nitrate and the complexities of aquifer quality modeling, the water quality of the aquifer is usually less considered in the development of operational policies of groundwater. Typically, these models use an aquifer simulation model (which can be an analytical solution or well-known codes such as MODFLOW) and an optimization tool to solve groundwater management problems. By implementing these groundwater management models, it is possible to provide scenarios of aquifer operation (Ebraheem et al. 2003;Karamouz et al. (2005Karamouz et al. ( , 2007; Esteban and Dinar (2013); Farhadi et al. 2016;Banihabib et al. 2019;Nazari and Ahmadi (2019); Sabzzadeh and Shourian (2020); Norouzi Khatiri et al. 2020). In this section, some of these studies are mentioned. Rejani et al. (2009) proposed a non-linear transient hydraulic management model (model 1) and a linear land allocation optimization model (model 2) for the Balasore coastal basin groundwater management. These models have been sued to optimize the pumping rate and determining the optimal cropping patterns maintaining groundwater levels within the desired limits. The results showed the groundwater levels improved significantly under the optimal conditions.
A GA-based simulation-optimization model was presented by Sedki and Ouazar (2011). In developed model, the optimal groundwater exploitation strategies in Rhis-Nekor Plain with the objective of maximizing groundwater withdrawals to supply water demands has been considered. This study shows that the proposed pumping strategy can capture an important amount of the fresh water discharging to the sea for local water supply. Ayvaz (2010) introduced a linked simulation-optimization model to determine the source locations and release histories together with the potential source numbers. In model, MODFLOW and MT3DMS packages were used to simulate the process of flow and transfer of groundwater; then, they integrated the models with the optimization model based on the heuristic harmony search. The results show that the prepared model can be used to solve the inverse pollution source identification problems. Majumder and Eldho (2016) examined the effectiveness of cat swarm optimization (CSO) for groundwater management using the combination of the analytic element method (AEM) and reverse particle tracking (RPT). In this study, three single-objective problems were considered where the objectives are defined as: maximization of the total pumping of groundwater from the aquifer, minimization of the total pumping costs, and minimization of groundwater contamination by capture zone management. The results show the superiority of CSO in comparison with other optimization algorithms. Alizadeh et al. (2017) developed a fuzzy multi-objective compromise methodology based on MODFLOW and MT3D simulation models, NSGA-II, and fuzzy transformation method (FTM) to determine the socio-optimal and sustainable policies for hydro-environmental management of Kavar-Maharloo aquifer. The results indicate the proper performance of the proposed model for determining the most sustainable allocation policy in groundwater management. A linkage between DSSAT, an agronomic model, and MODFLOW on an annual time step in order to assess groundwater conservation strategies in groundwater-irrigated regions was prepared by Xiang et al. (2020).
For QQ modeling of groundwater systems within water resources management models, it is necessary to establish the link between the simulation and optimization models in a desirable way. In previous studies, due to the software limitations of developed groundwater simulation models, the following approaches for using the results of the aquifer simulation model was used in the optimization process: response matrix method (Rejani et al. 2009;Tabari and Soltani 2013;Tabari and Yazdi 2014;Rashid et al. 2014, and etc.), analytic element models (Gaur et al. 2011;Majumder and Eldho 2016, and etc.), meta-models to communicate between simulation and optimization models (Karamouz et al. 2007;Tabari 2015;Alizadeh et al. 2017), and create or modify MODFLOW and MT3DMS codes (Elci and Tamer Ayvaz 2014;Sreekanth et al. 2015;Ayvaz 2016;Ghaseminejad and Shourian 2019;Norouzi Khatiri et al. 2020;etc.).
Review of the previous studies on optimal QQ management of groundwater resources indicated the use of an aquifer-distributed simulation model instead of using simplified relationships to simulate QQ parameters are essential to better describe the behavior of aquifers. Also, adopting management strategies in the abstraction of the aquifer and its control can be very effective in improving the QQ conditions of the aquifer. Accordingly, in this study, with the aim of achieving sustainable development in the QQ operation of the Karaj plain aquifer, the multi-objective management model was developed based on the GMS simulation model and the NSGA-II multi-objective algorithm. In this model, the accuracy of QQ aquifer simulation has direct effects on the results of the optimization model. Therefore, by preparing code in MATLAB environment for direct connection of the multi-objective optimization model with the GMS simulation model by direct access to the input and output files of this simulation model, the aquifer was simulated and calibrated as a cell by cell.
One of the parameters that usually increases in groundwater and has adverse effect on humans is nitrate. This parameter enters the water sources through different ways such as contact of water sources with sewage or discharge of agricultural return flow into the river. Therefore, it is necessary to pay attention to the simulation of this parameter and provide solutions and policies to improve the quality of groundwater in the development of aquifer management models. The concept of optimal use of groundwater resources, considering the development of nitrate pollution, is a relatively new concept that has received less attention. Controlling the concentration of nitrate in the proposed model is another new aspect of this study, which has been considered simultaneously with the quantitative management of the aquifer. Based on this, the innovation of this research can be summarized as follows: • Development of a cellular multi-objective QQ management model of the aquifer. • Preparation of MATLAB code linking GMS and optimization model for extensive aquifer simulation and increasing accuracy in investigation of the QQ behavior of groundwater. • Provide optimal operation policies for each of the operation wells to create sustainable QQ operating conditions.
The results obtained from the implementation of the proposed approach show the high performance of aquifer operation policies in order to create QQ sustainability during a short-term planning period. Therefore, applying this method in other aquifers can reduce the operating costs of groundwater and compensate the water depletion of groundwater and provide the sustainable operation of aquifers.

Study area
Due to the importance of studying critical areas with special political, social, and regional sensitivities, in this study, the aquifer of Karaj plain has been selected as case study (Fig. 1). The area of case study is 255 km 2 and is located in the aquifer of Karaj plain with an average altitude of 1274.14 m above sea level and 48 km northwest of Tehran. The Karaj plain has experienced a large population density in recent year due to its proximity to the capital of IRAN, Tehran. This has led to an increase in abstraction from groundwater resources and as a result a significant decline in GWL and decrease in water quality of the Karaj plain aquifer. The hydrogeological and hydrogeochemical characteristics of the case study are fully described in the paper of Chitsazan et al. (2017).
In order to investigate the QQ behavior of the aquifer of Karaj plain, it is necessary to first prepare an aquifer simulation model and then be calibrated and validated. The time period considered for aquifer modeling is the 2010-2011 water year. The reason for choosing this water year is the completeness of aquifer hydrogeological and hydrogeochemical data. It should be noted that for QQ modeling of the Karaj plain aquifer under steady and unsteady conditions, QQ data related to the 2010-2011 and 2011-2013 water years have been used for calibration and verification, respectively.
In qualitative modeling of Karaj plain aquifer, due to the significant growth and development of urban and agricultural land use in this region, the nitrate parameter has been considered as one of the important and effective parameters in the qualitative degradation of groundwater resources. Spatial and temporal distribution of measured nitrate concentration in drinking wells located in the study area over a 14-year period (2000-2013) are presented in the paper of Chitsazan et al. (2017).

The structure of the proposed approach
Development of groundwater resources operation scenarios for sustainable groundwater abstraction and its quality improvement requires the definition of specific management objectives. For this purpose, a novel approach based on the most appropriate simulation and optimization tools were developed. In this approach, it is first necessary to properly define the objective functions and constraints and introduce them to the developed model. In this study, three objective functions are considered, which are as follows: minimize the sum of drawdown of GWL in drinking wells located in the study area during horizon planning (as first objective function), minimize the sum of nitrate concentration in cells containing operation wells during horizon planning (as a second objective function), and minimize the sum of withdrawal rate from wells during horizon planning (as a third objective function). Due to the nonlinear and complex relationship between groundwater level drawdown, harvesting from aquifer, and nitrate concentration in each of the active aquifer cells, the mentioned goals do not work in the same direction and are considered as conflict objectives.
A remarkable point in the proposed approach is the use of the distributed GMS model to simulate the QQ behavior of groundwater in Karaj plain, which unlike lump models, leads to an increase in the accuracy in calculating the GWL and nitrate concentration parameters.
The mathematical form of the objective functions and constraints are defined as follows: Objective function: where: ΔH tz : The rate of GWL drawdown associated with well z in month t (m).
Qwell tz : Amount of water withdrawn from well z in month t ( m 3 ∕day ) (as decision variable).
QCwell tz : Current status of withdrawal from well z in month t ( m 3 ∕day).
R t : The amount of natural recharge of the aquifer in month t ( m∕day).
Cwell tz : Simulated nitrate concentration in well z and in month t ( mg∕l).
H tz : Simulated GWL in well z and in month t (m). f : A function based on which the quantitative behavior of aquifer (based on GWL parameter) is modeled.
g : A function based on which the qualitative behavior of aquifer (based on nitrate parameter) is modeled.
nwell : Number of wells in the study area. m : Number of months in planning horizon. In this part, introduced equations will be explained. In Eq. (1), which is the first objective function, the control of GWL drawdown (Eq. (5)) in each well is considered. Based on this relationship, it is necessary to first simulate the GWL time series using a validated aquifer model. In order to simulate the GWL, the GMS model is used which is presented in the form of Eq. (6). This model is a graphical user interface for the MODFLOW simulation model. In fact, for each of the solutions provided by the optimization algorithm, it is necessary to run an aquifer simulation model to determine the GWL and its drawdown in each of the simulated cells. Based on Eq. (6), it can be seen that the parameters of recharge, discharge, and GWL of the previous month are needed as input to simulate the quantitative behavior of the aquifer in GMS model.
Satisfaction of the first objective function can be effective in controlling the monthly GWL drawdown in operation wells, reducing pumping costs, and improving the water quality of the aquifer in the long-term with increasing the saturated thickness of the aquifer. In other words, the first objective function plays an effective role in controlling the quantitative stability of the aquifer.
In Eq.
(2), the second objective function of the proposed approach is realized, which is to reduce the nitrate concentration in the operation wells. In order to determine the nitrate concentration in each of the wells, it is necessary to run a calibrated qualitative model, which in this study is MT3DMS and is one of the packages of the GMS model, for different situations of well extraction. This value, which is determined based on Eq. (7), indicates the qualitative behavior of aquifer in the face of stresses due to water abstraction. In Karaj aquifer, the concentration of nitrate has increased to more than the permissible values due to the remarkable withdrawal by wells. Therefore, this groundwater overdraft was controlled by Eq. (4).
The third objective function of this study, which plays an important role in the QQ stability of the aquifer, is to minimize the amount of abstraction from operation wells. Indeed, in this objective function (Eq. (3)), the water supply demands of the region are not considered and long-term operation of the aquifer and attention to the stability of the aquifer in order to improve the water quality of the aquifer are considered as the priority of water withdrawal from wells.
Using these three objective functions, which are also complementary to the sustainable development of the aquifer, can determine the optimal allocation values from each of the wells. Also, based on optimal allocation of the aquifer, other available water sources (such as surface water resources) can be used to supply water deficits.
In order to solve the developed management model, initially, it is necessary to define the problem variables which are known as decision variables or unknowns of the model. In this study, the decision variables considered for each month are the monthly amount of water extracted from existing wells. As there are 2453 active operating wells in the Karaj plain, the total number of decision variables will be 58,872 within 3 years planning horizon.
The reason for considering each of the exploitation wells as a decision variable is the independence in their exploitation, which is mainly managed by the private sector and it is not possible to integrate them regionally. In fact, in the case of aggregation and determination of the optimal amount of water allocation, the optimal amount cannot be properly distributed among the stakeholders.
Also, reducing the number of wells and presenting it in the form of a limited number of wells to apply to the distributed model of the aquifer due to not applying the exact position of the harvest, it can lead to errors in simulating the QQ behavior of the aquifer. It should be noted that this approach can be very effective in developing optimal allocation guidelines from each well for inclusion in the exploitation license and provide appropriate guidance to decision-makers in this area.
In this study, a short-term planning horizon (3 years) was considered to extract the best policies for the operation of wells. Providing optimal aquifer operating policies for a long-term period is not practical for reasons such as changing exploitation approaches as a result of the managerial instability in the organizations in charge of the operation of Iranian aquifers and lack of adequate monitoring of aquifer resources.
To determine the optimal amount of these decision variables, we used the NSGA-II that is one of the most suitable optimization tools. According to this multi-objective optimization algorithm, first, an initial population of random values of water withdrawal from each well (as set of solutions) is generated. Then, the values of the objective functions are calculated for each of the solutions. Using the operators defined in the NSGA-II algorithm, the generated solutions are improved during successive iterations to satisfy the developed objective functions. This process continues until there is no change in the optimal trade-off curve. Under these conditions, it can be stated with great probability that the value provided for the decision variables is near to global optimal and can be used as scenarios for the exploitation of wells in the study area. Based on these optimal allocation amounts from different wells, the optimal operation policies on a monthly or seasonal basis can be formulated. To study the NSGA-II algorithm further, one can refer to Tabari et al. (2020) article. Figure 2 shows the flowchart to achieve an optimal tradeoff between the objective function and how to extract an aquifer operation policy based on the proposed methodology. According to this figure, to achieve the goals of this study, a structure based on the hybrid of GMS simulation Fig. 2 The structure of the proposed multi-objective simulation-optimization model model and the NSGA-II multi-objective optimization algorithm are proposed. Initially, the QQ data sets of the studied aquifer are collected and monitored.
Then, the QQ simulation models of groundwater based on GMS tools were prepared and calibrated under steady and unsteady conditions to determine the status of variation in hydrodynamic coefficients of the aquifer. To ensure the accuracy of the prepared aquifer simulation model, the model was validated based on new data which have not been used in the calibration process.
Since the proper simulation of the QQ behavior of the aquifer in the optimization process is of great importance, therefore, to increase the precision of the aquifer simulation results, a code that can be used to simulate the fully distributed of aquifer within the optimization model was prepared in the MATLAB2018b application environment. This code can call all input and output files of GMS software in MATLAB environment and is able to QQ modeling of aquifer for a short time (approximately 4 s) during a period of 24 months.
In developed MATLAB code, the user will be able to change the status of the stresses in the groundwater system (like recharge and discharge) using a file with the h5 extension in the GMS model and observe variation in GWL and nitrate concentration after the implementation of GMS model. In this code, it is necessary to use a calibrated simulation model of an aquifer, which indicates the actual behavior of the groundwater system in Karaj plain.
By implementing the proposed coupled simulation-optimization model, the optimal trade-off curve between the objective functions is extracted. Each trade-off curve contains numerous optimal scenarios for the operation of wells. Therefore, in this study, the MCDM methods were used to extract the most appropriate exploitation of aquifer scenarios in terms of objectives.
In this study, the MCDM methods used to determine the superior scenario are weighted aggregate sum product assessment (WASPAS), complex proportional assessment (COPRAS), a technique for order preference by similarity to ideal solution (TOPSIS), compromising programming (CP), and modified TOPSIS (M-TOPSIS). Due to the different ranking of scenarios (points located on optimal trade-off curves) in each of the MCDM methods, the BAM method was used to aggregate the results of these five MCDM methods of ranking the solutions generated by the NSGA-II algorithm and determine the final ranking of each scenario. Details of the decision-making methods used can be found in the Banihabib et al. (2017) paper.
Based on the best scenario, the QQ behavior of aquifer is analyzed under optimal operation. In fact, according to the optimal values obtained from this scenario, deciding on existing operating conditions becomes easier and operating and decision-making managers can appropriately present short-term and long-term plans for sustainable operation of the groundwater system.

Aquifer simulation using prepared MATLAB code
In this study, due to the high application of the GMS model and its very appropriate accuracy in investigating the QQ behavior of the aquifer, code was prepared in MATLAB environment to establish a link between this software and multi-objective optimization algorithms. It should be noted that there are other methods such as using mfLab and directly coupling the compiled MODFLOW Fortran code instead of calling GMS in the developed multi-objective optimization management model that can be used for distributed modeling of the aquifer. In this study, due to the simplicity of using GMS model, its widespread use by researchers, easy communication with MATLAB coding environment, and proper execution speed, an approach based on simultaneous calling of the GMS model in the multi-objective optimization model has been used. The application of this method on a hypothetical aquifer has been investigated by Majumder and Eldho (2015).
The steps of aquifer simulation using the proposed MAT-LAB code are presented in Fig. 3 flowchart. According to this figure, before establishing a link between the simulation model and optimization, it is necessary to be prepared a conceptual model of the studied aquifer based on aquifer geometry and hydrogeological data.
In the conceptual model of this study, a single-layer, unconfined aquifer with an area of 176 km 2 , which is located in the alluvial fan of the Karaj River, is considered. The main constituents of aquifer alluvial sediments are composed of different proportions of destructive elements such as clay, sand, gravel and rubble, and generally, coarser-grained elements make up a larger percentage of soil texture at the inlet and northern areas of the plain. The particle diameter gradually decreases by distancing from the deposition axes towards the outlet and south of the plain. The results of geophysical studies show that the thickness of alluvium in the northern areas is about 50 to 100 m and in the southern and central parts is about 250-300 m. According to the results of pumping tests, the amount of aquifer transmissibility coefficient in the northwestern parts of the plain is about 2500 m 2 ∕day , which gradually decreases towards the middle and outlet sections. The average amount of storage coefficient in Karaj aquifer is equal to 6%. Based on groundwater level variation, the slope of the hydraulic gradient in the northern half of the plain is higher than the southern half of the plain and this slope decreases towards the south. It should be noted that the average hydraulic slope of the aquifer is 1%, which is relatively high.
The boundaries of the modeling area are determined mainly by the spatial distribution of the observation wells (piezometric wells). The northern borders of the study area adjacent to the southern heights of Alborz and its eastern border cover the Karaj River. In order to determine the direction of groundwater movement and also to study the possible inflow and outflow fronts of groundwater, boundary conditions were determined using the hydrogeological feature and prepared the GWL map in ArcGIS10.2 software. In terms of boundary conditions, considering that the trend of groundwater movement from north and northwest to south and southeast and in the direction of Karaj River, so the north and northwest borders were considered as the inflow to the groundwater and the southeast and south borders as the groundwater outflow (Fig. 4). The eastern and western borders of the plain coincide with the streamlines and act like impenetrable boundaries due to the groundwater flow cannot move perpendicular to these lines.
After determining the boundaries of groundwater inflow and outflow, these boundaries were introduced to GMS tools Fig. 3 The proposed MATLAB code structure to establish a link between GMS model and multi-objective optimization algorithm as General Head Boundary (GHB) due to the uncertainty of the volume of inflow and outflow from the boundaries of aquifer. In fact, using this package and according to the changes in the groundwater level inside and outside of boundary and the calibrated hydraulic conductivity from steady conditions, the flow values at the boundaries are calculated automatically.
For qualitative modeling, based on previous studies, the value of the effective molecular diffusion coefficient for nitrate was considered equal to 0.00005 m 2 ∕day . To assign the initial values of nitrate concentration, nitrate point data related to the 2010-2011 water year were used. Also, in order to determine the quality boundary conditions, since the source of nitrate is from the absorption wells of wastewater disposal, so no flow containing nitrate concentration enters the border of aquifer and the boundary conditions were considered as boundaries without quality flux. According to the available data and in order to adapt to the quantitative model, qualitative data measured in 38 wells that had a good distribution in the plain were used. In order to investigate the qualitative behavior of nitrate concentration with observational values, the calibration process of the qualitative model was performed to determine the longitudinal diffusion coefficients and the amount of nitrate concentration entering the aquifer.
After collecting the required information and preparing a conceptual model, it is necessary to design a groundwater flow model for implementation. In this stage, the modeling boundary, number of stress periods and time steps, type of aquifer boundary conditions, cellular amounts of recharge and discharge, and other aquifer hydraulic parameters are determined. In order to groundwater system simulation with the finite difference method, it is necessary to divide the aquifer area into a smaller number of zones which is called a cell. According to the geological condition, topography, area of case study, and the amount of available data from the Karaj plain, the grid with cells 500 × 500 m and 44 rows and 39 columns, and containing 1017 active cells in the UTM geographical coordinate system was prepared.
In order to investigate the effects of reducing cell size on the accuracy of the results of the aquifer simulation model, the model cell size was reduced from 500 to 250 m. Based on the monthly simulated GWL under different cell dimensions, it can be found that the GWL in active operation wells (used in the first objective function of the developed optimization model) does not show significant variation with decreasing model cell size (Fig. 5). This is due to the lack of hydrogeological and hydrological data appropriate to the finer mesh.  Also, since the GW simulation model is used directly in the optimization model, in this study, reducing the mesh size can increase the computational burden without leading to a significant increase in the accuracy of the results.
Other parameters that have been introduced to the GMS tools for simulation of Karaj aquifer are the digital aquifer layers such as topography, bedrock, aquifer hydrodynamic coefficients, land use, monthly piezometric head, location and value of well discharge, initial nitrate concentration, recharge amounts by surface water resources, agricultural return flow, domestic absorption wells and precipitation, the hygrometry and meteorology stations data, etc. As an example, the topographic cell map, bedrock, initial GWL, the location of the operation wells, and the position of the piezometers are presented in Fig. 6.
By performing calibration process, the hydrodynamic coefficients of the aquifer and the longitudinal diffusion coefficients are calibrated and used for validation. In order to use the validated model during the optimization process, it is first necessary to call and store data related to the location and amount of withdrawal from wells based on the validated simulation model. For this purpose, it is necessary to use the following code in MAT-LAB environment: Dwells = hdf5read('VerificationKarajAquifer.h5','/ Well/07. Property'); In this code, "VerificationKarajAquifer" is the name of the validated groundwater simulation model.
The amount of decision variables undergoes many variations during the optimization process to achieve the optimal value. For this purpose, it is necessary for each variation in the decision variables, the original value (Dwells) replaced with new withdrawals values from wells and run a validated groundwater simulation model. H5D.close(dset_id); H5F.close(fid).
With the implementation of QQ simulation model using the following command, the required parameters to calculate the value of the three defined objective functions, namely, monthly GWL and nitrate concentration, are extracted for each active cell: • Command to call the executable file of the aquifer quantitative simulation model: command = 'D:\KarajModel_MODFLOW'; [status,cmdout] = system(command); system('mf2k_h5_parallel.exe VerificationKarajAquifer.mfn'); • Code to extract the amount of simulated GWL: fid2 = fopen(VerificationKarajAquifer.hed','r'); GWL = fread(fid2,'float32 = > float32'); • Command to call the executable file of the aquifer qualitative simulation model: command = 'D:\KarajModel_MT3DMS'; [status,cmdout] = system(command); system('mt3dms53.exe erificationKarajAquifer.mts'); • Code to extract the amount of simulated nitrate concentration: fid = fopen('D:\ KarajModel_MT3DMS\MT3D001. UCN','r'); No3_concentration = fread(fid,'float32 = > float32'); According to the described approach, the coupling process between the flow and transport with the optimization model, which needs to be performed continuously and sequentially, is as follows: • Determining random values of decision variables generated by the NSGA-II optimization model and within the allowable range (as a feasible solution). • Applying the decision variables in MODFLOW calibrated code to simulate groundwater level. • Nitrate concentration simulation using MT3DMS code based on the results of aquifer quantitative simulation model. • Calculation of objective functions based on cellular simulated groundwater level and nitrate concentration, and investigation of the proposed model constraints. • Execute the above three steps for all generated feasible solution in the decision space and determine non-dominate solutions as the optimal trade-off curve between the objective functions. Fig. 7 The map of calibrated hydraulic conductivity ( m∕day)

Results
In this section, initially, the results of the aquifer QQ simulation model and then the analysis and discussion of the optimization management model results are presented.

Results of aquifer simulation model
By implementing a prepared quantitative simulation model under steady conditions and based on the GWL of observation wells for September 2010, first by manual method (trial Fig. 8 The bar chart diagram (a) and scatter plot (b) between the observed and simulated GWL in piezometric wells under steady condition and error), the initial hydraulic conductivity was somewhat calibrated. Then, according to the appropriate capabilities of GMS software in the calibration of numerical models, for calibration of hydraulic conductivity (K) values, a number of pilot points are defined in the model, and based on the values obtained from manual calibration, an initial value of K was assigned to each of these points.
In the next step, the software calculates the hydraulic conductivity values for all model cells by interpolating the initial values given for the pilot points and simulates the GWL distribution in the study area by implementing the model using these values.
Finally, by comparing the observational and computational GWL values at the piezometric wells, the computational error of the model is determined and the model tries to provide a better description of the distribution of this parameter in the study area by modifying the hydraulic conductivity values at the pilot points.
The interpolation method used for the hydraulic conductivity of pilot points is the kriging method, which has more capabilities than other existing interpolation methods (such as IDW) and provides better control over the output results of the interpolation process. The calibrated values of the hydraulic conductivity for the modeling region are shown in Fig. 7. Also, the scatter plot and the bar chart diagram between the observed and simulated GWL in piezometric wells in order to assessment the accuracy of the calibrated model results are drawn in Fig. 8. As shown in this figure, overall, relatively good agreement between the observed and simulated GWL are found in all piezometers.
For the simulation in unsteady condition, it is necessary to construct an unsteady model so that the temporal variation of the aquifer is assessed. The time period considered in this model is one water year (2010)(2011). In this model, the number of stress periods is 12. All hydrological parameters for different stress periods are assigned to cells of aquifer according to the data available in different months. By implementing a quantitative simulation model under unsteady conditions, the spatial distribution of the calibrated specific yield coefficient can be presented in Fig. 9. To evaluate the accuracy of the unsteady calibration results, the GWL hydrograph in the studied piezometers are presented in Fig. 10.
In order to evaluate the prediction accuracy of the calibrated GMS model, the following statistical performance indices were used: where x mi and x ci are the measured and simulated values, respectively. Also, x mi is the average of measured values. The normalized root mean square error (NRMSE) facilitates the comparison between models with different scales. The NRMSE relates the RMSE to the observed range of the variable. Thus, the NRMSE can be interpreted as a fraction of the overall range that is typically resolved by the model. In all the above error indicators, the values closer to zero show that the model performance is more appropriate.
According to Fig. 10 and Table 1, it can be found that the aquifer parameters have been well calibrated in order to model the real conditions governing the groundwater system of the Karaj plain.
For modeling of nitrate concentration's monthly variations, the MT3DMS model was used. Using the manual method and changing the longitudinal dispersion parameter and the amount of nitrate entering the aquifer, the qualitative model was calibrated. Figure 11 shows the calibrated value of the aquifer longitudinal dispersion coefficient. According to this figure, the value of longitudinal dispersion coefficient in the aquifer varies from 0.05 in the central and western parts to 0.26 in the southeastern parts of the aquifer. Fig. 9 The map of calibrated specific yield coefficient (dimensionless)

Results of optimization model
In this study, an NSAGA-II multi-objective algorithm has been used to achieve the optimal trade-off curve between objectives. Parameters related to crossover and mutation operators were determined using trial and error method. Also, the tournament operator was used to select the parent's chromosomes. To reduce the execution time of the optimization model, an initial population with feasible solutions was used instead of completely random solutions. In the NSGA-II algorithm, the population size was considered to be 150. By implementing the developed three-objective simulation-optimization model on a computer with 16 gigabyte RAM and CPU core i7-9700, the optimal trade-off curve was determined.
Regarding the computational costs of implementing the proposed model, it can be stated that based on the According to optimal trade-off curve, the minimum and maximum value of first objective function were 3554.8 m and 3742.5 m (equivalent to an average drawdown of 1.45 m and 1.53 m per well) in total planning horizon, respectively. The minimum and maximum values of a second objective function were 3,347,638.5 mg∕l and 3,352,001 mg∕l (on average, 56.87 and 56.94 mg∕l nitrate concentration per well and monthly), respectively. Also, the range of third objective function is estimated between 226.11 and 230.92 MCM, respectively. These minimum and maximum values of objective function are used to determine the priority of each of the solutions on optimal trade-off curve using MCDM methods.

Extraction of the superior solution
To determine the rank of each solution located on the optimal trade-off curve, seven MCDM methods called WAS-PAS, COPRAS, TOPSIS, CP p=∞ , CP p=2 , and CP p=1 have been used. By applying these decision-making methods to the optimal solutions, the rank of each solution based on each MCDM methods was determined as Table 2. According to this table, the ranking for the solutions is very different with miscellaneous MCDM methods, and it is not possible to provide the final rank. For this purpose, the BAM method was applied to select the final rank that has a higher Berda scoring.
Based on NSGA-II result, it can be found that 139 and 11 solutions (of the 150 generated solutions) are considered as non-dominate and dominate solutions. Therefore, ranking is done on 139 non-dominate solutions and if a solution has a rank of one, its Berda scoring will be equal to 138. Similarly, this Berda scoring can be easily calculated for other solutions. By performing this process for each MCDM method and extracting the sum of Berda scoring obtained for each solution, the final Berda scoring of solutions is determined. By sorting descending of scores, the final rank of each solution in the form of Fig. 12 can be specified. Based on this figure, it can be seen that solution 108 is in the first rank (selected scenario) in terms of satisfy three proposed objective functions simultaneously. Accordingly, the optimal decision variables in proportion to this solution, which contains optimal values of withdrawing from operation wells, are evaluated as a desirable alternative compared to another alternative located on the optimal trade-off curve. Based on the results of alternative number 108, optimal amounts of groundwater abstraction from the aquifer for sustainable QQ development can be extracted. Also, based on these values, the QQ analyses are carried out on the status of each well to determine the effects of the proposed approach on improving the QQ status of the aquifer.

Investigation of the aquifer quantitative status under optimal abstraction conditions
Based on the selected scenario and considering the optimal allocation of existing wells, it can be concluded that for establishing stability in the QQ status of the aquifer, it is necessary to reduce the current abstraction (471.55 MCM) to 228.49 MCM (with a 51.55% reduction) over the planning horizon (Fig. 13). By applying the proposed approach, the quantitative behavior of the aquifer has dramatically improved, so that the reduction of the pumping of the wells has led to an increase of 4.6 m in GWL (an average of 19 cm per month) over 3 years (Fig. 14).
The response of the aquifer to the reduction of pumping is indicative of the high sensitivity of the aquifer to the stresses on it. Therefore, in order to supply the demands of water in the study area, it is necessary to plan other available surface water resources such as increasing the allocation of Karaj dam, increasing the amount of water transferred from Taleghan Dam. Also, in order to achieve aquifer sustainable development, water consumption must be decreased in different sectors and water use efficiency increased in the agricultural section. To observe the status of GWL rise under optimal abstraction conditions, the GWL hydrograph in the operation wells has been drawn in different positions of aquifer according to Fig. 14.
The results show that in most parts of the aquifer, the GWL is rising, and this is especially significant in situations where the number of wells is high. This is due to increase in the saturation thickness of the aquifer as a result of reduced withdrawal.
Histogram analysis of GWL variation in 2453 operation wells shows that with the implementation of optimal withdraw policies, 24.13% of wells (592 wells) experience an increase of 4.5-5 m in GWL during the study period. Also, 13.82% and 13.98% of operation wells, which are   Fig. 15, will have a GWL rise. It is worth mentioning that after applying the optimal allocation results from wells, in 7 wells the GWL increases to more than 12 m in 3 years.
Since the Karaj aquifer area is 254.25 km 2 , an increase of 4.6 m in the GWL, including a specific yield of 0.163, will lead to the annual addition of 95.3 MCM of water to the saturation thickness of the aquifer. This amount, which is equivalent to 41.7% of the total optimal withdrawal from wells, can lead to an increase in the static reservoir volume of the Karaj plain aquifer in the long term.

Investigation of the nitrate concentration under optimal abstraction conditions
The results of the extracted from the best scenarios are analyzed to investigate nitrate concentration. Initially, the general process governing the qualitative status of the aquifer is described after applying the optimal pumping policies, and then details related to the qualitative variations made on wells during the planning horizon are presented.
A study on the time series of nitrate variations over the 3 years (2012-2014) shows that despite a significant reduction in the pumping of wells, the average reduction in nitrate concentration was about 3% (Fig. 16). This is due to the severe pollution of the Karaj plain aquifer as a result of the entry of urban and agricultural wastewater. If these optimal operating conditions persist, this reduction in nitrate concentration can be intensified by utilizing the municipal wastewater collection system over a period of 10 years and can be reduced to less than the permissible limit of nitrate in drinking water (that is, 50 mg∕l according to the World Health Organization (WHO) guideline). In other words, optimal withdrawal policies can lead to a significant improvement in the aquifer quantitative sustainability in the short term, but to create the desired qualitative condition, more time is needed.
Qualitative zoning of the nitrate parameter under the conditions of optimal allocation and continuation of the current withdrawal process was drawn to evaluate the effectiveness of the proposed approach in improving the qualitative aquifer conditions (Fig. 17). Investigation of nitrate concentration variation after applying the optimal operation policy shows the nitrate concentration in the northern, western, and eastern parts of the aquifer has been greatly reduced and it is in a more favorable condition. Continuation of optimal allocation policies can lead to quantitative stability of aquifer short-term and improve the nitrate concentration status. By calculating the levels covered by each of the nitrate concentration ranges (based on the nitrate zoning map shown in Fig. 17), it can be found that the zones with high concentrations of nitrate in the current conditions gradually replaced by zones with lower concentrations of nitrate and the general qualitative conditions of the plain are moving towards zones with low concentrations of nitrate. This is especially evident in areas with a large number of operation wells. For example, using the optimal values allocated to each well, the area of aquifer with a nitrate concentration of more than 110 mg∕1 has been reduced from 5.38 km 2 to 0.67 km 2 (Fig. 18).
In order to evaluate the proposed approach efficiency, the time series of nitrate concentration for the two different operating conditions (continuation of the current situation of aquifer operation and applying proposed model) were drawn in Figs. 19, 20, and 21. Based on these figures, the effectiveness of the developed structure in creating the qualitative stability of the aquifer is quite evident. Recommendation to continue the optimal operation of wells can be significantly effective in reducing the concentration of nitrate, due to the decreasing slope of this parameter.

Conclusion
This paper developed a coupled optimization-simulation model for QQ sustainable management of groundwater resources in an arid and semiarid area using GMS and NSGA-II tools. In the proposed model, three objective functions were formulated to minimize the sum of drawdown of GWL in drinking wells, minimize the sum of nitrate concentration in cells containing operation wells, and minimize the sum of withdrawal rate from wells during planning horizon. Due to appropriate accuracy and widespread use of the GMS model in investigating the QQ behavior of the aquifer, code was prepared in MATLAB environment with the aim of establishing link between this software and multi-objective optimization algorithm (NSGA-II). With implementing the developed model, aquifer optimal operation policies were extracted using MCDM methods.
Analysis of optimal allocation values of wells shows that in order to create sustainability in the QQ status of the aquifer, it is necessary to reduce the total amount of aquifer withdrawal from 471.55 to 228.49 MCM over the planning horizon. This reduction in abstraction has led to an average increase of 4.6 m in GWL, which adds 95.3 MCM of water to the static reservoir volume of the Karaj plain aquifer. Also, the nitrate concentration variation after the implementation of the proposed approach shows the area of aquifer zones with high nitrate concentration has decreased and the quality status of the aquifer have improved. For example, the area with a nitrate concentration above 110 mg∕1 has decreased by 87.5% and reached less than 0.67 km 2 . Findings obtained from the application of the proposed approach shows the simultaneous application of the NSGA-II and GMS models, and MCDM methods using the developed MATLAB code can be successfully used to manage complex aquifer systems. In fact, this approach can be used for any area, but areas whose main demands are groundwater resources and there is no sustainable access to surface water, the results of the developed model can make significant changes in the modification of groundwater level decreasing trends and improve the quality of the aquifer. In order to increase the accuracy in preparing well operation policies, it is recommended to define a smaller mesh size especially around the well field for aquifer modeling.