Groundwater budget deficit caused by drought and overexploitation

Detecting the relative contributions of climate and anthropogenic impacts on hydrological systems is a necessity for the management of water resources (Knowling et al. 2015). Drought, a result of climate change, and overexploitation, with human causes, are the main factors affecting groundwater depletion. Several studies have either separately or simultaneously analyzed these two terms (Taylor et al. 2009; Sheffield & Wood 2012). Recently, understanding the relative contributions of stressors to hydrological systems has been the subject of intense research (e.g. Hao et al. 2008; Dai et al. 2010; Li et al. 2014; Knowling et al. 2015; Asoka et al. 2017).

Several research groups have made efforts to separate drought impacts from human influences on groundwater systems. For instance, Van Loon & Van Lanen (2013) proposed an observation-modeling framework to quantify the effects of natural (drought) and human (overexploitation) impacts on the anomalies of groundwater. Knowling et al. (2015) focused on the relative contributions of climate and human stresses to aquifer depletion using the MODFLOW model. Some studies have employed a combination of the models (seasonal version of WetSpass and MODFLOW) to estimate groundwater recharge via water balance (Batelaan & De Smedt 2007; Ghouili et al. 2017). Based on these two models, there has been some research on the effects of climate scenarios on the level and recharge of groundwater (Woldeamlak et al. 2007; Ajjur & Mogheir 2012) and some studies have estimated the amount of pumped groundwater required for irrigation (El Idrysy & De Smedt 2006). They concluded that combining these models better contributes to groundwater management on both seasonal and annual scales. Mustafa et al. (2017) employed WetSpass and MODFLOW to analyze the role of different factors influencing the depletion of groundwater levels in northwestern Bangladesh. They showed that the groundwater level deficit was significant and had a minimal correlation with the Standardized Precipitation Index (SPI), which they ascribed to the increasing trend in groundwater abstraction. Previous studies have utilized different methods such as statistical analysis and model-based investigations, such as Lumped, and one- to three-dimensional models to estimate groundwater recharge and identify groundwater stressors. In the present research, a spatially distributed monthly water balance model (WetSpass-M) was combined with a groundwater model (MODFLOW) to separate the groundwater budget deficit caused by drought and overexploitation in Aleshtar aquifer, a plain located in Lorestan province, western Iran. Another novelty of this study is the quantification of the groundwater budget deficit by stressors in three clustered locations of the aquifer. These clusters were classified based on similar behaviors of groundwater level measured by hierarchical K-means algorithm and hydraulic diffusivity index. Overall, the main contribution of the current study is its use of water balance as an indicator to separate the natural impacts of groundwater budget deficit from man-made impacts or overexploitation. The ultimate objective is to enable managers to plan their risk management tools for mitigating or controlling the drought and overexploitation, separately.

The study area covers the watershed and aquifer of Aleshtar in Lorestan province, west of Iran (Figure 1). The watershed area is approximately 798 km², of which about 110 km² is covered by aquifer. According to the latest geologic maps from the Geological Survey and Mineral Explorations of Iran (GSI) and National Iranian Oil Company (NIOC), Aleshtar watershed includes a range of limestone mountains with high permeability rocks in the north and northeastern part and conglomeratic reliefs with low permeability in the west. A combination of different materials with very low permeability such as sandstone, marl, and clay further covers the southern part. The average annual precipitation in the watershed and aquifer is estimated at approximately 533 mm and 512 mm, respectively. The annual average temperature and total pan evaporation in Aleshtar plain are around 12.9 °C and 2,238 mm, respectively.

