**Potential Impacts of Future Climate Change Scenarios on Ground Subsidence**

**Antonio-Juan Collados-Lara 1,\*, David Pulido-Velazquez 1, Rosa María Mateos <sup>2</sup> and Pablo Ezquerro <sup>3</sup>**


Received: 12 December 2019; Accepted: 8 January 2020; Published: 13 January 2020

**Abstract:** In this work, we developed a new method to assess the impact of climate change (CC) scenarios on land subsidence related to groundwater level depletion in detrital aquifers. The main goal of this work was to propose a parsimonious approach that could be applied for any case study. We also evaluated the methodology in a case study, the Vega de Granada aquifer (southern Spain). Historical subsidence rates were estimated using remote sensing techniques (differential interferometric synthetic aperture radar, DInSAR). Local CC scenarios were generated by applying a bias correction approach. An equifeasible ensemble of the generated projections from different climatic models was also proposed. A simple water balance approach was applied to assess CC impacts on lumped global drawdowns due to future potential rainfall recharge and pumping. CC impacts were propagated to drawdowns within piezometers by applying the global delta change observed with the lumped assessment. Regression models were employed to estimate the impacts of these drawdowns in terms of land subsidence, as well as to analyze the influence of the fine-grained material in the aquifer. The results showed that a more linear behavior was observed for the cases with lower percentage of fine-grained material. The mean increase of the maximum subsidence rates in the considered wells for the future horizon (2016–2045) and the Representative Concentration Pathway (RCP) scenario 8.5 was 54%. The main advantage of the proposed method is its applicability in cases with limited information. It is also appropriate for the study of wide areas to identify potential hot spots where more exhaustive analyses should be performed. The method will allow sustainable adaptation strategies in vulnerable areas during drought-critical periods to be assessed.

**Keywords:** ground subsidence; climate change; Vega de Granada aquifer

#### **1. Introduction**

In many agricultural regions, as well as areas with a rapid urbanization and population growth, prolonged groundwater exploitation due to increasing water demand is causing land subsidence impacts [1–7]. In many coastal and delta cities, such as Ho Chi Minh, Bangkok, Manila, Tokyo, and Jakarta, land subsidence considerably exceeds (by up to 10 times) absolute sea level rise, which increases flood vulnerability and triggers severe, damaging impacts [8]. Sinking cities—within the framework of global change—is and will be a transnational threat.

Land subsidence is related to falling groundwater levels in (generally) unconsolidated alluvial or basin-fill aquifers with a significant proportion of compressible fine-grained materials. Increases in effective stresses—caused by headwater declines—determine the aquifer system compaction at local or regional scales [9]. This deformation is typically elastic (reversible) and results in small vertical displacements, but when the aquifer is subjected to head declines that exceed the critical levels, much of the compaction is related to an inelastic deformation and the accompanying subsidence is permanent [10].

Differential interferometric synthetic aperture radar (DInSAR) techniques are satellite-based remote-sensing tools that have been applied successfully for monitoring land subsidence thanks to their high spatial and temporal coverage, fast data acquisition, and low cost [11]. They are capable of measuring mm-scale ground displacements at a spatial resolution of 5–10 m over large regions (hundreds to thousands of square kilometers). They have been widely applied in numerous recent studies of land subsidence triggered by intense groundwater withdrawals [12–15]. A prominent case is Jakarta, the capital city of Indonesia (around 10 million people), where subsidence values of 28 cm/year have been registered in some locations. The impacts in Jakarta have been seen in several forms: not only cracking and damage in buildings and infrastructure, but subsidence also enlarges the (tidal) flooding inundation areas and makes the coast more vulnerable to sea-level-rise phenomena [1]. Indonesian authorities think the capital should be relocated.

In the present work, we applied the persistent scatterer interferometry technique (PSInSAR), an operational tool for precise ground deformation mapping, acting as a geodetic network [16]. This technique allows quantification of deformation measurements, combining them with geological and hydrogeological data in a geographical information system (GIS). The application of PSInSAR has improved rapidly in the last decade, and it is now a valuable tool for monitoring seasonal and long-term aquifer-system responses to groundwater pumping.

In the literature, we found many different approaches that have been employed to model the impact of groundwater level depletion on the subsidence. Some of them are based on physically based distributed hydro-geomechanical models [17,18], but we also found other approaches such as machine learning [19] and conceptual [20] or simple regression approaches [21], like the one that was employed in this work.

In order to assess potential future impacts of climate change (CC) in a rational way, we need to use the climatic scenarios generated by simulating with climatic models the emission scenarios identified by the Intergovernmental Panel of Climate Change (IPCC). The most recent of these are those included in the 5th assessment report [22], the RCP scenarios. These scenarios are not predictions of future climate; they are internally consistent pictures of plausible future climates that constitute a basis for other workers to evaluate the possible impacts of CC [23]. In order to make this information relevant to assessment of impacts in specific case studies, they have to be translated to regional and local scales by applying statistical correction techniques that take into account historical data [24].

Despite the spread of uncertainties involved in the assessment of future CC impacts, there is no excuse for delays or inaction in assessing/identifying adaptation strategies, taking into account that there are environments that could be very vulnerable [25]. The market for technologies for adaptation to CC is growing rapidly, given that "the cost of repairing damages is estimated to be six times greater than adaptation costs" (H2020WATER-2014/2015). In the literature, we found multiple examples of the propagation of future CC scenarios in order to assess hydrological impacts at different scales, including continental [26,27], country [28,29], river basin [30,31], and aquifer [28,32] systems. Nevertheless, the literature on the assessment of the potential future impacts of CC on land subsidence is very limited [33]. Only a few works have been developed on identification and assessment of potential adaptation strategies, like the work published by Brouns et al. [34], in which a bottom-up approach was applied to define the scenarios.

The main objective of this research was to propose a new, parsimonious approach [35,36] with which to perform a first assessment of potential CC impacts on land subsidence related to groundwater-level depletion in detrital aquifers. We also applied and validated it in the Vega de Granada aquifer (southern Spain). We have proposed a general approach applicable to any type of aquifer including in coastal areas, which are more vulnerable to subsidence [37,38]. The main innovative aspect of the proposed methodology is its applicability in cases with limited information through a combination of several parsimonious techniques: (1) a statistical approach for CC generation, (2) a simple hydrological approach to assess the impacts of CC on groundwater levels, and (3) simple linear regression models to propagate the impacts to land subsidence. The model might be also useful for a preliminary assessment of different adaptation strategies. It does not require a distributed groundwater flow model of the system [39,40]. Taking into account the uncertainties around future potential CC scenarios, the model provides a useful first approach, even in wide areas (e.g., country or continental scale), that can identify potential hot spots where more exhaustive analyses should be performed.

#### **2. Case Study**

The Vega de Granada is a flat region located in the metropolitan area of the city of Granada (southern Spain, 530,000 people) (see Figure 1). With an extension of 200 km2, the Vega de Granada is a traditional agricultural region with increasing urban growth during recent decades. Both agricultural and urban demands exert a high pressure on the aquifer. The river Genil flows (from SE to NW) through the center of the basin, being the aquifer's main drainage axis [41].

**Figure 1.** Location of the case study. The Vega de Granada aquifer is of regional importance. It is located within the Granada Basin domain (southern Spain) and in the central sector of the Betic cordillera (Sierra Nevada).

From the geological point of view, the Vega de Granada is within the Granada Basin domain, in the central sector of the Betic cordillera. It is a tectonic depression delineated by normal faults which control the basin infill. The Vega de Granada detrital aquifer presents a multilayer and heterogeneous structure with quaternary levels of gravel, sand, silt, and clay, with a maximum thickness of 250 m in its central part [15,42,43]. Both the river sedimentation and the surroundings alluvial fans (related to the normal faults) determine the facies distribution.

The Vega de Granada aquifer is of regional importance. With renewable water resources of about 160 hm3/yr [43], groundwater exploitation has intensified considerably in recent decades because of urban sprawl, particularly during the severe droughts that periodically affect the region. The heterogeneous sedimentary structure controls the great spatial variability of the hydraulic

