Development of an integrated multi-objective approach to formulate optimal harvesting policies with the aim of sustainable management of groundwater resources: study area: Varamin Plain

Lack of water in many parts of the world has created significant problems for the provision of safe drinking water, the production of agricultural crops and the general cycle of human life. Also, over-harvesting, drilling of unauthorized wells and insufficient monitoring of the amount of harvesting outside of the exploitation license have also led to a drop in the aquifer level and land subsidence in many plains of the world (Rejani et al. 2009; Zhang et al. 2014; Xiang et al. 2020).

Considering the climatic conditions governing most of the plains of Iran, farmers’ uncontrolled harvesting of underground reserves along with the return period of severe droughts, especially in the Varamin Plain (Azizi et al. 2021; Azizi & Nejatian 2022), caused irreparable damage to the groundwater resources of the plains of Iran. It should also be noted that, due to the many droughts and the increase in water demand all over the world, it is not possible to create sustainability in the quantitative conditions of groundwater resources in many arid and semi-arid regions. However, it is possible to minimize the unfavorable effects of over-harvesting with the optimal management of water resources and the formulation of appropriate policies.

Nowadays, transitional optimization methods for the use of available resources have been introduced to create sustainable conditions for water resources and maximize benefits by using simulation–optimization tools that can provide managers and planners with optimal solutions (Maier et al. 2014). These models use an aquifer simulation model and an optimization tool to solve the management problems of groundwater resources, so in these studies, the optimization tools are based on the number of decision-making variables, the complexity of the problem and the number of goals is different. By implementing these models, it is possible to formulate scenarios for the optimal management of the groundwater table (Esteban & Dinar 2013; Farhadi et al. 2016; Banihabib et al. 2019; Norouzi Khatiri et al. 2020; Sabzzadeh & Shourian 2020). In this section, we review some of these methods.

Rejani et al. (2009) proposed a non-linear management model and a linear optimization model to optimize the pumping rate of exploitation wells for groundwater management in the Balasore coastal basin. Both of the models were developed to optimize the amount of pumping and determine the optimal cultivation patterns to maintain the aquifer level in the desired range. The results obtained based on the policies formulated for wet, normal and dry years showed that the pumping programs and cultivation patterns in three scenarios had a significant difference and the level of groundwater in optimal conditions compared to the current conditions significantly increased by 257% (in the wet period), 167% (in the normal period) and 112% (in the dry period) annually.

Rezapour Tabari & Soltani (2013) presented a multi-objective optimization model for integrated operation management using a non-dominate sorting genetic algorithm (NSGA-II) and steady state genetic algorithm (SSGA). In this study, the modeling is done in such a way that the violation of the reservoir capacity during operation is minimized. The results of the formulated model indicate that the NSGA-II model can provide optimal allocation values in a much shorter period of time. Joodavi et al. (2015) formulated the development and application of the integrated cultivation and exploitation optimization model with an emphasis on the uncertainty parameters of the aquifer to prevent the pressure on the groundwater resources in the Firozabad plain of Fars. The results revealed that the current cultivation pattern should be changed and the cultivation levels of some crops should be significantly reduced to control the groundwater resources.

Alizadeh et al. (2017) developed a multi-objective fuzzy approach based on MODFLOW and MT3D simulation models, the NSGA-II algorithm and the fuzzy transform method (FTM) to determine optimal policies for water management and environment of the Kavar-Maharloo Aquifer. The results revealed the appropriate performance of the proposed model to determine the most sustainable allocation policy in the management of groundwater resources. Sabzzadeh & Shourian (2020) developed an approach using the MODFLOW mathematical code and the particle swarm optimization (PSO) optimization algorithm to maximize the net profit of crop production, considering the sustainability of the aquifer in the Asman Abad Plain in the west of Iran. Their results revealed that by using the proposed approach, the net profit of agriculture will increase significantly compared to the current situation in the plain, while the aquifer changes remain within a sustainable range.

Song et al. (2020) introduced the simulation–optimization approach in Northwest China with the MODFLOW-NWT model with the aim of preserving groundwater resources and minimizing the water transfer cost through the ε-MOMA algorithm. Their results revealed that the reduction of runoff will have a great negative impact on management goals, so the integrated management of surface and groundwater in northwest China will be of great importance. Majedi et al. (2021) simulated the aquifer and surface waters of the Karun Basin with the help of the Modflow mathematical model and the water evaluation and planning (WEAP) programming model with the aim of minimizing the drop in the aquifer of the basin with regard to meeting the maximum needs of the area by the NSGA-II optimization algorithm. Their results revealed the proper performance of the proposed approach in meeting the desired goals until the time horizon of 2040.