According to the hydrogeological studies (Zhrfab Payesh 2005), the type of aquifer is unconfined. Pumping tests and geo-electrical investigations have shown that the transmissivity of aquifer varies between 897 (m²/day) in the west and 1,650 (m²/day) in the eastern north-southern strip. Storage coefficients range from 0.032 in the south to 0.051 in the northern and eastern parts of the aquifer. Data for 176 exploitation wells and 15 observation wells are fully available on a monthly scale from the regional water company Lorestan (http://www.lsrw.ir/). Figure 1 shows the lithological map of the region and the distribution of exploitation and observation wells. The authors prepared this map based on the geologic maps from GSI, NIOC, and data related to Lorestan Regional Water Company.

In this research, a two-step methodology was applied: monthly groundwater recharge was simulated using the WetSpass-M model, and groundwater storage changes were simulated using recharge and estimated through combining the WetSpass-M and MODFLOW models. To quantify the relative contributions to groundwater budget deficit, the groundwater model was run under impacts of drought and overexploitation. Drought anomalies were investigated using a threshold level method and standardized groundwater-level index (SGI) (Soleimani Motlagh et al. 2017). This stage helps estimate groundwater budget deficit caused by drought. The following sections provide a detailed description of the steps in the present study.

Monthly groundwater recharge was simulated using WetSpass-M model through water balance, developed by Abdollahi et al. (2017).

The WetSpass-M model utilizes four static layers, namely soil, digital elevation model (DEM), slope, and land use, and monthly layers of seven dynamic hydro-physical variables. In this research, data for precipitation, temperature, pan evaporation, wind speed, snow cover, groundwater depth, and leaf area index (LAI) were obtained over a period of 276 months from October 1991 to September 2014. The grid cell resolution in this study was 50 m × 50 m with 829 columns and 596 rows. In addition to the above layers, the number of rainy days and snowmelt degree-day factors per month were introduced to WetSpass-M model.

The soil map and DEM layer with 50-m resolution were obtained from the Soil Conservation and Watershed Management Research Institute of Iran. Topography was in the elevation range of 1,486–3,613 m and 0–218% slope. Based on the Food and Agriculture Organization classification, the soil of the study area was classified into seven classes of soil texture, namely loamy sand, sandy clay loam, silty clay loam, clay loam, sandy clay, silty clay, and clay.

The land-use map was prepared based on the Landsat-5 TM and Landsat-8 OLI images for May and June of 1991, 2000, and 2013, which were downloaded from http://earthexplorer.usgs.gov. The radiometric and atmospheric corrections were done by Fast Line-of-sight Atmospheric Analysis of Hypercubes (FlAASH) tool. In addition, the old maps of land use and the field observations in control points were utilized to update the new land-use/cover maps. Finally, the maximum likelihood algorithm method in ENVI software, with a high accuracy reported in other studies (Al-Ahmadi & Hames 2009), was employed to classify the land use. The accuracy of this method, which employed kappa coefficients in contrast to control points, was over 80% in different images.

Using geostatistical tools (Kriging and Inverse Distance Weighting), monthly dynamic hydro-physical variables such as precipitation, temperature, pan evaporation, wind speed, and groundwater depth were converted into a grid. Data of these variables, belonging to the synoptic and climatology stations inside and outside the study area, were collected from Meteorological and Water Organizations of Lorestan Province, Iran. In addition, long-term monthly LAI images were taken using a moderate resolution imaging spectroradiometer (MODIS, http://earthexplorer.usgs.gov) and an advanced very-high-resolution radiometer (AVHRR, ftp://ftp.glcf.umd.edu/glcf/GLASS/LAI/AVHRR/) products; these images were then resampled using the nearest technique based on the grid cell dimensions (50 m). In addition, long-term monthly snow cover images have been provided by MODIS products since 2000. As these images had not existed before 2000, they were estimated based on the correlation equations among snow cover and precipitation and temperature variables during the snowy months (Table 1).

Calibration (1991–2008) and validation (2009–2014) of WetSpass-M model were carried out seperately using river flow data at the outlet of the study area (Sarab-Seyedali station). Calibration was manually done via trial and error through changing the model's input parameters including ∅︀, β, x, a, surface runoff parameter (LP), α, and snowmelt factor (MF) for recharge contribution, storage, runoff delay, interception, surface runoff, evapotranspiration, and snowmelt, respectively. Sensitivity analysis using a method facilitating model calibration was performed on three sensitive parameters (a, LP, and α) (Abdollahi et al. 2017). Simulated data during calibration and validation periods were compared with the observed data (direct runoff and base-flow) at the outlet of the basin. To separate the stream-flow hydrograph, a non-tracer base-flow separation (filtering and recursive methods) called Web-based Hydrograph Analysis Tool (WHAT) system was used. Eckhardt (2005) and Gonzales et al. (2009) corroborated the good performance of these methods in water balance assessment and calibration/validation of hydrological models. In the present research, the optimal method of river flow separation was selected based on the extreme values of base-flow and direct runoff, such that in dry seasons, the direct runoff was close to zero and the river flow came close to the base-flow.

WetSpass-M model was assessed by Nash–Sutcliffe Efficiency (NSE) coefficient (Nash & Sutcliffe 1970), R², and root mean square error (RMSE) along with Mann–Whitney test for a better assessment.

In the conceptual model, the hydrogeological features of the groundwater were identified using all available data for running the mathematical model. Based on alluvial deposits, pumping test analysis, geo-electric surveys, and drilling logs provided by consulting companies, the aquifer is identified as unconfined (Zharfab Payesh 2005). In the MODFLOW model, this area was vertically divided into one layer, discretized into a grid of 333 × 237 cells with a spacing of 50 m. The model grid of Aleshtar aquifer consists of 43,879 active cells. General groundwater in this region flows from north and northeast to south. Accordingly, three types of boundaries were considered in the aquifer margin based on the maps of the geologic and groundwater contours. The northern and northeastern boundaries of the aquifer were treated as lateral inflow boundaries, where the area is surrounded by permeable geological formations mainly containing dolomites and limestones, defined as the active cells of the groundwater received in the general head boundary (GHB) package. The initial values of these boundaries were calculated according to Darcy's Law and further calibrated in the model. The sources of intra-model lateral feeder flows belong to calcareous mountains in the north and northeastern parts for Cluster 1 and the north part for Cluster 3. As far as Cluster 2 is concerned, these sources mainly belong to the downstream part of Cluster 1 and a small part in the downstream of Cluster 3.

A small part to the south of the aquifer, including alluvial materials, was considered as the outlet boundary. These cells take groundwater out from the aquifer. The remaining borders of the aquifer were identified as no-flow boundaries due to impermeable or low permeable deposits containing Pliocene Bakhtiari conglomerate, Eocene conglomerate and sandstone, and Eocene marly limestone, in the west, northwestern, and southeastern of the aquifer, respectively (Figure 1). As inactive cells in the GHB, these boundaries are parallel to flow direction, through which no water flows into the model domain. Based on Darcy's Law, the initial values of the inflow and outflow were manually estimated at about 20 and 4 million cubic meters (MCM), respectively, in 1991. In November 1991, the hydraulic heads were selected as initial conditions in the steady state due to fewer fluctuations. The top layer of the aquifer was obtained by DEM at a cell size of 50 m. The bottom layer was made via interpolating the depths of the four exploratory wells and exploitation wells drilled close to the bedrock of the aquifer with a spacing of 50 m. The thickness of the aquifer was calculated between 38 and 111 m in the northwest and south of the aquifer, respectively. Storage coefficient was obtained by interpolating the results of pumping test in four exploratory wells for all aquifer networks. The map of the hydraulic conductivity was created based on the transmissivity amounts and the thickness of the aquifer.

In the groundwater model of the studied area, the initial value of the riverbed hydraulic conductance was defined as varying from 32 to 498 m²/day and then optimized during calibration. The amount of initial evapotranspiration from groundwater level was also obtained using a White graph (White 1932) about 1.13 MCM and then calibrated in the model. This graph shows a linear relationship between the evaporation percentage of the pan and the depth of groundwater table. The pumping of the wells was simulated using monthly pumping rates of 176 exploitation wells in Aleshtar plain, which was about 24 MCM in 1991. In addition, the estimated recharge based on water balance by WetSpass-M model was used in the same grid with groundwater model for MODFLOW simulation. Based on Zhrfab Payesh (2005), initial infiltration values via the return-flow of the irrigation (about 15% of the total consumed water), industry, urban, and domestic water (about 80% of total consumed water) were defined for the model in the initial year of the simulation and then calibrated.

The groundwater levels of the 15 observation wells in November 1991 were used to calibrate the model under a steady state. Moreover, the data from November 1992 to October 2007 were used for the transient calibration. The model was validated with the available data from November 2007 to October 2014. In this model, 30 and 14 stress periods were defined for the transient calibration and validation, respectively. The stress periods were assigned based on the periods of recharging, pumping, and surveying the trends of monthly hydrograph per hydrological year. The performance of the model simulation was assessed with R², NSE, and RMSE.

Natural impacts (climate-oriented or drought) and human impacts (pumping or overexploitation) on groundwater budget deficit were characterized by a spatially distributed model based on MODFLOW for each hydrogeological zone in the studied aquifer. To categorize the hydrogeological zones in the studied aquifer, the approximate range of each zone was determined using the hierarchical classification method of K-means. Further details of this method were explained by Soleimani Motlagh et al. (2017). Afterwards, the boundary of each zone was optimized based on hydraulic diffusivity (transmissivity/storage coefficient) (Figure 1) because it plays a significant role in assessing the importance of the transience reaction to various types of stress (Barlow & Leake 2012).

Since the objective was to determine the annual groundwater budget deficit caused by stresses, their relative contribution to groundwater was calculated by automatic modeling. In this regard, the calibrated groundwater model was used in the absence of pumping to quantify the impacts of climatic stresses, assuming that natural stresses affected the groundwater deficit alone and the groundwater pumping had no roles. The relative contribution of overexploitation to groundwater by wells was obtained by subtracting the computed deficit from annual groundwater storage under both stresses (climatic and human stresses). Negative annual storage values of the groundwater are the budget deficit induced by stressors.

Figure 2 presents the separation results of base-flow and direct runoff in the outlet of the basin using one parameter digital filter on a daily basis. According to Figure 2, both daily flow hydrograph and its components (direct runoff and base-flow) followed daily precipitation. According to the results, about 85% of the river flow was in the base-flow form. Subsequently, the daily values were converted to monthly mean.

The results of the WetSpass-M model sensitivity analysis showed that LP was the most sensitive parameter (followed by α and a); therefore, by increasing the amount of LP, the runoff decreased and the recharge increased. The order of parameter sensitivity in the present watershed is different to that of Abdollahi et al. (2017), possibly due to the local condition of the study area. The optimum results of the WetSpass-M parameters are shown in Table 2.

Results of the WetSpass-M model calibration (for 75% of the total data from October 1991 to December 2008) and validation (for 25% of the total data from January 2009 to September 2014) showed that the simulated data (runoff and base-flow) are in agreement with the filtered data by WHAT (Figure 3(a) and 3(b)). The evaluation of the model further corroborates these results (Table 3).

The simulated recharge maps associated with each month per year were summed up to compute the annual recharge. Figure 4 presents the annual potential groundwater recharge for the hydrological normal (1993–1994) and dry (1998–1999, 1999–2000, and 2007–2008) years. As observed, with the increase in the severity of droughts, the groundwater recharge rate was reduced.

The hydraulic conductivity of the aquifer was optimized in the range of 3–9 m/day during steady-state calibration. The optimum values concerning the hydraulic conductance of the Aleshtar rivers ranged from 52 to 592 m²/day. Moreover, the hydraulic conductance of the permeable lateral aquifer boundaries was optimized in a range of 200–508 m²/day.

Figure 5(a) shows the scatter plot of the values simulated versus observed in a steady state with an R² of 0.99. Figure 5(b) shows the volumetric water balance for the steady state. Figure 5(c) presents the results of simulations for 15 observation wells in the transient period with an R² of 0.9, NSE of 0.89, and RMSE of 0.72 m.

As shown in Figure 6, the inflow components of the groundwater, including recharge by the river, lateral inflow, and recharge by return-flow had a rising trend (the slope of linear equations is positive), but precipitation and recharge from precipitation did not. According to this figure, at the beginning of the simulated period, the inflow components of the balance almost corresponded to precipitation changes up to the hydrological year of 1998–1999. Subsequently, the recharge through the river and the aquifer boundary (lateral inflow), despite their increasing trends, did not completely follow the precipitation changes. The amount of return-flow entering the aquifer was low at the beginning of the study period, increasing from 1998–1999. In this figure, all outflow components of the groundwater budget had a falling trend, except for pumping from wells (consisting of the river discharge from the aquifer, the lateral outflow from the aquifer boundary, and the evapotranspiration from groundwater level). Such changes in the groundwater balance components result from groundwater overexploitation and drought (for instance, precipitation shortage in 1998–1999, 1999–2000, 2007–2008, 2008–2009), shown in Figure 6(c) and 6(d).

More details can be observed in the graphs grouped in Figure 6. As shown in Figure 6(a), at the starting point of simulation (1998–1999), because of low pumping and high groundwater levels (groundwater table is above the riverbed surface), values of the river drainage was more than its recharge; however, after that year, the recharge by the river increased more than its discharge as a result of droughts and overexploitation of wells. Therefore, the contribution of recharge and discharge by the river is related to the influence of stress factors in the same years. These changes entailed the loss of the flow from the river to the aquifer, possibly creating an imbalance in the hydrologic system of the watershed. Sanz et al. (2011) and Stefania et al. (2018) reported similar results regarding the impacts of groundwater pumping on the abduction of river flow to aquifers and reduction in river discharge. Figure 6(b) also shows that from the late 1990s onwards the lateral inflow and outflow values of the aquifer were significantly varied. These changes may be associated with the increase in the groundwater pumping and recharge shortage from precipitation. The fact that these factors have increased may be due to the change in the hydraulic gradient and increased lateral inflow to the aquifer caused by the groundwater depletion.

To better evaluate the effects of stress factors on groundwater budget deficit, the aquifer area was divided into three zones (Figure 1). Division was carried out based on the identification of similar behaviors in the water level measured in the observation wells using a hierarchical K-means clustering algorithm, the hydraulic diffusivity, the density of exploitation wells, and the volumetric flow fed or drained by the rivers. Table 4 shows the hydrogeological properties of each cluster.

According to Figure 1 and Table 4, the highest and the lowest percentages of the plain area belonged to Clusters 1 (about 45%) and 2 (around 24%), respectively. On average, Cluster 2 had the highest hydraulic diffusivity (high transmissivity and low storage coefficient). In this cluster, the average groundwater extraction was very low in comparison with the other two clusters. In addition, the rivers in this cluster drained more water from the aquifer. In Cluster 1, the hydraulic diffusivity index was very close to that of Cluster 2. The average discharge from the wells in this area was significantly high (about 30 MCM). In Cluster 3, the mean discharge from the well was lower than that in Cluster 1. Furthermore, in this cluster, the recharge and discharge values of the river were very small, and the average hydraulic diffusivity was the lowest.

Table 5 presents the relative contribution of the stress factors to the groundwater budget deficit in dry years on the entire aquifer and three clusters. According to the results of the SGI index by Soleimani Motlagh et al. (2017), the entire aquifer experienced drought at the onset of the 2000s and the late 2000s, specifically. These droughts were associated with the meteorological droughts (demonstrated by SPI) with a lag time of 12–24 months. In addition, these droughts were closely matched with the predictions of groundwater models in natural conditions when groundwater storage was lower than the mean threshold level. According to Table 5, on average, groundwater balance deficit by stressors in Cluster 1 was the highest because of the over-withdrawals of the groundwater by wells. Interestingly, despite the high extraction in this area, the hydraulic diffusivity and recharge rate of the aquifer by the river were high (Table 4), making up for the storage deficit in the aquifer. Despite the lower area of Cluster 3 compared to Cluster 1, the groundwater budget deficit was substantial, probably due to a low hydraulic diffusivity, high groundwater abstraction, and low groundwater recharge in drought years. Therefore, the compensation of its deficit needs to be done cautiously.

On average, the relative contribution of drought and overexploitation to groundwater storage deficit during the dry years was 2.44 and 3.32 MCM, respectively, meaning the share of overexploitation was about 36% more than drought. These results are in line with Van Loon & Van Lanen (2013) and Knowling et al. (2015) who concluded that due to climate impacts on hydrologic systems, the groundwater abstraction, on average, had a greater impact compared with the drought and depletion of the aquifers.

It is indispensable to quantify the relative contribution of drought and overexploitation to groundwater budget deficit for groundwater management. The present research integrated two spatially distributed models, namely WetSpass-M and MODFLOW, to provide an appropriate basis for the separation of stressors impacts on the groundwater storage deficit. Furthermore, to distinguish the role of stressors, hydrogeological specifications were analyzed, and the hydraulic diffusivity factor was used in different clusters of the aquifer. The relative contribution of drought and overexploitation to groundwater budget deficit of Aleshtar aquifer was specified by means of the presented methodology. In order to include the spatial variability of the groundwater, recharge was estimated using the WetSpass-M model and calibrated with the MODFLOW model. Results obtained from stream-flow simulations showed that WetSpass-M had a high efficiency in simulating groundwater recharge. The outputs were used to assess groundwater budget and identify the contribution of the stress factors to the groundwater budget deficit, where all inflow parameters (except for recharge) had a rising trend while all outflow parameters (except for groundwater withdrawals) had a falling trend. The causes of such anomalies can be attributed to either overexploitation or natural droughts. Assessment of the relative contribution of the causes in the studied area showed that the deficit caused by overexploitation (pumping stress) was about 36% higher than natural drought (climate stress). However, this conclusion may be incomplete if we ignore the compounding effects of intensified aquifer abstraction during the dry years. Attention to the origin of the stressor may provide a more valid basis to seek better risk management plans for mitigating or controlling the drought and overexploitation, separately.

This work, as a PhD thesis, was financially supported by the University of Kashan, Iran, for which the authors are grateful.