parameters in the aquifer (permeability, transmissivity, etc.) and the distribution of the most productive wells. All these variables in the aquifer can be used to explain its spatial and temporal response to hydraulic head changes and the subsequent vertical ground movements.

#### **3. Data and Methods**

A flowchart of the proposed method has been represented in Figure 2.

**Figure 2.** Flow chart of the proposed methodology.

#### *3.1. Data Employed and Their Origins*

The historical climate data for our case study were taken from Version 04 of the Spain02 dataset [44] (see Figure 3a,b). The Spain02 dataset includes daily temperature and precipitation estimates from observations (around 2500 quality-control stations) of the Spanish Meteorological Agency.

We also used data from individual global future climate projections (see Table 1) produced by regional climate models (RCM) nested within different global circulation models (GCM). This information was retrieved from the Coordinated Regional Downscaling Experiment (CORDEX) project [45]. The Spain02 project uses the same grids as the EURO-CORDEX project, which has a spatial resolution of approximately 12.5 km.

In addition, we also employed some hydrological information, including hydraulic head evolution in different piezometers obtained from the Spanish official network for monitoring the quantitative state of groundwater (see their location in Figure 1) and the historical recharge and pumping rates within the aquifer, which were taken from the information included in the Guadalquivir River Basin Plan (2015–2021).


**Table 1.** Regional climate models (RCMs) and global circulation models (GCMs) considered.

For those wells, we also used information about the land subsidence rates obtained by applying PInSAR techniques. The average vertical displacement based on the time series was estimated by using buffer areas with a radius of 1000 m for each piezometer and considering all PSI (persistent scatterer interferometry) data included in the area [15]. The percentage of fine-grained material (clay and silt) for each piezometer was obtained by exploiting geological data recorded in 38 boreholes drilled in the area. They were interpreted [15] to provide isolines of fine-grained sediment percentage, taking into account the clay and silt content in the first 50 m of the borehole, as the borehole depths were quite variable (from 50 to 122 m).

The results showed higher clay content in the piezometers located in the northern and southern extremes of the aquifer, as well as in its central part (see Table 2).

**Table 2.** Content of clay and silt in the considered piezometers.


The historical data and the future trends of population within the area were obtained from the Spanish Statistical Office <https://www.ine.es/>.

#### *3.2. Assessment of Subsidence from Satellite Data*

In the present study, spatial and temporal ground surface deformation assessment was conducted by exploiting 70 SAR images from three different satellites: ENVISAT (2003–2009, C band), Cosmo-SkyMed (2011–2014, X band), and Sentinel 1A (2015–2016, C band). The processing methodologies applied for each dataset were: (1) PSIG Cousins [46] for the ENVISAT and COSMO-Sky-Med datasets, and (2) direct integration approach [47] for the Sentinel 1A dataset. More detailed specifications can be found in Mateos et al. [15]. This combination allowed a thorough assessment of the ground deformation pattern in the aquifer and both the temporal and spatial dimensions of the subsidence.

PSInSAR measurements were obtained in the aquifer of the Vega de Granada (southern Spain), covering a large temporal span of 13 years (from 2003 to 2016). In this time, a severe drought affected the area during the ENVISAT period (2003–2006), and greater groundwater withdrawals took place.

PSInSAR data were correlated (temporally and spatially) with hydraulic head changes in the aquifer along the monitoring period, and with geological data (from boreholes) regarding the clay and silt content in the aquifer. Based on the borehole information, isolines of fine sediments percentage were obtained by Mateos et al. [15]. Clay and silt content is key information which can explain the spatial response of and aquifer system to hydraulic head changes and the subsequent vertical land movements.

#### *3.3. Definition of Local CC Scenarios*

A statistical method was used to define local future global change scenarios for the pilot, based on the historical information for the adopted reference period (1986–2015) and the available RCM simulations.

These scenarios were derived from RCM simulations available in the CORDEX project [45] for the most pessimistic emission scenario, RCP 8.5, and the temporal horizon 2016–2045. The future series were generated by applying the first- and second-moment correction technique under the bias-correction approach [24]. The bias-correction techniques applied a perturbation (transformation function) to the control series of the RCM simulations to obtain another series with statistics more similar to the historical series. The transformation function in the first- and second-moment correction technique is defined by focusing on the mean and standard deviation of the climate series. The same transformation function was applied to the future simulations of the RCM to obtain the climate change projections. An equifeasible ensemble of the individual climate change projections have been proposed in order to define more robust climate projections that are more representative than those based on a single model [32,48].

#### *3.4. Hydrological Impacts of Climate Change on Groundwater Levels*

A simple approach proposed by Scott [49] was applied to assess future CC impacts on global lumped drawdowns due to the future potential rainfall recharge and pumping. Following this approach, simple balance equations were applied in order to assess the global lumped change in hydraulic head (Δ*ht*) from aquifer storage (Δ*St*), which was calculated as follows.

$$
\Delta S\_t = R\_t - E\_t \tag{1}
$$

where *Rt* is the aquifer rainfall recharge and *Et* is the aquifer extraction, which can be obtained from Equation (2).

$$E\_t = Eag\_t - Enomag\_t \tag{2}$$

where *Eagt* represents the agricultural extractions and *Enonagt* the non-agricultural extractions. They were calculated using Equations (3) and (4), respectively.

$$\text{Eag}\_{t} = \text{Eag}\_{t-1} \left[ 1 + \frac{ET\_{t} - ET\_{t-1}}{ET\_{t-1}} \right] \tag{3}$$

$$Enonag\_t = Enonag\_{t-1} \left[ 1 + \frac{Pop\_t - Pop\_{t-1}}{Pop\_{t-1}} \right] \tag{4}$$

where *ET* is the evapotranspiration calculated using the Blaney–Criddle method (*ET* = *p*(0.46Tmean + 8)) on the basis of monthly temperature in ◦C (Tmean) and latitude-derived sunshine-hour fraction (*p*), and *Pop* is the population of the area.

Finally, the hydraulic head (Δ*ht*) was calculated using Equation (5).

$$
\Delta h\_t = \frac{\Delta S\_t}{A \times S\_y} \tag{5}
$$

where *A* is the aquifer area and *Sy* the specific yield.

The initial conditions used to simulate the recharge (R*t*−1) and pumping (agricultural or non-agricultural, *Eagt*−<sup>1</sup> and *Enonagt*−1) evolutions were taken from the information included in the Guadalquivir River Basin Plan (2015–2021).

The lumped approach proposed by Scott [49] was also applied to estimate the lumped hydraulic head drawdowns in the reference historical period (1986–2015). The method allowed estimation of the delta change (percentage increase) in the lumped aquifer drawdowns, taking into account the relative difference between the maximum lumped drawdowns in the historical and the future periods (2016–2045). These results were obtained under the assumption that a business-as-usual management scenario will be maintained in the future. The future potential hydraulic head in each piezometer was obtained by applying a delta change correction, using the lumped change to modify the historical evolution (for the reference period) of this variable within the piezometer.

#### *3.5. Propagation of Hydrological Impacts to Subsidence*

Simple linear regression models have been defined in order to approximate subsidence as a function of hydraulic head drawdowns in the selected head observation wells. We tested different transformations (Tr(X)) of the explanatory and target variables (logarithm, inverse, square, and square-root mathematical transformations) in order to identify the one that provided the best approximation to the empirical data for this problem (see Table 3).

The models used assume that there is a linear relation between both variables, the dependent variable and the explanatory variable and its transformations, which is reasonable if the deformation is elastic. An analysis of the linear correlation depending on the percentage of clay and silt content in the ground was also proposed in order to identify and discuss when linear regression might represent a better approach.

**Table 3.** Regression models and transformation of variables applied. The symbol \* represents the tested combinations of models and transformation of variables.


#### **4. Results and Discussion**

For the case study, the future equifeasible ensemble series obtained by applying the bias correction approach showed a mean global temperature increase of 7.72% and a mean global reduction of precipitation of 6.24%. Figure 3 shows the mean historical and future yearly series of precipitation and temperature.