Chakrai et al. (2021) used the NSGA-II optimization algorithm and the WEAP planning model to optimize the allocation of surface and groundwater resources in the Zayandeh River basin. The results of their proposed approach showed the sustainability of the groundwater table up to 16%. Saghi-Jadid & Ketabchi (2021) managed and planned to improve the Namdan Aquifer in Fars province using a simulation–optimization approach by the restricted Boltzmann machine (RBM) method and the MODFLOW simulation model. The results of their proposed approach showed an increase in the groundwater level up to 10.61 m in a 10-year time interval.

The review of previous studies in the area of optimal management of groundwater resources shows that in the developed models, the use of an aquifer simulation model to better describe the behavior of groundwater resources and adopting management strategies in aquifer harvesting and control can be highly effective in improving the quantitative conditions of the aquifer. However, the review of the studies shows that most of the studies have been conducted with a single objective: to determine the optimal harvesting of water resources. Also, the use of the NSGA-II multi-objective optimization algorithm and the groundwater modeling system (GMS) aquifer simulation model, along with the use of multi-criteria decision-making methods (MCDMs) that simultaneously determine the value of the optimal harvesting scenarios for each active well in the direction of sustainable exploitation of groundwater, have not been reported in the records.

Accordingly, providing an approach with the aim of sustainable management of groundwater resources and meeting the water requirements of its sub-sectors is necessary and has been considered the main goal of this study. Accordingly, the multi-objective management model was developed based on the simulation–optimization approach. It should be noted that in this model, the accuracy of the quantitative simulation of the aquifer has a direct effect on the results of the optimization model. Thus, the MATLAB environment was used to write the code to directly connect the optimization model with the simulation model with direct access to the input and output files of this model. Finally, five MCDMs were used to determine the rank of each solution and the Borda aggregation method (BAM) was used to select the best scenario to provide optimal policies for aquifer management.

It should be noted that, based on the research approach and to investigate the quantitative behavior of the Varamin aquifer, first, a simulation model of the aquifer should be prepared and then calibrated and validated. The time period considered for aquifer modeling is 2014–2015. The reason for selecting this water year is the completeness of the hydrogeological data of the aquifer. It should be noted that for the quantitative modeling of the Varamin Plain aquifer in steady and unsteady conditions, quantitative data related to the 2014–2015 water years were used for calibration and a 30-month period was used for validation after the calibration period.

The preparation and development of scenarios for the exploitation of groundwater resources for sustainable harvesting of these resources require the definition of goals for optimal management as a first step. For this purpose, in this study, a new practical approach based on the most suitable simulation and optimization tools was developed. To achieve the goals defined in this study, first, it is necessary to present the objective function and limitations of the problem mathematically. The objective functions of the study are to minimize the total drop in the groundwater level in each of the cells of the Varamin Plain aquifer modeling range during the study period and to minimize the amount of water harvested from the exploitation wells in the aquifer.

Given the importance of this plain in the supply of agricultural products in Tehran Province and the necessity of adopting correct policies in line with the sustainable exploitation of the aquifer and preventing destructive processes such as subsidence and erosion, these goals are considered with the policy of simultaneously controlling the amount of harvesting and the reduction in the saturation thickness of the aquifer. It should be noted that the simulation using the GMS distributed model is used to determine the quantitative behavior of the groundwater table in the Varamin Plain. Unlike lump models, it leads to an increase in the accuracy of calculating the spatial and temporal levels of groundwater. The objective functions and limitations of the problem in the desired approach are mathematically as follows:

The variables presented in the above equations are defined as follows: Z₁ refers to the total drop in the level of groundwater in each of the cells of the Varamin Plain aquifer modeling range during the study period (m); Z₂ refers to the total amount of water harvested from exploitation wells in the study area (m³ per day); ΔH_tc refers to the amount of drop in the groundwater level related to cell c in the month t (m); Qwell_tw refers to the amount of water harvested from well w in the month t (m³ per day) (decision variable); QCwell_tw refers to the current situation of harvesting from the well w in the month t (m³ per day); R_tc refers to the amount of natural feeding of the aquifer in cell c from month t (m per day); H_tc refers to the level of groundwater related to cell c in the month t (m); f is a function based on which the quantitative state of the aquifer is simulated. In fact, to calculate the level of groundwater, the MODFLOW distributed model is used as the level simulation function inside the optimization model; nw refers to the number of wells in the study area; m refers to the number of months of the planning period; nc refers to the number of active cells located in the simulation range of the Varamin Plain aquifer.

