Modeling impacts of climate change and human activities on groundwater resources using MODFLOW

Over the last decades, groundwaters are considered as substantial water resources in many parts of the world, especially in developing countries (Shah 2007; Acharyya 2014). In arid and semi-arid regions, the intensive use of groundwater resources has often affected groundwater levels and water quality (Foster & Louks 2006; Lachaal et al. 2010). The International Panel of Climate Change (IPCC) demonstrated that climate change is occurring on a global scale (Mackay 2008). These changes are expected to include increasing temperature and precipitation patterns (Seiler et al. 2008; Green et al. 2011). There is wide agreement that climate change has caused noticeable fluctuations in the hydrological cycle and will impact on the availability of water resources, either surface water or groundwater. There is evidence that groundwater is an essential parameter of water supply systems in arid and semi-arid areas of the world (Scanlon et al. 2006). Climate change on different time scales, including inter-annual to multi-decadal time scales, appears to affect groundwaters in terms of levels and recharge (Tremblay et al. 2011; Goderniaux et al. 2015). Taylor et al. (2012) demonstrated that groundwater recharge is highly related to projected precipitation changes. Thus, future changes of the groundwater under climate change should be considered more seriously by politicians and managers (Calbo 2010; Boithias et al. 2014; Garrote et al. 2015). Sensitivity of groundwater response to climate change is associated with the aquifer size and groundwater depth, with shallow and small aquifers being more variable to climate change (Lee et al. 2006). To manage groundwater in a sustainable manner, it is essential to conserve groundwater resources for infinite time without any unsuitable economic or social outcomes (Alley et al. 1999). In order to achieve this goal, adequate knowledge of groundwater systems and evaluation of the groundwater resources with regard to the impacts of climate change is necessary (Carrera-Hernández & Gaskin 2006; Chenini & Ben Mammou 2010). Direct measurement of groundwater recharge is difficult due to the heterogeneity and complexity of recharge processes (Kinzelbach et al. 2002). Hydrological modeling is a robust indirect method to quantify groundwater recharge (Jayakody et al. 2014) which is necessary for assessing the effect of climate variability on aquifer recharge. Statistical and numerical models have been considered as crucial tools in identifying groundwater systems and their fluctuations with climate change. For example, the work presented by Lemieux et al. (2015) assessed the impact of climate change scenarios until 2040 on groundwater resources in Magdalen Islands in Canada. Simulations were conducted along a 2D cross-section until 2040 to assess the individual and combined impacts of sea-level rise, coastal erosion, and groundwater recharge on the position of the saltwater–freshwater interface. For simulation of the sea-level rise, the surface boundary conditions were modified to account for the progressive inland displacement of the sea shore. Results of this study showed a decrease in groundwater recharge and increase in the sea level over a 28-year period. Scibek & Allen (2006) applied MODFLOW to investigate impacts of climate change on groundwater recharge in an unconfined aquifer. To assess impacts of climate change, the Canadian Global Coupled Model (CGCM1) and LARS-WG stochastic weather were used. Results showed increases in aquifer recharge from spring to the summer season. In winter, snowmelt infiltration does not happen because of the frozen layer of the ground. The autumn season has moderate recharge compared to the early summer. Investigation of the climate change effects on groundwater storage and recharge has been studied by Jackson et al. (2011) in the UK Chalk aquifer. In this study, the A₂ emissions scenario of general circulation models (GCMs) are used to project climate variations. To simulate the Chalk aquifer, a distributed recharge model and a groundwater flow model were applied. Based on the results of this study, it is concluded that future groundwater recharge ranges between a 26% reduction and a 31% increase until 2080. Döll et al. (2009) used a global scale groundwater recharge model which demonstrated that by 2050, groundwater recharge in north Brazil, the southern edge of the Mediterranean Sea and southwest Africa will be decreased to 70% of potential groundwater recharge under A₂ and B₂ emission scenarios. Crosbie et al. (2013) investigated both decreases and increases in groundwater in the USA. Results of this study showed that the magnitude and direction of the groundwater are uncertain and variable between future dry and wet climate scenarios. A statistical approach named HARTT was applied by Emelyanova et al. (2013) to project future groundwater trend under climate change in the Northern Perth Basin (NPB) in southwest Australia. Results demonstrated that until 2030, groundwater levels were expected to have a slightly rising trend. Shrestha et al. (2016) investigated climate change impacts on groundwater of the Mekong Delta. They used MODFLOW to estimate groundwater level and storage which showed declining trends in groundwater levels towards the future.

