Groundwater sustainability under land-use and land-cover changes

In recent decades, agricultural activities have increased water withdrawals from the Shazand plain in the Markazi Province (Iran). In this study, the Soil & Water Assessment Tool (SWAT) model was used to estimate recharge content as an essential component of groundwater models. MODFLOW2000 in the GMS10.5 software was used for groundwater modeling, and the extent of land use change in the Shazand plain was investigated. The results showed that the agricultural sector allocated the largest change with an 18% increase. During 2009–2016, the water table declined by 5 m due to the decrease in recharge and the increase in exploitation. Therefore, the management scenarios of removal of irrigated crops, optimization of cultivated area, and reduction of the cultivated area by 10, 20, and 30% were applied. The results showed 3, 0.28, 0.49, 1, and 1.5 m increases in the water table at the end of the 7-year study period. Every year, reduction of the cultivated area by 10, 20, and 30%, almost 31, 20.5, and 10.3 million m3 less will be extracted from the aquifer.


Introduction
Water stress is among the greatest issues affecting many societies on a global scale. Water deficit is mainly caused by groundwater extraction for irrigation, which can deplete aquifers (Acero Triana et al. 2020;Dalin et al. 2017). As a result of extensive groundwater extraction in the second half of the twentieth century for agriculture, industry, and domestic water supply, groundwater resources are depleting at an alarming rate worldwide (Bartolino 2003). For instance, excessive groundwater depletion and water stress have been reported in northern and southern Africa, the central United States, northeastern China, Australia, and south and central Asia (Acero Triana et al. 2020). This decline causes disparities and complications in different areas, which can be avoided by implementing an integrated water resources management strategy. For an effective program, understanding the hydrological system and the interaction between groundwater and surface water should be integrated with robust modeling tools to make decisions (Izady et al. 2015).
Arid and semi-arid regions are characterized by a lack of precipitation and extensive groundwater harvesting. The interaction between groundwater and surface water is particularly important in semi-arid areas because management changes in one can affect the availability and use of the other (Fleckenstein et al. 2010). Groundwater is the main water source in arid and semi-arid regions of the world. Overall, groundwater provides 43% of irrigation water, 36% of domestic water, and 24% of industrial water (Döll et al. 2012). Iran is a developing arid and semi-arid country facing considerable population growth. As a result, groundwater resources are overused in this country. In recent years, extensive studies have been done for the integrated management of water resources with the help of SWAT integration with other models (such as SWAT-REMM, to improve simulations of riparian buffer zones; SWAT-SWMM, storm water management; SWAT-SOBEK, sediment and hydrodynamic flow simulation) and SWAT with MODFLOW. This study has used the Soil & Water Assessment Tool (SWAT) (Arnold et al. 1993(Arnold et al. , 1998Arnold and Fohrer 2005) and MODFLOW (McDonald and Harbaugh 1988) as modeling tools to verify and calculate groundwater levels. Dowlatabadi and Ali Zomorodian, (2016) examined SWAT and MODFLOW simultaneously. They estimated the effect of temporal and spatial factors on infiltration values 147 Page 2 of 12 using the SWAT model. These values were then passed to the MODFLOW model in two steady and unsteady states. Overall, their integration results were found to be acceptable. Similarly, Nasiri et al. (2022) used SWAT and MODFLOW to evaluate the effects of irrigation return flow on groundwater recharge and found that irrigation affected groundwater interactions in the basin. They also indicated that the combined use of these two models could be employed in areas with deficit groundwater sources where it is difficult to evaluate the interaction between surface and groundwater sources. Unlike the SWAT model, which cannot accurately represent groundwater, this model is accurate in estimating the recharge rate calculated using flow data and cultivated area. In groundwater modeling, the accurate recharge rate in the input data is one of the most important components that MODFLOW cannot calculate accurately. Therefore, the results of groundwater flow simulated by this model are subject to significant uncertainties (Izady et al. 2015;Kim et al. 2008).
Based on current withdrawal rates, groundwater resources' sustainability is being questioned. Land use/land cover (LULC) changes have been reported to have a greater impact on the hydrological cycle than climate changes (Gyamfi et al. 2017;Vörösmarty et al. 2004). Understanding the impact of groundwater recharge on land cover is critical to pursuing accurate policies to ensure groundwater sustainability under LULC change (Gyamfi et al. 2017). In this regard, the water table is expected to improve with the conversion of crops to pasture as the consumption of irrigation water and groundwater extraction is reduced (Acero Triana et al. 2020); therefore, this study hypothesized that the conversion of crops to pasture has a positive effect on water table in this region too. The present study aims to develop surface and groundwater models and the integrated management of the watershed with different scenarios to revitalize the aquifer by changing the cultivated area. The research objectives are described below: (1) Investigating the changes extent in the cultivated area over the last 10 years and their effect on infiltration levels.
(2) Investigating the water table decline extent over the last 10 years. (3) Investigating the management scenarios' effect on the water table, considering the economic situation of farmers.

Study area
The Shazand watershed, having an area of about 1100 km 2 , is located in the southwest of Markazi province (Iran) with geographical coordinates of 49° 27′-49° 30′ E and 33° 54′-33° 57′ N (see Fig. 1). It has a semi-arid to cold semi-arid climate with a mean annual rainfall of 420 mm (Davudirad et al. 2016). This area's maximum and minimum elevations are 3218 and 1775 m, respectively. Farmlands and gardens have been impressively established in this fertile plain in recent years. According to the census data in 2016, there are 3 cities and more than 150 villages in this region, the largest rural population in this province, moreover main occupation is agriculture and agriculture-related sub-sets. Many important national factories have also been established in this region, such as Imam Khomeini's oil refinery, petrochemical plants, power plants, and other private industries. As a result of industrial growth and urban expansion, water consumption has increased in all three sectors-industrial, domestic, and agricultural. Although the refinery has been a good incentive to increase education, the economic situation is below average because the natives do not have full access to refinery jobs. If there is a chance, it is only for Shazand city and 22 villages at a distance of 2-18 km from the refinery. Some people work in agriculture and industry at the same time, so that portion of participation in the agricultural sector in Shazand city is 17.5%. Since this watershed depends on groundwater resources, the increasing extraction of aquifers has led to a serious decline in these resources. Besides, the water balance in the watershed is strongly influenced by the cultivation pattern.

An overview of the SWAT and MODFLOW models
SWAT is a semi-distributed conceptual model provided by the United States Agricultural Research Service. In large and complex watersheds, this model can estimate various processes such as the effects of land use on water resources, sediments, and fertilization doses and dates over long periods (Arnold and Fohrer 2005;Arnold et al. 1998 where SW t = final soil water content (mm), SW 0 = initial soil water content on a day i (mm), R day = amount of precipitation on day i (mm), Q surf = amount of surface runoff on day i (mm), E a = amount of evapotranspiration on day i , Q gw = amount of return flow on day i (mm), and W seep = amount of water entering the vadose zone from the soil profile on day i (mm).
MODFLOW is a three-dimensional simulation model of groundwater flow presented at the American Geological Survey Center (Abbaspour et al. 2004). This model solves the equations using the finite difference method in the saturation state and establishes Darcy's law (McDonald and Harbaugh 1988). The following partial differential equations (PDEs) describe groundwater flow in three dimensions: where K xx , K yy , and K zz are, respectively, the hydraulic conductivities along the x, y, and z coordinate axes parallel to the major axes of hydraulic conductivities; h is the potentiometric head; W is a volumetric flux per unit volume representing sources (W is negative) and/or sinks (W is positive) of water; S s is the specific storage of the porous medium; and t is time. K xx , K yy , K zz , and S s are functions of space (x, y, z), and W is a function of space and time (x, y, z, t) (Todd and Mays 2005). This model was simulated in two stages. First, steadystate calibration is performed to determine optimal parameters, and then an unsteady simulation is run to estimate the water table and more accurately determine the parameters. MODFLOW2000 was used in this research on the GMS 10.5 software.

SWAT
Data included digital elevation models (DEMs), soil maps, land cover/land use (LULC) maps, and weather information. The required 12.5 × 12.5 m 2 resolution map DEM was obtained from the NASA database (https:// search. asf. alaska. edu). The model created the stream grid and 33 subbasins based on the DEM. It also provides two methods for  (Green and Ampt 1911;Neitsch et al. 2011). Also, surface runoff was calculated using the curve number (CN) method. Soil data were obtained from the FAO soil map with the scale of 1:5,000,000 (Unesco 1970). The texture was mainly clay-loam and loam. Land use information was extracted from Landsat 7 ETM + imagery with supervised classification using the maximum likelihood method. A 30-m resolution image was downloaded from the USGS database (https:// glovis. usgs. gov/), and five classes were extracted according to Anderson's classification scheme (Anderson 1976). The Arak and Shazand synoptic stations provided daily rainfall and minimum and maximum temperatures for the period 2009-2019. Also, the other data such as wind speed, solar radiation, and relative humidity were simulated using the weather generator in SWAT. The model is calibrated using a hydrometric station at the basin outlet. To this end, the data of the Pol-e Doab hydrometric station ( Fig. 2c) were used. Based on the available data, we also used the Hargreaves method (Hargreaves and Samani 1985) for evapotranspiration. This study used Arc SWAT2012 during a simulation period from 2009 to 2019. In 2009, the model was warm-up for simulation and removed from the analysis period. Based on the crops available in the area, the estimated water requirement, and the length of the growth period according to the climatic conditions, flood irrigation was used for crop management. The SUFI-2 algorithm was used to analyze sensitive parameters and calibrate and validate SWAT-CUP software. Calibration and validation were conducted for time periods of 2010-2016 and 2017-2019. Using the P-factor and R-factor, the best simulation results were obtained when the P-factor tends to 1 and the R-factor to 0 (Dowlatabadi and Ali Zomorodian 2016).

Modflow
In this study, the aquifer was considered a one-layer unconfined aquifer with an area of 287 km 2 . The modeling was performed using a grid of 500 × 500 m cells with 55 columns and 52 rows. Also, a DEM and bedrock were used to extract the upper and lower layers by interpolation. Pumping tests were conducted to determine the hydraulic conductivity and specific yield. The steady-state simulation was conducted in October 2009. Also, calibration and validation were done from November 2009 to October 2017 and from December 2017 to October 2019, respectively. Since the water table is relatively deep compared to the ground surface, evaporation was not considered. According to the output of the SWAT model, MODFLOW was supplied with a daily recharge based on the SWAT outputs. To obtain optimal values, 21 piezometers were used to calibrate the model by trial and error.

Performance indicators
In this study, R 2 and NSE criteria were used to estimate the accuracy of the SWAT model, and R 2 and RMSE criteria were used to estimate the accuracy of the MODFLOW model.

SWAT sensitivity analysis
The first step in calibrating the model is adjusting the input parameters so that the simulated results match the observed values (Kumar et al. 2017;Zeckoski et al. 2015). Sensitivity analysis is used to determine the parameters that can most effectively control the Shazand plain river discharge. In Table 1, these parameters have been determined based on previous studies and hydrological knowledge of the region. These parameters were determined based on the highest value of t stat and the lowest value of p-value (Abbaspour 2007). The most influential parameter was recorded by the soil evaporation compensation factor (ESCO), followed by the threshold depth of water in the shallow aquifer required for return flow to occur (GWQMN), initial SCS runoff curve number for moisture condition II (CN2), and snow melt base temperature (SMTMP). In contrast, base flow alpha factor (ALPHA BF) had the smallest effect. The parameters of the method for changing the parameter values listed in Table 1 are as follows:v parameter value is replaced by given value or absolute change; r parameter value is multiplied by (1 + a given value) or relative change (Abbaspour 2007).   Nash and Sutcliffe (1970), an NSE value greater than 0.75 represents a good simulation, while a value greater than 0.36 is satisfactory. Based on the results of this study, the NSE and R 2 values were 0.73 and 0.77, respectively, and the R-factor and p-factor were 0.11 and 0.67, respectively, indicating uncertainty in the conceptual model, parameters, and input data and good model performance. Using validation allows for predicting future conditions without resetting the calibration parameters (Zheng et al. 2012). This validation was performed with the algorithm used in the calibration for a period of 3 years (2017-2019). The performance of the validation period is shown in Table 2. Thus, the R-factor and p-factor were obtained to be 0.0 and 0.13, respectively. Meanwhile, the NSE and R 2 were estimated to be 0.57 and 0.65, respectively, indicating the validity of the developed model.

MODFLOW calibration and validation
The aquifer's recharge rate was calculated using the SWAT model, and these values were used as inputs to the MODFLOW 10.5 model in the GMS package for steady and unsteady states. Calibration was based on a monthly recorded water table. According to Chiang and Kinzalbach (1998), MODFLOW requires an initial water table to begin modeling. For this reason, we assumed a steady state for one month (October-November 2009) in which certain initial parameters were changed by trial and error and then constructed an unsteady state for the remaining months. Markazi Regional Water Authority records 1134 pumping wells in the Shazand plain (Fig. 4). After determining K and zoning hydraulic conductivity (Fig. 5a), the groundwater model was calibrated in an unsteady state for 84 stress periods (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), and the specific yield value was estimated (Fig. 5b). Changes occurred mainly in the northern parts (outlet of the basin), which may be attributed to the diverse geological facies and non-uniform erosion resistance of the basin (Markazi Regional Water Authority 1998). As shown in Fig. 5c, the 10-year mean recharge values for the aquifer vary from 1 to 550 mm. The highest recharge value was recorded in March due to spring rains. In addition, Fig. 6 shows the hydrograph of some unsteady state observation wells. The hydrographs are accurate in most of the observation wells, although some simulated values differ from the observed values in some wells. It may be caused by unauthorized wells, a lack of attention to recording statistics, or a model error in estimating values. In addition, the model was validated for the years 2017-2019. In Figs. 7 and 8, simulated and observed water tables are shown in calibration and validation modes, respectively.

Management scenarios
First, the results of land use maps produced for 2009 and 2019 and their impact on infiltration calculated by the SWAT model were examined. Next, the effects of reducing infiltration and increasing groundwater extraction due to cultivated area growth were explored. Finally, the effects of different scenarios for improving the water table were considered. 1. Land use changes Figure 9 presents land use changes for the 2009-2019 period. There have been changes in all 5 classes of land use, with most of these changes occurring in urban, agricultural, and pasture areas. As shown on this map, urban and cultivated areas have increased by 11% and 18%, respectively, since 2009, while pastures have decreased by 27.1%. Also, the industry sector and barren have faced slight changes. Most of these changes have taken place with the increase in emigration due to the existing industries in the region and food supply. Based on SWAT modeling and a comparison of output results for these 2 years, it is concluded that the annual groundwater recharge has decreased by 1 cm on average. Gyamfi et al. (2017) reported that the annual aquifer recharge level has decreased by a total of 1.271 cm from 2000 to 2013. In this regard, the agricultural sector experienced the greatest change, with a 20.1% increase during the same period. These findings are consistent with the results of this section.
2. Figure 10 depicts the hydrograph of the aquifer studied by piezometers in the region over 7 years (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). The figure demonstrates a downward trend with a 5 m decline. According to the Markazi Regional Water Authority (2017), the withdrawal from this aquifer in the agricultural sector is 231 million m 3 per year. Meanwhile, the needs of this sector are met by only 5 million m 3 of surface water per year. This phenomenon shows that this region is dependent on groundwater.

Removal of irrigated crops
Wheat, barley, beans, and alfalfa, which had the largest area under cultivation, were eliminated by applying this scenario. According to the results, the water table increases by 42 cm every year. As a result, after 7 years, a 3-m increase in the water table is expected. In addition, aquifer extraction declines due to water demand and cultivation area of almost 71 million m 3 annually. Considering that the Shazand plain is dependent on groundwater, there was no significant change in surface water as expected. De Almeida et al.
(2018) evaluated the effects of soil tillage and land cover on water infiltration into the soil. They concluded that soil water infiltration is more influenced by vegetation cover than by tillage. Our results also highlight the effect of vegetation on infiltration rate.

Modification of cropping patterns by reducing water consumption and considering the economic situation
In MATLAB, the area under cultivation was optimized based on each plant's water consumption and income. The best result was obtained in the range of ± 50. Thus, water consumption was reduced, and the income level was decreased by the least extent. Table 4 presents the area under the initial and optimized cultivation. According to this scenario, the water table will grow 28 cm in 7 years, and extraction from the aquifer will decrease by 2 million m 3 over that period. Despite the low increase in the water table, Dowlatabadi and Ali Zomorodian (2016); Kouchakzadeh and Saleh (2014) have shown that the recharge values obtained from the SWAT model used in MODFLOW have good certainty. The explanation for this result is that the total cultivated area is considered constant in this scenario. Ghaffari et al. (2010) studied the impact of  land use change and concluded that pasture and agricultural land have decreased by 34.5% and 14.3%, respectively. Subsequently, these changes have led to a 22% decrease in groundwater recharge. These findings are consistent with the results of this study. Likewise, Sun et al. (2018) identified that land use changes are often considered the main factor for soil infiltration. Also, they showed that soil infiltration rate decreased by 45.23% when pasture land was converted to farmland.

The scenario of removing 30%, 20%, and 10% of the cultivated area
Practical recommendations should be made on changing the cultivated area to be successful with the proposed models. It is, therefore, important to make these changes gently to be applicable (Alipour et al. 2019;Ghodspour et al. 2021). This section aimed to examine the effect of a 10-30% reduction in the cultivated area on the water table. According to Scenarios 30, 20, and 10%, the water table will probably rise by 1.5, 1, and 0.49 m, respectively (see Fig. 11). Every year, almost 31, 20.5, and 10.3 million m 3 less will be extracted from the aquifer. Kibii et al. (2021) found that with the reduction of forest cover and change in land use to urban, aquifer recharge increased by 17% and surface flow decreased by 9%. Overall, it can be argued that the observed increase in the water table and no change in surface runoff are consistent with our findings. According to Lerner and Harris (2009), groundwater provides an excellent water supply for domestic use, and land use significantly affects groundwater resources through changes in recharge rate.

Conclusions
Recharge value is one of the most important components of groundwater. In this study, SWAT and MODFLOW models are used for simultaneous stimulation of the surface and groundwater of Shazand plain. To this end, topography, soil, land use, and crop management were considered essential variables for understanding the status of groundwater resources. In the MODLFOW model, recharge values were estimated using the SWAT model. Both models were calibrated and verified satisfactorily. Overall, the results of the combination of the two models were successful. Since the cultivated area has increased in Shazand plain between 2009 and 2019, the groundwater resources have suffered a negative balance due to a decrease in infiltration and an increase in exploitation. Based on the results of the scenarios, the water table can be increased after a couple of years.