The hydrological approach described in Section 3.4 was employed to propagate the generated future local climatic series to assess future hydrological impacts in terms of future recharge and pumping (Figure 4).

The mean reduction of the recharge and increase of pumping expected for the potential future horizon contemplated were 1.4 Hm<sup>3</sup> year−<sup>1</sup> and 1 Hm3 year<sup>−</sup>1, respectively. In the global budget, the impacts of CC on the recharge will have a higher influence on the future hydraulic head drawdowns.

Taking into account the projected changes in the future recharge and withdrawals within the aquifer, the global lumped hydraulic head drawdowns were obtained by applying the Scott [49] approach. The maximum lumped drawdown obtained in the future was 3.3% greater than the one obtained in the historical period. This relative change was employed to apply a delta change to correct the historical drawdowns in the selected head-observation wells (see Figure 1), obtaining for the future period the values presented in Figure 9.

In order to propagate the impacts of these hydraulic head drawdowns on the subsidence, different regression approaches were tested (see Figure 2). In each piezometer, the coefficient of determination of the calibrated models depended on the selected transformation (see results for the tested models in Figure 5). The two best combinations of model and transformation were S = a × P + b and <sup>S</sup><sup>2</sup> <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup><sup>2</sup> <sup>+</sup> b, where S represents the subsidence and P the hydraulic head (Models A and F in Figure 5). The mean R2 values of these models for the five considered piezometers were 0.60 and 0.73, respectively. These models were employed to predict the impacts of the potential future CC scenario on

subsidence. An example of the fit of the regression models for the target variable (subsidence) and the explanatory variable (hydraulic head) is included in Figure 6 for Piezometer 3 and the selected models.

**Figure 3.** Yearly mean historical and generated future series of precipitation (**a**) and mean temperature (**b**) for the periods 1986–2015 and 2016–2045. The historical data were obtained from the Spain02 dataset and the future series was generated by applying the proposed methodology.

**Figure 4.** Yearly mean historical and future series of recharge (**a**) and withdrawals (**b**) for the mean year in the periods 1986–2015 and 2016–2045. The historical data were obtained using the methodology proposed by Scott [49] and historical recharge and pumping rates from the Guadalquivir River Basin Plan (2015–2021). The future series were generated by applying the proposed methodology.

**Figure 5.** Coefficient of determination (R2) of the tested models for the five considered piezometers.

**Figure 6.** Relationship between the target variable and the explanatory variable for the two best linear regression models (S <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup> <sup>+</sup> b; S2 <sup>=</sup> <sup>a</sup> <sup>×</sup> P2 <sup>+</sup> b) for Piezometer 3.

The selected models were compared in terms of mean error and mean squared error too (see Figure 7). The two models showed a mean error of 0.0 mm; the model S = a × P + b showed a mean squared error of 20.63 mm2 and the model S2 <sup>=</sup> <sup>a</sup> <sup>×</sup> P2 <sup>+</sup> b showed a mean squared error of 23.44 mm2.

**Figure 7.** Mean values of the mean error and mean squared error for the five considered piezometers, obtained for the two best linear regression models (S <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup> <sup>+</sup> b; S2 <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup><sup>2</sup> <sup>+</sup> b).

The granulometry of subsurface sediments has a significant impact on the ground subsidence. We analyzed the influence of the percentage of clay and silt (in the surrounding area of the piezometer) on the historical subsidence rate and the linear behavior of the subsidence, assessed in terms of the R2 of the regression models.

Figure 8a shows the relationship between the percentage of fine-grained material and the historical subsidence rate, and Figure 8b the coefficient of determination obtained for the best approach for each piezometer vs. the percentage of clay and silt. In general, a higher percentage of fine-grained material was related to a lower subsidence rate, but the correlation of this relationship was poor (R<sup>2</sup> = 0.24). On the other hand, higher coefficients of determination were related with a more linear behavior, which was observed for the cases with lower percentage of clay and silt. Note that in general, the relationship between the percentage of fine-grained material and the coefficients of determination of the linear models had a good correlation (R2 higher than 0.8). In fact, PSInSAR results showed an inelastic deformation in the aquifer where a higher clay–silt content was identified [15]. Percentage of fine-grained material, thickness, and distribution of lenses significantly affected spatiotemporal subsidence patterns, which was in agreement with the results observed by other authors [10,50].

**Figure 8.** Relationship between the percentage of clay and silt and the historical subsidence rate (**a**) and the correlation coefficients obtained for the two best linear regression models (S = a × P + b; S2 <sup>=</sup> <sup>a</sup> <sup>×</sup> P2 <sup>+</sup> b) (**b**).

For Piezometer 6, which was the one with highest percentage of clay (70%), we obtained very low coefficient of determination (0.33 and 0.57 for the models S <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup> <sup>+</sup> b and S<sup>2</sup> <sup>=</sup> <sup>a</sup> <sup>×</sup> P2 <sup>+</sup> b, respectively) and the linear regression approach was not appropriate. For this reason, we have not included the results for the assessment of the future subsidence for this observation well.

The propagation of the impacts of the potential future CC scenario on subsidence (Figure 9) showed important increases of the maximum subsidence (55.3% and 52.7% for the models S = a × P + b and S<sup>2</sup> <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup><sup>2</sup> <sup>+</sup> b, respectively) with respect to the historical maximum observed values. The highest increase of the maximum subsidence (68.3% and 65.7% for the models S <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup> <sup>+</sup> b and S2 <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup><sup>2</sup> <sup>+</sup> b, respectively) with respect to the historical maximum observed values occurred in Piezometer 5, which was located in the western area where the percentage of clay was 40%.

In terms of mean subsidence (see Figure 10) the mean increase of subsidence was 4.1 mm and 4.5 mm for the models S <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup> <sup>+</sup> b and S<sup>2</sup> <sup>=</sup> <sup>a</sup> <sup>×</sup> P2 <sup>+</sup> b, respectively. The piezometer with the highest increase of mean subsidence was Piezometer 3, with 5.4 mm and 5.9 mm for the models S = a × P + b and S<sup>2</sup> <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup><sup>2</sup> <sup>+</sup> b, respectively. This piezometer showed the highest variability of subsidence in the historical and future periods, and had the lowest content of clay and silt.

**Figure 9.** Historical and future subsidence and hydraulic head for the two best linear regression models (S <sup>=</sup> <sup>a</sup> <sup>×</sup> <sup>P</sup> <sup>+</sup> b; S<sup>2</sup> <sup>=</sup> <sup>a</sup> <sup>×</sup> P2 <sup>+</sup> b).

**Figure 10.** Historical and future subsidence for the two best linear regression models (S = a × P + b; S2 <sup>=</sup> <sup>a</sup> <sup>×</sup> P2 <sup>+</sup> b).

#### *Hypotheses Assumed and Limitations*

The assessment performed assumed some hypotheses or made simplifications as follows.

Climate change scenarios


Propagation of climate change in terms of hydrological impacts


Propagation of hydrological impacts to subsidence


#### **5. Conclusions**

A method to assess impact of potential future CC scenarios on land subsidence related to groundwater level depletion in detrital aquifers was described and applied in a case study. It does not require a previous distributed groundwater flow model of the aquifer, which is an important advantage for its applicability in cases with limited information. Taking into account the uncertainties around future potential CC scenarios, it could provide a useful first assessment of its impacts on subsidence. It allows analyses of wide areas to be performed in order to identify potential hot spots that require more exhaustive analysis. It will help to assess sustainable adaptation strategies in identified vulnerable areas, taking into account subsidence issues during drought-critical periods. The methodology was applied in the Vega de Granada aquifer (Granada, SE Spain). Good correlation between groundwater level depletion and the subsidence was obtained in the wells where the percentage of clay was below 50%. The analysis of results showed that, assuming a business-as-usual management scenario, the impacts of CC on subsidence would be very significant for the case study. The mean increase of the maximum subsidence rates in the considered wells for the future horizon (2016–2045) and the RCP scenario 8.5 was 54%. In order to avoid undesirable consequences/risks as observed in other regions worldwide (Jakarta, Ho Chi Min City, and Bangkok, among others) where land subsidence is causing severe impacts—permanent inundation of land, aggravated flooding, changes in topographic gradients, rupture of the land surface, structural damage to buildings and infrastructures, and reduced capacity of aquifers to store water—some adaptation strategies should be applied to control and minimize the land subsidence caused by groundwater withdrawals.