The above studies demonstrate that groundwater recharge predictions are variable and depend on the selected groundwater simulation model, geographical situation, and GCMs. To manage water resources in a sustainable way, reliable estimates of groundwater recharge are essential. To achieve this information, groundwater simulation should be considered using appropriate models. Analytical methods do not require many data, but their application is limited to simple problems. Numerical solutions can handle more complicated problems than analytical solutions. With the rapid development of computer processors and increasing speed, numerical modeling has become more effective and easy to use. In our work, a three-dimensional finite-difference groundwater flow model (MODFLOW) developed by Jyrkama et al. (2002) was used to simulate the aquifer. Due to its acceptable ability in predicting and determining aquifer recharge, this model has a reliable capability to assess future impacts of climate change on aquifers. Examples of MODFLOW applications in groundwater flow and its response to climate change have been reported by Scibek & Allen (2006) in southwestern British Columbia, Canada, Zhang & Hiscock (2010) in Nottinghamshire, UK, as well as Wang et al. (2008) in Yangling District, Shaanxi Province in China, and Ali et al. (2012) in Australia.

The Yazd-Ardakan aquifer is the main source of water supply for Yazd province, Iran. This region is the major industrial region in central Iran (Malekinezhad 2009). The present study was carried out to determine long-term sustainable groundwater availability under current and future climate change and human pressures. Compared to previous studies, this study investigates these problems from some new perspectives, mainly: (i) using a numerical model along with a regression model to model changes in groundwater recharge; and (ii) modeling the effect of Haloxylon vegetation in the region on groundwater recharge. The output will provide a comprehensive view of groundwater systems and deterioration under natural and artificial changes. The results of this study provide a robust and essential tool for sustainable groundwater management in the study area.

Different steps applied in this paper to simulate aquifer and find groundwater level fluctuations under climate change are presented in Figure 1.

The study area is located in an arid climate in Iran with average annual precipitation of 55 mm year⁻¹. Groundwater is the main source of water supply in the study area. Effective precipitation in the mountain region of Yazd province is the main source of the groundwater recharge. Population growth along with industrialization of the area has led to overextraction of groundwater. Groundwater overextraction of about 71 Mm³ year⁻¹ in the study area has caused the groundwater level to decline by an average of 0.58 m year⁻¹ (Yazd Regional Water Organization 2014). Water stresses in recent years have increased the need for water imports from other provinces such as Isfahan and Chahar Mahal Bakhtiari. However, increasing water requirement in the water exporting regions has caused uncertainty about the future of such projects. Since groundwater is the main water supply resource in the study area, a comprehensive assessment of it under future changes is fundamental for the sustainable management of aquifers.

MODFLOW can simulate external flow stresses such as wells, areal recharge, evapotranspiration, drains, and rivers with a set of stress packages (Harbaugh et al. 2000; Harbaugh 2005).

To simulate groundwater by MODFLOW, the flow region must be subdivided into blocks in which the average characteristics are generally assumed to be homogeneous. The grid model of Yazd-Ardakan aquifer is located between 31° 35′ and 32° 30′ north and 53° 55′ and 54° 30′ east (Figure 2). It covers an area of 2,140 km². The grid model of the study aquifer consists of 100 rows and 100 columns where there are 2,583 active cells. The hydraulic gradient of the water table is in the south–north direction.