In Equations (1) and (2), our objective functions are presented in mathematical form. In Equation (1), the first objective, which is the total drop in the level of groundwater in each of the cells of the Varamin Plain aquifer modeling range during the study period, is considered. For this purpose, it is necessary to first simulate the time series of the groundwater level using a calibrated aquifer model. This rate is based on Equation (5) (this equation works based on the MODFLOW simulation model). For each of the solutions provided by the optimization algorithm, it is necessary to run the aquifer simulation model to determine the level of groundwater and the amount of subsidence in each of the simulated cells. In other words, it represents the quantitative behavior of the aquifer in the face of exploitation resources (wells) and its control by a defined constraint (Equation (4)) can play a significant role in the sustainability of the groundwater table. In other words, controlling the monthly drop of the groundwater level in each of the active cells defined in the aquifer level prevents the reservoir deficit and the reduction of the saturation thickness of the groundwater resources, in addition to reducing the pumping costs. In Equation (2), the second objective function of the proposed approach, which is to minimize the amount of water harvested from the exploitation wells in the aquifer, is realized.

Since the current situation of the exploitation of the aquifer has caused an exponential drop in the aquifer, this over-harvesting is controlled by Equation (3). This objective function plays a significant role in the quantitative sustainability of the aquifer. In other words, in this objective, maximum provision of the area requirements is not considered and long-term exploitation of the aquifer and paying attention to the sustainability of the water table to improve the quantitative situation is the priority of the harvesting policies because there are significant surface water resources in this area that if restrictions are imposed on harvesting from underground reserves, these resources can be considered to compensate for the reduction in the harvesting. By using these two defined objectives, which are complementary in the sustainable management of the aquifer, it is possible to determine the optimal amounts of allocation from each of the wells in the study area and based on that, the necessary plans can be adopted for meeting the water requirements of the area from other water sources (such as surface water reserves).

To solve the developed optimization model, it is necessary to first define the variables of the problem, which are the decision variables or unknowns of the model. In this study, the decision variable considered for each month is the amount of water harvested from the wells in the aquifer level of Varamin Plain. Since there are 2,058 active wells in this plain, the total number of decision variables of the problem will be 24,696 considering 1 year study period. It should be noted that the reason for considering each of the exploitation wells in the plain as a decision variable is an independence of their exploitation, which is mainly managed by the private sector and cannot be integrated regionally.

Then, the quantitative simulation model of the aquifer was prepared and calibrated by GMS 10.3 software. It should be noted that the input data of the GMS simulation model can be called from the GIS (geographic information system) software. Thus, in the present study, the digital layers of the aquifer, such as topography, bedrock, piezometer data, hydrodynamic coefficients of the aquifer, land use and climatic data were prepared in a GIS environment and introduced to the GMS model.

Since the correct simulation of the quantitative behavior of the aquifer is very important in the process of optimization and determining the optimal values of the decision variables, a code was developed in the MATLAB software environment to increase the accuracy of the simulation results of the aquifer. This code, which can call all the input and output files of the GMS software, can perform quantitative modeling of the aquifer for a period of 12 months in a short period of time (approximately 1.25 s). In other words, the user can change the state of stresses in the groundwater system and observe the level changes by implementing the GMS model by the developed code. In this code, it is necessary to use the calibrated simulation model that represents the real behavior of the groundwater table in Varamin Plain. For this purpose, the GMS model was first calibrated under steady and unsteady conditions to determine the state of changes in the hydrodynamic coefficients of the aquifer. Then, to ensure the accuracy of the prepared simulation model, the model was validated based on the observation data that were not used in the calibration process.

Due to the complexity of the prepared structure, the dimensions of the decision variables and the need to increase the speed of achieving the optimal solution (the amount of monthly harvesting from each well), it is necessary to use the most appropriate optimization tool. A number of multi-objective evolutionary algorithms (MOEAs) have been suggested in the last two decades and the non-dominated sorting genetic algorithm II (NSGA-II) is one of the promising MOEAs and has been successfully applied in many fields. In this study, NSGA-II, which was proposed by Deb et al. (2000), was used to solve the problems of the classical genetic algorithm. The NSGA-II is an improved version of the NSGA developed to address issues of computational complexity as well as to provide an explicit mechanism for diversity preservation.