**Author Contributions:** D.P.-V. and R.M.M. conceived and designed the research; P.E. analyzed the subsidence data; A.-J.C.-L. analyzed the data and conducted the experiments. All authors contributed to writing the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially supported by the research projects SIGLO-AN (RTI2018-101397-B-I00) from the Spanish Ministry of Science, Innovation and Universities (Programa Estatal de I+D+I orientada a los Retos de la Sociedad) and the GeoE.171.008-TACTIC project from GeoERA organization funded by European Union's Horizon 2020 research and innovation program.

**Acknowledgments:** We would like to thank the Spain02 and CORDEX projects, the ENVISAT, Cosmo-SkyMed and Sentinel satellites, the Spanish Statistical Office, the Spanish official network for monitoring the quantitative state of groundwater and the Guadalquivir River Basin authority for the data provided for this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Approaches for Estimating Daily River Leakage from Arid Ephemeral Streams**

#### **Leilei Min 1,2,**†**, Peter Yu. Vasilevskiy 3,**†**, Ping Wang 2,4,\*, Sergey P. Pozdniakov <sup>3</sup> and Jingjie Yu 2,4**


Received: 27 December 2019; Accepted: 9 February 2020; Published: 12 February 2020

**Abstract:** Despite the significance of river leakage to riparian ecosystems in arid/semi-arid regions, a true understanding and the accurate quantification of the leakage processes of ephemeral rivers in these regions remain elusive. In this study, the patterns of river infiltration and the associated controlling factors in an approximately 150-km section of the Donghe River (lower Heihe River, China) were revealed using a combination of field investigations and modelling techniques. The results showed that from 21 April 2010 to 7 September 2012, river water leakage accounted for 33% of the total river runoff in the simulated segments. A sensitivity analysis showed that the simulated infiltration rates were most sensitive to the aquifer hydraulic conductivity and the maximum evapotranspiration (ET) rate. However, the river leakage rate, i.e., the ratio of the leakage volume to the total runoff volume, of a single runoff event relies heavily on the total runoff volume and river flow rate. In addition to the hydraulic parameters of riverbeds, the characteristics of ET parameters are equally important for quantifying the flux exchange between arid ephemeral streams and underlying aquifers. Coupled surface/groundwater models, which aim to estimate river leakage, should consider riparian zones because these areas play a dominant role in the formation of leakage from the river for recharging via ET. The results of this paper can be used as a reference for water resource planning and management in regulated river basins to help maintain riparian ecosystems in arid regions.

**Keywords:** river-aquifer interaction; numerical simulation; sensitivity analysis; MODFLOW; Heihe River

#### **1. Introduction**

Surface water and groundwater are important components of the terrestrial water cycle, and their interaction forms the surface morphology, controls the material and energy fluxes in the subsurface zone, and affects the riparian ecosystem [1,2]. However, as indicated by Sophocleous [3], the interactions between surface water and groundwater are complex, and obtaining a deep understanding of these interactions in relation to the climate, landform features, geology, and biotic factors remains a great challenge. Clearly, the exchange between surface water and groundwater is likely to become even more challenging due to the impacts of human activities and climate change [4,5], which caused the disappearance of approximately 90,000 km<sup>2</sup> of permanent surface water between 1984 and 2015 [6].

Given that a streambed acts as the physical interface between the surface and subsurface of a stream [7,8], the hydraulic properties of streambeds mainly control the interactions between the stream and the underlying aquifer [4,9–11]. However, as confirmed by numerous field investigations (e.g., [12–14]), the hydraulic properties of streambeds usually exhibit large spatial and temporal variations, which are mainly caused by continuous changes in the streambed properties (e.g., the topography or hydraulic conductivity) during erosion and sedimentation processes [7]. Significant changes in the hydraulic properties of streambeds may even occur during short flooding events [15]. Additionally, the thermal dynamics of streambeds induced by diurnal and seasonal fluctuations in the stream water temperature also greatly influence the hydraulic properties of the streambeds [16–18]. For intermittent rivers, which constitute more than 30% of the total length and discharge of the global river network [2,19], the hydraulic properties of streambeds are even more variable due to constant alternation between dry and wet conditions [1,16].

The states of the connections between streams and underlying aquifers exert another important influence on the flux exchange between surface water and groundwater [4]. For losing streams, when stream-aquifer systems are transformed from connected to disconnected systems, the lateral flow induced by capillarity or heterogeneity plays a vital role in the stream water and groundwater interaction [20–22]. Brunner et al. [23] provided a theoretical criterion for justifying the connection/disconnection states between a stream and the underlying aquifer and suggested that the disconnection problem could be solved via a fully coupled, variably saturated flow model. Therefore, in addition to the broad range of field methods (e.g., [24–27]) and associated analytical solutions e.g., [17,28,29], numerical simulations have been widely applied to investigate stream–aquifer interactions at different scales because these simulations can analyze the influences of transient flows and streambed heterogeneity on surface–groundwater exchanges [4,30].

Recently, interest in the interactions between intermittent streams and groundwater in arid and semi-arid regions has been continuously growing due to the unique role of these interactions in shaping fragile riparian ecosystems (e.g., [27,31–35]). It is clear that stream water leakage is the dominant recharge mechanism in such regions; however, the infiltration processes during various stream discharge patterns and the factors that control the stream–aquifer interactions in typical losing connected/disconnected river systems remain unclear. In this context, it is critical to quantify the potential impacts of the variations in the stream width and leakage coefficient [36,37], which vary greatly for intermittent rivers in arid regions [1], on stream–aquifer interactions.

The lower Heihe River Basin represents a typical extremely arid region in north-western China, where the mean annual precipitation is less than 50 mm, while the mean annual evaporation can reach 1500 mm [38–40]. The lower Heihe River is characterized by intermittent streams, and the streambeds usually remain dry from April to June [41]. The hydraulic property dynamics of these streambeds were investigated in detail in recent studies (e.g., [13,16,41]), and the monthly river leakage was approximately estimated using the River (RIV) package of MODFLOW-2005 [42,43]. While this estimation was based on regional groundwater modelling, it lacked a detailed analysis of the influence of the stream/streambed dynamics on the water exchange between the stream and aquifer. Therefore, to fill this gap, the objectives of this study were to (1) quantify the daily river leakage rates by numerically simulating the flux exchange between the rivers and the underlying aquifers and (2) identify the predominant factors that control the river–aquifer interactions in intermittent dryland rivers using parameter sensitivity analysis.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area is located in the lower reaches of the second largest inland river in north-western China (Figure 1) and is characterized by a hyper-arid climate with an annual precipitation of only 35 mm and an annual potential evaporation of approximately 1500 mm [38,40]. Over the period of 1961–2015, the mean annual air temperature was +9.09 ◦C, with a minimum monthly mean air

temperature of −11.23 ◦C in January and a maximum monthly mean air temperature of +27.05 ◦C in July [16]. The topography of this area gradually declines from the southwest to the northeast, with an average slope of 1–3‰, and the elevation is between 1127 and 820 m [44]. The dominant landscape is the Gobi Desert, which is composed of wind-eroded hilly land, desert, and alkaline soils [45].

**Figure 1.** Study area and simulation domain. LHS represents the Langxinshan hydrometric station.

The study area is located in a regional tectonic basin where the bedrock is composed of Sinian (Z) and Late Jurassic (J3) formations. From the southwest to the northeast of the study area, the regional Quaternary (Q) aquifer system varies from a zone that consists of unconfined gravel and pebbles to a multi-layered zone that consists of sand and silt with a depth of several hundred metres [46,47]. In general, the phreatic aquifer is recharged by river water, and the groundwater flows from the southwest towards the terminal lakes and is then discharged via evaporation [39,48].