To simulate three-dimensional unsteady groundwater in this study, adequate information and a proper description of the groundwater system are required. This database includes hydraulic parameters (hydraulic conductivity, specific yield, transmissivity, etc.), budget components (recharge from rivers, precipitation, well pumping, evapotranspiration, return water of irrigation, etc.), geological data (digital elevation model (DEM) layer, aquifer thickness, etc.), and boundary conditions (impermeable boundaries and general or constant head). To assemble this information, some stages must be pursued. First, a spatial database needs to be developed, which, whatever the domain of application, permits gathering all the requisite data in one structure, making achievement easier and offering facility updates when needed (Bonomi 2009; Chesnaux et al. 2011). The characterization and the conceptualization of the Yazd-Ardakan aquifer were done using ArcGIS 10.2 software. Then, these layers were imported into the GMS 8.3 software to develop the conceptual model using multiple layers. The conceptual model was converted to MODFLOW model using a 3D grid and GMS software (Figure 3).

Basic data and information were collected and imported in GMS 8.3 to develop the groundwater model. This information is presented in different sections.

Aquifer geometry: The assessment of the geophysical and geo-electrical documents resulted in the iso-thickness map of the aquifer. Investigation of this map reveals that aquifer thickness varies from 25 to 150 m depending on the tectonic status of the region (Figure 4). Subtraction of the thickness and top elevation maps led to bed elevation layer.

Hydrological parameters: Hydraulic conductivity of Yazd-Ardakan aquifer has been extracted from available documents and studies available at the Yazd Regional Water Organization. Based on these studies, horizontal hydraulic conductivity ranges between 10 and 18 m/day (Yazd Regional Water Organization 2014). The special distribution of K_h is presented in Figure 5. This range variability is explained by the nature of sediments and their variations from silt to coarse sand. The vertical component of the hydraulic conductivity was obtained using trial and error processes as a portion of horizontal conductivity. At the end of this process, the ratio of K_v/K_h was equal to 0.8.

Water withdrawal information: Based on statistical sampling in 2009 and the last sampling in 2014, there are 945 wells in the Yazd-Ardakan basin, 821 of which with an annual extraction of 275 Mm³ are located in Yazd-Ardakan aquifer and others being located in subsidiary and small mountainous aquifers. Moreover, there are 22 deep wells with an annual extraction of 15 Mm³ for supplying drinking water. Also, 52 Mm³ of the groundwater of Yazd-Ardakan basin is drained through qanats, but only 7 Mm³ of this extraction is related to the main aquifer (Yazd Regional Water Organization 2014). These data were up to date and used in modeling. Figure 6 shows the location of water resources in the study area.

Piezometeric wells: There are more than 100 piezometric wells in Yazd-Ardakan basin, but only 59 of them are applicable for groundwater modeling. Others are unusable because of some problems including drying, statistical deficiency, or location outside of the aquifer boundaries. Data of the 59 observation wells were applied to simulate the steady state groundwater model (Figure 7). These data are useful for model calibration and validation as well.

Boundary conditions: To build the numerical model, different boundary conditions must be defined. In the case of Yazd-Ardakan aquifer, fixed-head boundary conditions in specific head and general head were applied to describe boundary conditions. Boundary conditions in the study area were simulated using general head. In some cases, groundwater recharge from the mountainous area to the plain area was simulated by well package.

Evapotranspiration from the aquifer by Haloxylon aphyllum: Haloxylon aphyllum is one of the main shrubs used in a desertification control project in Iran. There is a large area of about 424 km² of Haloxylon plantation in Yazd-Ardakan basin.

Some studies have indicated that the Haloxylon plantation would utilize groundwater. For example, Rad et al. (2011) declared that Haloxylon aphyllum water requirement is about 2.4 m³ per year. According to average and effective precipitation of the study area, which is 55 mm and 26 mm, respectively, they concluded that suitable density for Haloxylon aphyllum in the study area is approximately 112 shrubs per hectare. If density exceeded this limit, it would result in groundwater consumption.