To determine the decision variables of the model using the NSGA-II tool, random values are considered for the amount of harvesting from each well (decision variables). Based on these values, the value of the target functions is calculated and in the case of validating the limitations considered for it, the produced solution is revised until the amount of violation is zero. During the process of improving the values of the decision variables, the objective functions are also examined to have the lowest amount among the solutions. This process continues until there is no change in the values of the objective functions. In this case, it can be stated with a very high probability that the value provided for the decision variables is optimal and can be used as a plan for the exploitation of wells in the studied area.

According to these values, it is possible to formulate optimal policies monthly or quarterly. To get familiar with the details of the steps of running this algorithm, you can refer to the article by Deb et al. (2002). Also, for the use of this optimization algorithm for the exploitation of water resource systems, the articles presented by Rezapour Tabari & Soltani (2013), Hajiabadi & Zarghami (2014) and Bozorg-Haddad et al. (2016) are recommended. In this study, the NSGA-II optimization algorithm code was prepared and run in MATLAB 2018b environment to apply the objective functions and the presented constraints. By running the optimization–simulation model, the optimal trade-off curve between the objective functions is extracted. In the optimal trade-off curve obtained using the NSGA-II algorithm, each point is considered as a resource exploitation scenario (exploitation wells) and the decision maker can select one of the scenarios and extract the optimal values based on it. Since the points located on the optimal trade-off curve cannot meet both goals simultaneously, selecting the right point using MCDMs can be helpful to a large extent in selecting the desired scenario (Bolouri-Yazdeli et al. 2014; Ahmadisharaf et al. 2015; Bozdağ 2015; Bozorg-Haddad et al. 2016). It should also be noted that a short-term planning horizon (1 year) was considered in this study to extract the best policies for the performance of wells. It is not practical to provide optimal policies for the exploitation of groundwater resources for a long period of time due to reasons such as changing the exploitation approaches in the organizations in charge of managing Iran's aquifers and insufficient monitoring of aquifer resources.

Since the MCDMs are very diverse and their use may create different results, an approach was presented in this study using BAM to aggregate the results of several MCDMs. For this purpose, five MCDMs called compromise programming (CP), complexity proportional assessment (COPRAS), a technique for order preference by similarity to ideal solution (TOPSIS), modified-TOPSIS (M-TOPSIS) and weighted aggregates sum product assessment (WASPAS) were used. Then, the results of these five methods regarding the ranking of the solutions produced by the NSGA-II algorithm are aggregated using the BAM method and the final ranking of each point located on the optimal trade-off curve is determined. More details of the used method can be found in the article by Banihabib et al. (2017).

In the calibration stage, the input parameters of the model, including hydraulic and hydrodynamic data were adjusted to achieve an acceptable match between the groundwater level observed in the piezometers and the groundwater level calculated by the model in two sustainable and unsustainable ways. The aim of calibration is to adjust the specific discharge coefficient of the groundwater table, calibrate various parameters of the model and minimize the error in successive time steps. In this study, the maximum acceptable error between the observational and simulated water levels of the observational wells is ±2 m. Based on the calibrated model, the quantitative behavior of the aquifer can be simulated for a short period of time. In fact, based on the calibrated parameters, it becomes possible to predict the level of groundwater in future conditions. For this purpose, due to the lack of proper discharge and feeding information from the aquifer in this study, assuming the prevailing conditions related to the 2014–2015 water years, the model was validated for a period of 30 months after the calibration period (according to the existence of information from an observational well) was simulated and its results were evaluated in the form of a hydrograph of the groundwater level.

For the sustainable management of groundwater resources using the proposed model, it is necessary to first calibrate the Varamin aquifer simulation model. Then, by implementing the multi-objective optimization model and harvesting the exploitation scenarios of the groundwater system, it is possible to determine the optimal performance policies of the wells. In the first step, we first examine the results of the quantitative simulation model of the aquifer and then analyze the results of the optimization model.

In addition, based on the values of the simulated groundwater level, the modeling error can be determined on a monthly basis using error indices and the correlation coefficient (Table 1). According to this table, it can be seen that the prepared model can simulate the aquifer well and can be used to predict the future conditions of the aquifer according to the amount of harvesting on the plain.