The lower Heihe River, which is divided into two losing streams at the Langxinshan hydrometric station (LHS), i.e., the Donghe and Xihe rivers, flows through the Gobi Desert before entering terminal lakes (the East and West Juyan lakes) (Figure 1). At the LHS, the surface flow to the Donghe and Xihe rivers is regulated by a system of sluices. These river systems are the primary sources of shallow groundwater recharging via riverbed infiltration [39,49] due to the relatively high vertical hydraulic conductivity [13]. The limited vegetation in the region is distributed along the rivers and relies on surface water and shallow groundwater in the riparian zone for sustenance [50–52]. More detailed descriptions of the study area are presented in Wang, Yu, Zhang, and Liu [48]; Wang, Yu, Pozdniakov, Grinevsky, and Liu [39]; and Yao, Zheng, Liu, Cao, Xiao, Li, and Li [42].

#### *2.2. Simulation of River Water and Groundwater Interactions*

#### 2.2.1. Simulation Domain and Boundary Conditions

As noted by Wang, Pozdniakov and Vasilevskiy [16], approximately 71% of the total runoff was allocated to the Donghe River at the LHS from 1988 to 2015 (Figure 1). Additionally, the characteristics of the streambed sediment formation in the Donghe River are typical of the study area [13]. Therefore, the Donghe River was selected to analyze the water exchange processes between the river and the underlying aquifer. The simulation domain was determined according to the surface water and groundwater interaction zone during river flow events. Based on previous studies (e.g., [39,42,48,53]), the eastern and western boundaries of the model were determined by the regional flow direction, which was generally parallel to the river channel and approximately 10 km from the Donghe River. Thus, the eastern and western boundaries were generalized as no-flow boundaries. The southern and northern boundaries were at the LHS and Angcizha, respectively. The simulation domain is presented in Figure 1, and the total area of the simulation region was 2306.25 km2.

To focus on the river–aquifer interactions, the upper unconfined aquifer was simulated as spatially heterogeneous single-layer aquifer system. The top and bottom elevations of the aquifer are shown in Figure 2b,c. Based on geochemical data [48,54] and the regional modelling of the groundwater flow system [42], the southern and northern boundaries of the model were designed as general head boundaries (GHBs) [55] in accordance with lateral groundwater flow from adjacent groundwater basins into the study area. The flux at the top boundary was determined according to the meteorological conditions, surface water and irrigation water infiltration, and evapotranspiration (ET) processes. The bottom of the aquifer was considered impermeable; therefore, the bottom boundary was designed as a no-flow boundary.

**Figure 2.** Model setup: (**a**) parameter division; (**b**) land surface elevation; (**c**) bottom elevation of the aquifer; (**d**) initial water levels; and (**e**) river segment division and monitoring sites.

2.2.2. Field Observations and Model Parameterization

The daily groundwater levels were monitored at five observation wells (G1, G2, G3, G4, and G5) located along the riverbank (Figure 2e). The daily river water levels were also observed at wells R1, R2, R3, and R4, which were installed along the Donghe River (Figure 2e). The water levels were recorded by Schlumberger Mini-Diver pressure transducers (Eijkelcamp, EM Giesbeek, The Netherlands). The water level measurements compensated for barometric changes, which were measured by Schlumberger Baro-Divers. The uncertainty of the water level measurements was ±5 mm. In addition, the daily streamflow data at the LHS and daily E-601 pan evaporation rates at the Ejina meteorological station (Figure 1) were available for the period of 2010–2012.

Based on previous hydrogeological investigations [46] and numerical modelling [43], four aquifer zones with different hydraulic parameters were defined (Figure 2a, Table 1). The initial values of hydraulic conductivity varied from 6 to 21 m/day, and the initial specific yield was set to 0.15.

**Table 1.** Initial and calibrated values of the hydraulic conductivity (*K*) and specific yield (*Sy*) of the aquifer.


Based on the riverbed hydraulic conductivities measured in previous studies [13,16], we divided the Donghe River into eight segments, as shown in Figure 2. The parameters of each river segment, including the river width, thickness, and hydraulic conductivity, are listed in Table 2.