Google Earth image survey: To determine the effect of the Haloxylon plantation on groundwater, Haloxylon density was evaluated in different regions of the study area. In order to perform this, Haloxylon vegetation boundaries were determined using Google Earth and field reconnaissance (Figure 8).

Then, using optical supervision of Google Earth images by JMicroVision software and field survey, the farmland area was categorized based on plantation densities. Some preprocessing stages were done to improve these farmland Google Earth images before using them in JMicroVision software.

JMicroVision software was used to enhance optical vision and obtain better classification (Figure 9).

To measure Haloxylon farmland evapotranspiration in the study aquifer, it was categorized into different groups according to the shrub densities using Google Earth images, JMicroVision software, and field reconnaissance. The results of this stage are summarized in Table 1.

Based on the findings presented in Table 1 and regarding the suitable density for Haloxylon in the study area (about 112 shrubs per hectare), excessive density of Haloxylon aphyllum was calculated (column 5, Table 1). Considering Haloxylon aphyllum water requirement of about 2.4 m³ year⁻¹ (Rad et al. 2011) and multiplying this number in excessive density and cultivated area, it is concluded that groundwater fed transpiration from the Haloxylon farmland is around 2.4 Mm³ year⁻¹. Recent studies showed that Haloxylon has the ability to use soil water vapor contents and so it can accelerate groundwater evaporation and depletion even in a deep aquifer (Liu et al. 2002; Rad et al. 2011; Zhu & Jia 2011). Average precipitation of the study area is about 55 mm year⁻¹ and water depth of the Yazd-Ardakan aquifer is more than 70 m. In such a situation, Haloxylon water extraction is not from soil water content, it is from soil water vapor. Using EVAP package to estimate evapotranspiration from such a depth increases uncertainty of the modeling. Therefore, the estimated extraction volume related to Haloxylon was represented in MODFLOW model as a pumping well.

MODFLOW model and Haloxylon water extraction: Water extractions of Haloxylon farmland from the aquifer were calculated based on the results of Haloxylon densities compared to the suitable density of Haloxylon farmland in the area. To consider this water removal in the MODFLOW model, the results were considered in the form of a hypothetical pumping well. Two steady state and transient MODFLOW models were run and calibrated in GMS software.

Using PEST package and trial and error process, calibration of the recharge rate and hydraulic conductivity of the steady state model was done. PEST package includes an extremely powerful nonlinear predictive analyzer. This allows the user to minimize changes between estimated and observed values of hydraulic parameters considering model constraints. Model calibration is stopped at the end of the simulation when reasonable matches between the observed and calculated hydraulic head are achieved. After each run, differences between calculated and observed heads were estimated with the aim of minimizing differences. The residual between observed and calculated heads was used to measure the root mean squared error (RMS) and the mean absolute error (MAE) presented in Table 2.

Based on results presented in Table 2, either PEST or trial and error process are reliable in the steady state model calibration. Results of the PEST process were used to calibrate the steady state model.

The transient modeling period (2005–2011) was based on the available monitoring data. The simulation period was divided into 60 monthly steps. In these time steps, groundwater level changes in response to recharge quantity fluctuations and water pumping variations.

The hydraulic conductivities and the geometry of the aquifer for the transient model are the same as those used for the steady state model. Hydraulic head changes in response to recharge, and porosity variations with regard to the time series data of the water level were used for transient model calibration. PEST package of the GMS software and trial and error process were applied to calibrate recharge rate and porosity yield of the aquifer. In order to minimize differences after each run, the residual between calculated and observed water level was estimated. The residual between observed and calculated heads was used to measure the RMS and the MAE. These results are presented in Table 3.

It is obvious from Table 3 that transient model calibration using trial and error process and PEST process concluded in approximately the same results, but the results driven from the PEST process are more reliable.

Steady state model validation was analyzed using observed and estimated data from October 2012. Using groundwater data from 1998 to 2004, transient model validation was analyzed too.