Based on the objectives defined in the proposed structure, the NSGA-II algorithm was run with 200 chromosomes and 500 iterations (it continues until the minimum difference is reached). The run time for each of the iterations of this model on a system with the specifications of Core i7-9700k and RAM-16G is estimated to be 330 s and for 500 iterations, it is necessary to run the program in about 46 h. By running the proposed approach by the NSGA-II multi-objective optimization algorithm, the optimal value of the decision variables, which is the optimal harvesting rate from 2058 exploitation wells of the Varamin Plain aquifer, was extracted.

In the optimal trade-off curve obtained by using the algorithm (NSGA-II), each point is considered as a scenario of aquifer exploitation and the decision maker can choose one of the scenarios based on the two objectives and extract values based on that. Due to the fact that the points located on the optimal trade-off curve cannot satisfy both goals at the same time, choosing the right point using MCDMs can be of great help in choosing conflicting scenarios. Details of the decision-making methods used can be found in the paper by Banihabib et al. (2017).

In this study, in order to determine the rank of each solution located on the optimal trade-off curve, five MCDMs called WASPAS, COPRAS, TOPSIS, CP_p=3, 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 was determined in Table 2. The position of the chromosomes in the last optimal generation shows the basis of the numbering presented in Table 2. From this generation, the dominant answers are removed from the set of answers and the remaining chromosomes are numbered from beginning to end and are ranked using MCDMs. Therefore, ranking is done on 191 non-dominate solutions.

According to this table, the ranking for the solutions is very different with miscellaneous MCDMs 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 score, in such a way that if a solution has a rank of one, its Berda scoring will be equal to 190. Similarly, this Berda score can be easily calculated for other solutions.

The obtained results show that, on average, 23% of the levels with a drop in groundwater level have decreased monthly, with a change range of 16% in October 2014, and 27.8% in September 2015. Accordingly, by applying optimal exploitation policies and their continuation during a medium term period of several years, it is possible to compensate for the adverse effects caused by uncontrolled harvesting and to control the harvesting from the wells of the aquifer for the sustainable development of the aquifer.

In the present study, an optimization–simulation model in a semi-arid region (the Varamin aquifer) was developed with the aim of sustainable groundwater management. In the proposed model, two objective functions were formulated to minimize the total drop in the groundwater level and the amount of water harvesting from the exploitation wells in the aquifer during the planning horizon. Due to the appropriate accuracy and widespread use of the GMS model in investigating the quantitative behavior of the aquifer, a code was prepared in the MATLAB environment with the aim of establishing a connection between this software and the multi-objective optimization algorithm (NSGA-II). After the calibration and validation of the GMS simulation model in steady and unsteady conditions and using it in the multi-objective optimization model, the groundwater optimization model was implemented. In this study, five MCDMs were used to determine the rank of each solution and the BAM method was used to select the best scenario.

The analysis of the proposed approach indicates that to create sustainability in the quantitative situation of the Varamin aquifer, it is necessary to reduce the total amount of aquifer harvesting from 376.27 million m³ to 199.81 million m³ in the planning horizon. For this purpose, it is necessary that the permitted amount of harvesting from the wells in this area be reduced from 21 to 11 L/s. In other words, the optimal value of the obtained decision variables (harvesting rate from 2,058 wells) is used as the exploitation policies of each of the wells at the aquifer level, the violation of which can result in moving away from the goals of achieving the sustainable development of the Varamin Plain aquifer. These exploitation rules, which can be notified to each of the operators in the form of exploitation licenses and permissions, should be associated with incentive policies to reduce the amount of harvesting from groundwater tables so that they can be implemented by those in charge of groundwater resource management systems. For this purpose, it is necessary to implement this harvesting reduction and legally by including it in the exploitation licenses. Also, the projects of monitoring the wells in operation are necessary for the correct implementation of the exploitation values included in the exploitation license.

It is also possible to use incentives such as low-interest loans to improve their irrigation system to help farmers due to the changes and limitations in the amount of water harvested from wells and consequently affecting the area under crops. Investigating the results obtained from the application of the proposed approach in the quantitative management of the aquifer indicates that the simulation–optimization model shows the high performance of the harvesting policies from the aquifer to create quantitative sustainability in a short-term planning period. In other words, the simultaneous use of NSGA-II and GMS models and MCDMs using the developed MATLAB code can be used to improve the quantitative situation and achieve sustainable management of water resources practically in other aquifers. It should be noted that the proposed approach can be evaluated in the form of other optimization goals, such as the modification of the cultivation pattern of the area or the qualitative management of the aquifer, as well as other MCDMs.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.