**Table 2.** River width (*Lr*), riverbed *k*0/*m*0, and riverbed hydraulic conductance (*C*) of the river segments.

#### 2.2.3. Numerical Simulations

A three-dimensional finite difference-based groundwater flow model from MODFLOW-2005 [55] was used in the pre- and post-processing modelling environment of Processing MODFLOW [56] to simulate the saturated subsurface flow and surface water exchanges via the river–aquifer interface. In the present study, surface water leakage from the intermittent river was considered the main source of recharge for the aquifer [48]. The riverbed conductance (*C*, L2T−1) is a key parameter that controls the interaction between the surface water and groundwater and is represented in MODFLOW [55,57] by the following equation:

$$\mathcal{C} = D\_{\bar{r}\bar{w}} L\_{\bar{r}} k\_0 / m\_0 \tag{1}$$

where *Driv* is the length of the river reach within the grid cell (L), *Lr* is the river width (L), *m*<sup>0</sup> is the thickness of the riverbed (L), and *k*<sup>0</sup> is the hydraulic conductivity of the riverbed (LT<sup>−</sup>1).

The Streamflow Routing (STR) package [58] was selected to simulate the interactions between the river water and groundwater. This package was able to simulate the major features of the surface water in this study, i.e., changes in the flow along the river due to interactions with groundwater. In the STR package, the flow in a stream is instantaneously routed downstream. The streamflow routing is designed through a network of streams and always flows in the same direction along the streams. The stream stages (*Hr*, L) of a rectangular stream channel are calculated using Manning's equation as follows [56]:

$$H\_r = \left(\frac{\mathcal{Q} \cdot n}{L\_r \cdot S\_{riv}^{\frac{1}{2}}}\right)^{\frac{1}{2}}\tag{2}$$

where *Q* is the calculated river discharge (L3/T), *Sriv* is the slope of the river channel (L/L), and *n* is Manning's roughness coefficient (-).

The amount of recharge from precipitation was insignificant, and the total precipitation over the period of 2010–2012 was only 91 mm, based on observations at the local metrological station. For arid regions in north-western China, the direct recharging of rain-fed groundwater was estimated by the chloride mass balance method to be 1.5% of the mean annual precipitation [59,60]. In the study area, precipitation infiltration is most likely negligible, as indicated by the fact that single rainfall events of more than 5 mm were extremely rare between 2010 and 2012 and occurred predominantly during summer and autumn, when the potential evapotranspiration (PET) was extremely high (Figure 3). The PET was obtained from the water surface evaporation data observed by an E-601 evaporator at the local meteorological station (Figure 1) during the non-freezing period (April to October) and the calculated PET from Du, Yu, Wang, and Zhang [38] during the freezing period (November to the following March).

**Figure 3.** Daily potential evapotranspiration (PET, mm) and precipitation (P, mm) at the Ejina meteorological station.

Groundwater ET is a predominant discharge term in the water budget in the study area [39]. The depth of the water table in this area is mostly between 2 and 4 m [45]; therefore, the relationship between ET and water table depth can be simply assumed to be linear [61] with an extinction depth of 5 m [62]. For this reason, the Evapotranspiration (EVT) package [55], which assumes a linear relationship between ET and water table depth, was used to simulate the processes of ET. As shown in Figure 1, excluding the dominant landscape of the Gobi Desert, there is a large area of riparian oasis within the model domain. To simulate plant transpiration in the riparian oasis and soil evaporation in the Gobi Desert, we assigned different maximum ET rates (*ETmax*) for these two landscape types. The monthly *ETmax* values for the riparian oasis and Gobi Desert were independently calibrated to match the linear relationship between ET and water table depth within the interval of 2–4 m, while the

monthly variations followed the monthly measured and calculated PET rates at the Ejina meteorological station (Figure 3).

Groundwater extraction by pumping wells for agricultural irrigation occurs primarily in the southern and north-eastern areas of the study region (Figure 1). Based on field investigations, the farmland area irrigated by a single irrigation well was approximately 1.08 <sup>×</sup> 105 m2, and the daily water pumping rate of a farmland irrigation well was approximately 300 m3/day during the irrigation period from May to August. The fraction of water returned to the aquifer as a result of agricultural irrigation was estimated to be approximately 0.1 [43]. Therefore, the water withdrawn from the aquifer for irrigation was 625 m3/day within a grid cell of 2.50 <sup>×</sup> 10<sup>5</sup> m2. The Well package implemented in MODFLOW-2005 [55] was used to simulate this groundwater withdrawal from the aquifer. The total number of cells containing irrigation wells within the simulation domain was 182.

Previous studies indicated that vertical unsaturated flow is insignificant in the Gobi Desert based on the water exchange between rivers and aquifers [63,64]. In addition, the dominant vegetation in the riparian area comprises groundwater-dependent species (e.g., *Populus* and *Tamarix*), which mostly use groundwater for transpiration [50]. Therefore, unsaturated flow in the study area can be neglected and was not addressed in the simulations.

The numerical model consisted of 225 rows and 41 columns, with a total simulation domain of 112.5 × 20.5 km. The simulation period was from 21 April 2010 to 7 September 2012, and the temporal resolution was 1 day. We used the measured daily river flow at the LHS from 1 January to 21 April 2010 and ran the model to obtain the initial groundwater levels. The simulated groundwater levels on 21 April 2010 were set as the initial groundwater levels (Figure 2).

#### 2.2.4. Model Calibration

We used the combination of the measured river water (R1–R4) and groundwater (G1–G5) levels to calibrate the hydraulic conductivity, specific yield of the aquifer, riverbed hydraulic conductance, and GHB hydraulic conductance. The root mean squared error (*RMSE*, m) and correlation coefficient (*R*) were used to evaluate the goodness of fit of the observed and calculated levels [65]. The calibrated parameters of the aquifer and streambed are listed in Tables 1 and 2, respectively. Notably, the river in this study is intermittent with a wide riverbed, and the river width is highly dependent on the river stage. To account for the influence of the river width on the river–aquifer exchange, we used the temporal variations in the riverbed hydraulic conductance (see Table 2) associated with each river flow event in the simulations. In addition, the GHB hydraulic conductance was also calibrated, with values of 2000 m2/day for the southern boundary and 2200 m2/day for the northern boundary.

As shown in Figure 4, the simulated and observed water levels at the nine monitoring wells are generally consistent. The simulated water levels reflect the variations in the water levels at the observation points. The differences between the calculated and observed water levels are less than 0.5 m. The *RMSE* varies from 0.04 to 0.31, depending on the observation well. The correlation coefficients between the simulated and observed water levels at most monitoring wells are relatively high (0.77–0.95), with the exception of that at well G2, where no correlation is observed. This finding can be explained by the relatively large distance of well G2 from the Donghe River (Figure 2). Thus, the groundwater level at well G2 is less affected by the fluctuations in the river water level compared to the levels at the other monitoring wells, which are located closer to the river. The slight seasonal variations in the groundwater level at well G2 are more likely affected by the ET process. Subsequent research should aim to address simulating ET processes using a nonlinear relationship between *ETmax* and the groundwater level depth (e.g., [61,62]) to enhance the conformity of the observed and calculated levels in areas far from rivers.

**Figure 4.** Simulated versus observed water levels from 21 April 2010 to 7 September 2012.

#### **3. Results**

#### *3.1. Groundwater Budget*

During the period from 21 April 2010 to 7 September 2012, the total groundwater recharge within the simulated region was 4.26 <sup>×</sup> 108 m3. The river leakage through the streambed to the aquifer was 3.59 <sup>×</sup> 108 m3, which accounted for approximately 84% of the total groundwater recharge. The other 16% of the groundwater recharge was due to lateral flow through the southern and northern model boundaries (6.77 <sup>×</sup> <sup>10</sup><sup>7</sup> m3), which can be evidenced by hydro-geochemical analyses [48,54] and regional groundwater flow simulations [42].

The total groundwater discharge within the simulated region was 2.8 <sup>×</sup> 108 m3, and the groundwater ET was 1.64 <sup>×</sup> 108 m3, which accounted for 59% of the total discharge. Groundwater discharge to the river (5.53 <sup>×</sup> 107 m3), groundwater exploitation (4.19 <sup>×</sup> 107 m3), and water outflow through the southern and northern model boundaries (1.85 <sup>×</sup> 10<sup>7</sup> m3) accounted for approximately 20%, 15%, and 6% of the total discharge, respectively.

The groundwater budget analysis indicated that river water leakage was the main recharge source, and that groundwater ET was the predominant discharge type (Table 3). As noted by Wang et al. [66] and Wang et al. [67], these two processes determine the changes in water storage and control the spatial and temporal dynamics of the studied groundwater system. The difference between the total recharge and discharge (1.46 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup>3) can be explained by the increased groundwater storage in the model domain from 21 April 2010 to 7 September 2012. These results were consistent with previous studies (e.g., [39,45]), which demonstrated that, especially in the riparian zone, the groundwater level increased as a result of environmental flow controls aimed at delivering a set amount of surface water to the study area after 2000.


**Table 3.** Groundwater budget of the model domain from 21 April 2010 to 7 September 2012.

As shown in Figure 1, oasis area accounts for approximately 20% of the total simulation area, and the other 80% represents the Gobi Desert, where groundwater ET processes are negligible [63]. The growing season lasts from May to September and accounts for approximately 150 days per year. As a result, the average ET rate in the oasis area during the growing season is approximately 0.79 mm/day, which is close to *Tamarix*'s ET rate of 0.63–0.73 mm/day, as estimated using water table fluctuation methods during the growing seasons of 2010–2012 in this area [68].

#### *3.2. Daily River Leakage*

Daily river leakage through the streambed was highly dependent on river inflow, and river leakage generally followed river runoff processes (Figure 5). As shown in Figure 5, after high-flow events, e.g., from 18 August 2011 to 3 November 2011, the river leakage through the streambed was negative, which indicated that groundwater was discharged to the river. This finding can be explained by the fact that the riverbank stored surface water during the high stage, and when the river flow decreased rapidly, the groundwater stored in the riverbank was released into the river channel.

**Figure 5.** Observed daily runoff of the Donghe River at the LHS and simulated daily river leakage from 21 April 2010 to 7 September 2012.

As the total runoff of the Donghe River at the LHS was 9.19 <sup>×</sup> 108 m3, the river leakage rate, i.e., the percentage of river leakage divided by the total inflow, in the river reaches from the LHS to the

Angcizha was approximately 33% from 21 April 2010 to 7 September 2012. However, the average leakage rates during the eight individual flow events varied significantly from 24% to 99.8% (Table 4).


**Table 4.** Leakage rates of the river flow periods during the simulation period; correlation coefficients indicate relationships between the stream runoff and river leakage.

For events with short-term flow durations (e.g., from 21 April to 9 May 2010 and 13 to 17 July 2011), the river leakage rates were relatively high, with values of 42–99.8%. In contrast, for events with long-term flow durations (e.g., from 18 August to 3 November 2011 and 4 December 2011 to 24 May 2012), the river leakage rates were relatively low, with values of 24–38%. This difference is mainly associated with the decreasing hydraulic gradient between the surface water and groundwater during events with long-term flow durations.

Table 4 shows the correlation coefficients between the stream leakage and river runoff for eight individual flow events. These correlation coefficients are relatively high (from 0.52 to nearly 1, except for the flow event from 6 January to 13 May 2011 with a value of 0.25). In addition, there is a common tendency that shorter flow events during summer periods exhibit higher correlations. For example, short flow events, e.g., 13–17 July 2011, 21 April–9 May 2010, and 12–26 July 2010, have correlation coefficients of 0.79–0.99, while relatively long flow events, e.g., 6 January–13 May 2011, 3 August–7 September 2012, and 4 December 2011–24 May 2012, have correlation coefficients of 0.25–0.61. Therefore, it can be stated that if water needs to be transported from Langxinshan to terminal lakes with minimum leakage, the water should be transported during winter periods. However, to increase the groundwater level in the riparian zone and support riparian ecosystems, it is better to transport water predominantly during spring and summer periods.

From the perspective of ground and surface water interactions, the interaction regime is constantly connected, which means that even after continuous periods without runoff (for example, from April to June), the groundwater level remains above the bottom of the riverbed sediments. Thus, water from the stream percolates via saturated sediments, which contributes to relatively high leakage rates.

#### *3.3. Parameter Sensitivity Analysis*

To investigate the influence of the model parameters on the simulation results, an uncertainty analysis was conducted. The present study revealed the model sensitivity by analyzing the effects of the ET rate, aquifer parameters, and streambed parameters on the cumulative river leakage.

Sensitivity represents the effect of variations in one parameter on the calculation results and is generally represented by the following equation [69]:

$$\beta\_k = \frac{\Delta Q}{\Delta \alpha\_k} \approx \frac{Q(\alpha\_k + \Delta \alpha\_k) - Q(\alpha\_k)}{\Delta \alpha\_k} \tag{3}$$

where β*<sup>k</sup>* represents the sensitivity of the model variable (leakage (*Q*) in the present study) to the parameter, α*<sup>k</sup>* represents the actual parameter value, Δα*<sup>k</sup>* represents the change in the parameter value, *Q*(α*<sup>k</sup>* + Δα*k*) represents the leakage simulated by the model after the parameter variation and *Q*(α*k*) represents the leakage simulated by the model before the parameter variation.

A local sensitivity analysis method was used to analyze the sensitivity of the model parameters to river leakage. Specifically, we changed the hydraulic conductivity (*K*) and specific yield (*Sy*) of the aquifer, the maximum ET rate (*ETmax*), and the riverbed conductance (*C*) by −20% to 20%. The parameter range selected for the sensitivity analysis was based on in situ experiments of streambed hydraulic conductivity [13,16] and the observed potential evaporation in the study area (Figure 3). Only one parameter was changed each time to determine the leakage variation (Figure 6). The results of the sensitivity analysis showed that the hydraulic conductivity of the aquifer had the greatest effect on leakage. When the hydraulic conductivity changed by 20%, the leakage changed by approximately 11%. Thus, the sensitivity of leakage to the hydraulic conductivity of the aquifer was 0.57. The calculations indicated that the sensitivities of leakage to the maximum ET rate, specific yield, and riverbed conductance were 0.28, 0.08, and 0.02, respectively. Therefore, leakage was most sensitive to the hydraulic conductivity of the aquifer, followed by the maximum ET rate, specific yield of the aquifer, and riverbed conductance.

**Figure 6.** Sensitivity of river leakage to the model parameters.

The significant sensitivity of the river leakage to the maximum ET rate can be explained by the increase in the hydraulic gradient between the surface and groundwater when *ETmax* rises, which causes a decline in the groundwater level. Figure 7 demonstrates that if we enhance *ETmax*, the hydraulic gradient rises, causing the leakage rate from the river to increase. Thus, according to the simulation results, riverbank recharge processes are highly dependent on riparian ET processes.

**Figure 7.** Groundwater levels under different values of the maximum evapotranspiration rate (cross-section from west to east in the middle of the simulated area).

#### **4. Discussion**

To explain the high sensitivity of river losses to groundwater ET in the riparian area, we consider a supporting analytical model for the formation of groundwater flow from the river towards the ET zone in the riparian area adjacent to the river channel. For this purpose, let us consider the one-dimensional groundwater flow formed by losses from the river and discharged by the ET of groundwater in the adjacent riparian zone. The details of the formulation of the corresponding groundwater flow boundary problem are given in Appendix A. As shown in Appendix A, water losses from the river during stationary groundwater flow periods are controlled by two hydraulic parameters expressed in length units. The first parameter, which is the streambed hydraulic resistance length, Δ*L*, characterizes

the additional hydraulic resistance due to the bottom sediments of the river channel [36]. The length Δ*L* is expressed as follows:

$$
\Delta L = b\_r^{-1} \coth(0.5 L\_r b\_r);
\
b\_r = \sqrt{\frac{k\_0}{m\_0 T}}\tag{4}
$$

The second parameter is the characteristic width, *LET*, of the groundwater ET zone adjacent to the river channel. This length depends on the parameters characterizing the decrease in ET with the depth of the groundwater table and is expressed as follows:

$$L\_{ET} = \sqrt{\frac{d\_{\text{max}}T}{ET\_{\text{max}}}};\ d\_{\text{max}} = Z\_{\text{surf}} - Z\_{\text{crit}} \tag{5}$$

The expression for the specific flow rate of losses from the river, *q*, to the riparian area is as follows:

$$q = 2T \frac{H\_r - Z\_{\rm crit}}{L\_{\rm ET} + \Delta L} = 2T \frac{H\_r - Z\_{\rm crit}}{L\_{\rm ET}(1 + a\_r)}; \ a\_r = \sqrt{\frac{ET\_{\rm max} m\_0}{d\_{\rm max} k\_0}} \text{coth}(0.5 L\_{\rm r} b\_{\rm r}) \tag{6}$$

The factor of 2 in Equation (6) assumes that there will be symmetrical losses to both riverbanks with riparian zones. The notations for Equations (4)–(6) are listed in the Appendix A.

The equations shown above help to clarify the fact that the sensitivity of leakage to riverbed conductance is only 0.02, while that to the hydraulic conductivity of the aquifer is 0.57, which can be explained by the value of the streambed hydraulic resistance length (Δ*L*), which is between 60 and 80 m, depending on the number of river segments (Table 2). The parameter Δ*L* characterizes the hydraulic resistance of the bottom sediments (Appendix A). In the study area, the hydraulic conductivity of the riverbed sediments is relatively high, i.e., 1–40 m/day [41], which leads to relatively small bottom sediment hydraulic resistance values.

In addition, an analytical criterion, α*<sup>r</sup>* (Appendix A), can be calculated using typical values of *ETmax* during the vegetation period as a function of *k*0, *m*0, and *dmax* (extinction depth). The calculated value of α*<sup>r</sup>* is 0.01, which indicates almost no sensitivity of leakage to the streambed hydraulic conductance. However, the sensitivities of leakage to the hydraulic conductivity of the aquifer and *ETmax* are comparable. Thus, under the present conditions, river leakage is more sensitive to the hydraulic parameters related to the aquifer and ET rather than to the riverbed, and thus, it is more important to study aquifer parameters and ET parameters to accurately estimate river leakage from such intermittent streams under arid conditions. The present conclusion is valid for cases with high maximum ET rate values and relatively low Δ*L* values.

Equally important is that the EVT package [55] selected in this study to simulate ET processes is oversimplified and is highly based on previous empirical analysis of the dependence of ET on the water table depth between 2 and 4 m [45,68]. This approach is probably acceptable for analyses of surface–groundwater exchanges at regional scales. However, to address riparian ET processes, nonlinear ET models that account for plant types with different rooting depths (e.g., [70]), or even physical-based models with dynamic root optimization schemes (e.g., [71]), are required.

#### **5. Conclusions**

The present study conducted coupled simulations of surface water and groundwater exchanges in an arid ephemeral stream–aquifer system using the MODFLOW-2005 code with the STR package. The results showed that from 21 April 2010 to 7 September 2012, river water leakage accounted for 33% of the total river runoff for the simulated segments.

A sensitivity analysis demonstrated that the most important parameters of the studied system that influence river leakage are the hydraulic conductance of the aquifer and the maximum ET rate. Almost no sensitivity was obtained for the riverbed hydraulic conductance, which was explained by the relatively high hydraulic conductivity of the riverbed sediments. Thus, instead of studying

the hydraulic parameters of riverbeds, further research should focus on studying the ET parameters and selecting an appropriate ET model that reflects the eco-physiology of riparian ecosystems [70]. The present conclusion is valid only for cases with relatively high streambed hydraulic conductivities (compared to the hydraulic conductivity of the aquifer) in arid regions, as demonstrated in the studied case. Coupled surface/groundwater models, which are used to estimate river leakage, should consider riparian zones because they play a dominant role in the formation of leakage from rivers for recharging via ET.

As the model synchronously simulated the daily variation in the river water and groundwater levels that affected leakage, the simulation results are more reliable than those of previous models that used only groundwater level data collected over long periods for verification when simulating leakage (temporal resolution greater than 10 days). To the best of our knowledge, the present study is the first to simulate and analyze the daily river leakage process under the conditions of ecological water transport in the downstream region of the Heihe River. The study results can provide scientific evidence for further ecological water transport research.

**Author Contributions:** Conceptualization, J.Y., P.W., and S.P.P.; methodology, L.M., P.Y.V., and S.P.P.; validation, P.W.; data curation, L.M. and P.W.; writing—original draft preparation, L.M., and P.W.; writing—review and editing, P.Y.V., P.W., and S.P.P.; supervision, J.Y.; and funding acquisition, P.W., J.Y., P.Y.V., and S.P.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by grants from the National Natural Science Foundation of China (Nos. 41671023 and 41877165) and the NSFC-RFBR Programme 2018–2019 (Nos. 41811530084 and 18-55-53025 ГΦEH\_a), and the second author was supported by the Russian Foundation for Basic Research (No. 19-35-90014 Acпиpaнты).

**Acknowledgments:** We would like to thank Grinevsky Sergey of the Department of Hydrogeology, Lomonosov Moscow State University, for his helpful advice and Juntao Zhu, Yongliang Xu, Zhiyong Wang, and Dandan Wang for their participation in the fieldwork. Ping Wang and Sergey P. Pozdniakov are grateful for support from the Special Exchange Program of the Chinese Academy of Sciences 2019–2020. The authors also gratefully acknowledge the anonymous reviewers for their valuable comments and suggestions that have led to substantial improvements over an earlier version of the manuscript.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**

Groundwater flow from a river to a riparian evapotranspiration (ET) zone.

In an arid region with a groundwater ET zone located along a stream in summer (Figure A1), river water infiltrates the aquifer. This water is discharged within the riparian zone with shallow groundwater levels. Consider the situation shown in Figure A1, and suppose that the groundwater flow consists of two zones.

**Figure A1.** Schematic flow of groundwater discharge in an evapotranspiration (ET) zone located along a stream.

Zone 1 is the zone below the stream. In this zone, the surface–groundwater exchange rate, *vz*(*x*), is proportional to the head difference between the river, *Hr*, and groundwater, *H*1(*x*), as well as the hydraulic resistance of bottom sediments as follows:

$$v\_{\mathbf{z}}(\mathbf{x}) = \frac{k\_0}{m\_0} (H\_{\mathbf{r}} - H\_1(\mathbf{x})) \tag{A1}$$

where *k*<sup>0</sup> and *m*<sup>0</sup> are the hydraulic conductivity and the thickness of the bottom sediments, respectively, and *x* is the distance perpendicular to the river direction.

In the next zone, zone 2, the rate of groundwater discharge *ETgw*(*x*) can be described as a linear function of the groundwater depth as follows:

$$\left(ET\_{\rm gw}(\mathbf{x})\right) = ET\_{\rm max} \left(\frac{Z\_{\rm surf} - H\_2(\mathbf{x})}{d\_{\rm max}}\right) \; 0 \le Z\_{\rm surf} - H\_2(\mathbf{x}) \le d\_{\rm max} \tag{A2}$$

where *ET*max is the PET and *dmax* is the groundwater ET extinction depth.

Let us consider a steady-state groundwater flow such that the groundwater recharge from the river is balanced by the groundwater ET. Using the Dupuit precondition regarding the head hydrostatic distribution of the aquifer saturated thickness and constant transmissivity, we can express the following system of 1D steady-state equations for the groundwater flow:

$$\begin{split} \frac{d^2 H\_1}{dx^2} - b\_r^2 (H\_1 - H\_r) &= 0, \,\frac{d^2 H\_2}{dx^2} - b\_{ET}^2 (H\_2 - H\_{crit}) = 0, \\\ d\_{\text{max}} &= Z\_{\text{surf}} - Z\_{\text{crit}} \end{split} \tag{A3}$$

where *T* is the transmissivity of the aquifer, and *Zcrit* is the absolute position of the groundwater ET extinction depth.

The boundary conditions for this system assume symmetry of flow in the middle of the river and equality of heads and gradients at the boundaries of the zones as follows:

$$\begin{cases} \mathbf{x} = 0: \frac{dH\_1}{dx} = 0\\ \mathbf{x} = \frac{1}{2}L\_r: H\_1 = H\_2; \frac{dH\_1}{dx} = \frac{dH\_2}{dx} \\ \mathbf{X} \to \infty: H\_2 = Z\_{crit} \end{cases} \tag{A4}$$

The general solution of system (Equation (A3)) is:

$$\begin{aligned} H\_r - H\_1(\mathbf{x}) &= \mathbb{C}\_1 \exp(-b\_r \mathbf{x}) + \mathbb{C}\_2 \exp(b\_r \mathbf{x}) \\ H\_2(\mathbf{x}) - Z\_{\text{crit}} &= \mathbb{C}\_3 \exp(-b\_\varepsilon \mathbf{x}) + \mathbb{C}\_4 \exp(b\_\varepsilon \mathbf{x}) \end{aligned} \tag{A5}$$

The unknown coefficients can be determined using (Equation (A4)) as follows:

$$\begin{aligned} H\_{1}(\mathbf{x}) &= H\_{\text{r}} - (H\_{\text{r}} - Z\_{\text{crit}}) b\_{\text{ET}} \frac{\cosh(b\_{\text{r}} \mathbf{x})}{b\_{\text{ET}} \cosh(0.5 \mathbf{l}\_{\text{r}} b\_{\text{r}}) + b\_{\text{r}} \sinh(0.5 \mathbf{l}\_{\text{r}} b\_{\text{r}})} \\ H\_{2}(\mathbf{x}) &= Z\_{\text{crit}} + (H\_{\text{r}} - Z\_{\text{crit}}) b\_{\text{r}} \frac{\exp(-b\_{\text{ET}} \mathbf{x})}{\exp(-0.5 \mathbf{l}\_{\text{r}} b\_{\text{ET}})} \frac{\sinh(0.5 \mathbf{l}\_{\text{r}} b\_{\text{r}})}{b\_{\text{ET}} \cosh(0.5 \mathbf{l}\_{\text{r}} b\_{\text{r}}) + b\_{\text{r}} \sinh(0.5 \mathbf{l}\_{\text{r}} b\_{\text{r}})} \end{aligned} \tag{A6}$$

Consider the specific flow *q* from the river to the riparian zone as follows:

$$\begin{array}{l} q = -T \frac{dH\_1}{dx} \Big|\_{x=0.5L\_r} = T \frac{H\_r - Z\_{\rm crit}}{\Delta L + L\_{ET}}\\ \Delta L = b\_r^{-1} \coth(0.5L\_r b\_r); \ L\_{ET} = b\_{ET}^{-1} \end{array} \tag{A7}$$

For the more general case of a limited riparian zone with a width of 0.5*L*0, note that (Equation (A7)) is also valid, although *LET* should be calculated as follows:

$$L\_{ET} = b\_{ET}^{-1} \text{coth}(0.5 L\_0 b\_{ET}) \tag{A8}$$

Thus, the groundwater flow from the river used for ET depends on two characteristic lengths. The first length, Δ*L*, characterizes the hydraulic resistance of the bottom sediments. The second length, *LET*, characterizes the hydraulic resistance of the groundwater discharge due to evaporation, which changes linearly with the depth of the groundwater table. Thus, (Equation (A7)) can be rewritten as follows:

$$q = T \frac{H\_r - Z\_{crit}}{L\_{ET}(1 + a\_r)};\\a\_r = \sqrt{\frac{ET\_{\text{max}} m\_0}{d\_{\text{max}} k\_0}} \text{coth}(0.5 L\_r b\_r) \tag{A9}$$

For a small α*<sup>r</sup>* (α*<sup>r</sup>* < 0.1), the hydraulic resistance of the bottom sediments does not play an important role in the river water loss due to ET, while for a large α*<sup>r</sup>* (α*<sup>r</sup>* > 10), the losses of river water due to ET are controlled by the hydraulic resistance of the bottom sediments.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