Future climate of the Yazd-Ardakan plain was predicted by statistical downscaling outputs from GCMs (HADCM₃ for SRES A₂ scenario).

The LARS-WG model is a stochastic weather generator which can be applied to simulate weather data at a single site and predict climate change impacts (Semenov et al. 1998). The climatic parameters used in the LARS-WG model include precipitation (mm), maximum and minimum temperature (°C), and solar radiation (MJ m⁻²d⁻¹). In LARS-WG, during precipitation downscaling, precipitation of each month for a local station is observed and historical data are investigated to extract statistical information such as number of dry and wet days as well as mean daily precipitation in each month of a year. This knowledge is used to construct semi-empirical distributions for the lengths of wet and dry day series and daily precipitation amount (Semenov et al. 1998).

To investigate climate change in the recharge and the depletion area of Yazd-Ardakan aquifer, two climatologic stations were selected: one in the recharge area to predict future recharge and the other in the depletion area to predict future water consumption. Historical daily maximum and minimum air temperature (°C), precipitation (mm), and solar radiation (MJ m⁻² d⁻¹) were obtained from these climatologic stations. Table 4 details the characteristics of each climatologic station.

Climate change scenario used in the present study was based on the output from HadCM₃ included in the LARS-WG model. The HadCM₃ was used to predict climate change for the IPCC 3rd and 4th Assessment Reports, and has been widely used in other studies for impact assessment. To simulate future climate the IPCC SRES (Special Report on Emission Scenarios) A₂ scenario was selected, which assumes that the current socioeconomic situation will continue (Solomon 2007). Simulated climatic data for this study were calculated for the period of 2011–2030.

To evaluate the vegetation factor, detailed knowledge of the cultivated crops during the vegetation season under identical environmental conditions is required. This information was obtained through the Food and Agricultural Organization (FAO) Irrigation and Drainage Paper 56 guidelines for arid regions. Essential information and procedures (Equations (3) and (4)) were applied in CROPWAT 8.0 to determine actual evapotranspiration.

For the 1976–2010 period, crop water requirement (CROPWAT 8.0) software was applied to calculate reference ET_o using climatic parameters obtained from Yazd climatologic station. The climate/ET_o module of CROPWAT software includes calculations, producing radiation and ET_o data using the FAO Penman–Montieth approach for the following parameters: rainfall, minimum and maximum temperature, humidity, wind speed, sunshine hours. Calculations were done using the mean value of each parameter. The module is primarily for data input, requiring information on the meteorological station (country, name, altitude, latitude and longitude) together with climatic data.

Irrigated agriculture is considered the main groundwater extraction sector in the study aquifer. The share of agriculture in the aquifer's water extraction exceeds 90% (Yazd Regional Water Organization 2014). To predict human pressure on groundwater resources under climate change, it is essential to estimate water requirements of the cultivated area under future climate variations. The current agricultural water requirements (differences between crop evapotranspiration and precipitation) were calculated using CROPWAT 8.0 software. Future simulated climatic parameters were applied to estimate evapotranspiration and water requirements for the dominant cultivated crops under climate change.

The calibrated and validated groundwater model could be applied for different management and planning studies. To predict the aquifer response to future changes, the optimized parameters obtained through model calibration are useful. The simulations were performed to clarify the aquifer variations in response to climate change (recharge rate) and human pressures (pumping rate) on the aquifer. In order to have sustainable management of the aquifer in the study area, different scenarios mentioned in Table 5, were applied.

Different steps were applied in simulating climate change, Haloxylon farmland and groundwater. Results of these stages are presented in different sections.

Results of the groundwater modeling is presented in different stages including simulation, calibration, and validation.

Simulated groundwater level distribution is presented in Figure 10.

As is shown in Figure 10, groundwater level was accurately simulated by the MODFLOW model. Almost all prediction error is located in the reliable range.

In the calibration stage it is important to obtain an optimal fit between measured and calculated data. Also, the reliability of the operational model depends on this stage. In this study, two consecutive steps including model calibration and validation were used to develop the MODFLOW model of the aquifer.

Steady state model: Using PEST package, the steady state model was calibrated. Figure 11 shows a scatter diagram of the observed and calculated hydraulic head resulting from the steady state model calibration.

Results of the steady state validation stage are shown in Figure 12.

Based on Figure 12, there is a good fitness between the estimated and calculated groundwater level of piezometers. It is concluded that the steady state model is reliable.

Transient model: The comparison of the calculated and observed water level related to the transient model calibration is shown in Figure 13. This output confirms appropriate model calibration.

A systematic overview of all simulated versus observed groundwater level graphs (Figures 11 and 13) indicated that a reasonable level of model-to-measurement misfit is achieved. Comparison of misfit values of the steady state and the transient model shows further increase of the groundwater level, which is due to obtaining the parameters of the model through steady-state calibrations and adopting them in the transient calibration. Hydraulic parameters value obtained from the calibration stage including K_h, K_v, and S_s are considered to be reasonable compared to available field values of these parameters. These outcomes confirm MODFLOW's ability in modeling Yazd-Ardakan aquifer.

The results of transient model validation are presented in Figure 14. This figure shows unit hydrograph of the Yazd-Ardakan aquifer.

Based on validation results of the transient model presented in Figure 14 for a seven-year period from 1998 to 2004, there is good agreement between estimated and observed data.

By using the A₂ emissions scenario of the LARS-WG model, precipitation was predicted for the period from 2010 to 2030 using the GCM of HADCM3 (Figure 15).

Based on Figure 15, the precipitation distribution will be substantially changed in the future. This pattern change will result in precipitation increasing in spring and decreasing in autumn and winter, especially in the recharge area (Dehbala station). The other important factor of climate is the temperature parameter. The absolute changes in the surface temperature (°C) of two stations (recharge area and water extraction area) in the years from 2010 to 2030 at the A₂ emission scenario were calculated by the HADCM₃ model and are shown in Figure 16.

According to Figure 16, the maximum temperature will increase in all of the months for the two stations until 2030. The changes for the minimum temperature of the two stations are shown in Figure 17 using the GCM model of HADCM₃ and emission scenario A₂.

According to Figure 17, similar to the maximum temperature, the minimum temperature will increase in all of the months for the two stations until 2030.

From the outcome of Figures 15–17, climate change is going to happen in the study area; the study area will have less rain and hotter summers. This finding is equal to that of the studies of Green et al. (2011) and Seiler et al. (2008). These researchers also indicated that climate change leads to an increase in temperature and variations in the amount and frequency of precipitation.

To model aquifer response to climate change, statistical data including recharge and AARR for the recharge area (Dehbala station), were divided into two parts: 70% of data was employed as training the data and the remaining 30% for testing. Using recharge and AARR data along with different lag between rainfall and groundwater recharge, the best model was recognized through MATLAB 2014 software. Figure 18 represents the selected regression model. This model was processed between AARR data and recharge data with a 35-month delay.

Selected regression model validation was done using 30% of the historical data. Figure 19 demonstrates the validation results.

As is obvious from Figure 19, the rainfall-recharge model is confirmed as an appropriate model that can predict groundwater response to future climate change.

The predicted recharge rate was compared to the present rate in Figure 20.

Based on the results presented in Figure 20, aquifer recharge will increase in winter and autumn despite decreasing in summer. According to Figure 15, precipitation variations in recharge area increase in spring and decrease in winter. Differences between Figures 15 and 20 are interpreted by the regression model of recharge and precipitation. As mentioned before, there is a 35-month delay between AARR data and recharge data. Hence, these trends of parameters differ during various months. Based on the results of the current study, recharge rate will be increased to 3.5% annually. From previous studies it was concluded that in arid and semi-arid regions, increase in amount of rainfall events results in groundwater recharge reinforcement (Bates et al. 2008). This conclusion, which was explained by future climate change in the study area, resulted in rainfall distribution changes and increasing heavy rainfall events.

Climate change prediction for Yazd station was applied as basic data in CROPWAT 8.0 software to estimate evapotranspiration and water requirements for the dominant cultivated crops under climate change. These results can be considered as excellent representing human pressure on the aquifer. Figure 21 represents human pressure on the aquifer under climate change in the period of 2010 to 2030 with respect to the observed period based on A₂ emission scenario of the HADCM₃ model.

According to Figure 21, agricultural water requirement has an increasing trend in all months, especially in the spring and summer by the year 2030. Based on the results of this study, human pressure on the aquifer will increase up to 7% compared to the present time.

Groundwater recharge results from effective precipitation (that is precipitation minus evapotranspiration) which infiltrates into the subsurface from where hydraulic gradients are downward (Taylor et al. 2012). In Yazd-Ardakan aquifer, the groundwater recharge area is mainly located in a mountainous region beyond the aquifer boundary. This recharge flows to the main aquifer based on hydraulic gradient. Although changing in the future, evapotranspiration will lead to increasing agricultural water demand; the proportion of agricultural water returned to the aquifer from the irrigated region is negligible compared to the recharge resulting from the mountainous area. Thus, groundwater recharge will increase as a result of future climate change.

In order to fulfill the objectives of this work, the comparison of different scenario analyses of the groundwater management is represented. Table 6 summarizes the comparison of the results of the impact of the climate change and the impact of different groundwater management scenarios. In order to facilitate the comparison, combinations of different situations are presented grouped in the three basic categories appearing in column two of Table 6.

To evaluate groundwater response to climate change and human pressures, three combinations of climate change and pumping schemes were considered (Table 6). Scenario 1 includes the climate change imposed variations in recharge, and assumes pumping rates to increase by 7%. The groundwater flow modeling shows that in this scenario, the aquifer storage will decrease by 67 Mm³ year⁻¹ and average groundwater level will decline by 0.6 m year⁻¹. Scenario 2 assumes a reduction in the pumping rates of 5% which results in an aquifer storage decline of 38 Mm³ year⁻¹ and an average decline in groundwater level of 0.3 m year⁻¹. In Scenario 3, the current pumping rate will decrease by 25% to conserve and restore the groundwater resources in Yazd-Ardakan aquifer. In Scenario 3, the reduced pumping rate by 25% results from irrigation system modification and cropping pattern improvement. In this scenario, the average hydraulic head and annual storage of the aquifer will increase to 0.1 m year⁻¹ and +10 Mm³ year⁻¹, respectively. The results of these simulations indicated that groundwater conservation and aquifer restoration is possible through proper management of water resources in the study area.

Numerous studies show that global warming and climate change are occurring throughout the world. To deal with these changes, it is necessary to consider the effects of future changes on water resources. Investigation of climate change on groundwater resources in Yazd-Ardakan aquifer showed that future climate variation in the study area will result in increasing trends in both aquifer recharge and water extraction from the aquifer. To relieve undesirable effects of climate variation on the aquifer, three different scenarios were simulated using the MODFLOW model. The results of these simulations revealed that groundwater conservation and recovery is possible under irrigation system modification and cropping pattern improvement (using drought-resistant plants). Also, there is some debate about Haloxylon forests' groundwater use in the region. Recent studies showed that Haloxylon aphyllum has the ability to use soil water vapor contents and so it can accelerate groundwater evaporation and depletion (Liu et al. 2002; Zhu & Jia 2011; Rad et al. 2011). Looking at this concept, the amount of Haloxylon groundwater extraction was calculated and simulated in the MODFLOW model. It was the first time that this mechanism was simulated in modeling a deep aquifer. The results of this study showed that groundwater extraction of Haloxylon forests is negligible compared to human water depletion. Considering Haloxylon forests' advantages, especially in wind erosion control, it could be concluded that if their density was reduced to an acceptable rate, Haloxylon forests are beneficial in the study area.