Groundwater Depletion in the West Liaohe River Basin, China and Its Implications Revealed by GRACE and In Situ Measurements

Abstract

The West Liaohe River Basin (WLRB) is one of the most sensitive areas to climate change in China and an important grain production base in the Inner Mongolia Autonomous Region of China. Groundwater depletion in this region is becoming a critical issue. Here, we used the Gravity Recovery and Climate Experiment (GRACE) satellite data and in situ well observations to estimate groundwater storage (GWS) variations and discussed the driving factors of GWS changes in the WLRB. GRACE detects a GWS decline rate of −0.92 ± 0.49 km³/yr in the WLRB during 2005–2011, consistent with the estimate from in situ observations (−0.96 ± 0.19 km³/yr). This long-term GWS depletion is attributed to reduced precipitation and extensive groundwater overexploitation in the 2000s. Long-term groundwater level observations and reconstructed total water storage variations since 1980 show favorable agreement with precipitation anomalies at interannual timescales, both of which are significantly influenced by the El Niño-Southern Oscillation (ENSO). Generally, the WLRB receives more/less precipitation during the El Niño/La Niña periods. One of the strongest El Niño events on record in 1997–1998 and a subsequent strong La Niña drastically transform the climate of WLRB into a decade-long drought period, and accelerate the groundwater depletion in the WLRB after 1998. This study demonstrates the significance of integrating satellite observations, ground-based measurements, and climatological data for interpreting regional GWS changes from a long-term perspective.

1. Introduction

Groundwater, as an important component of total terrestrial water storage (TWS), plays a key role in the global water cycle. Overuse of groundwater has caused various environmental problems in many places globally, e.g., desiccation of Lake Urmia, groundwater depletion in Northwest India and in the California Central Valley, and land subsidence in the North China Plain [1,2,3,4,5]. However, monitoring groundwater storage (GWS) is limited in many regions [6]. The absence of in situ monitoring well networks, low-quality observation data, and uncertainties in storage coefficients (i.e., specific yield, S_y) for translating groundwater level changes to groundwater storage variations restrict our knowledge of GWS changes [7,8].

Since its launch in 2002, the Gravity Recovery and Climate Experiment (GRACE) mission has provided a unique way to monitor TWS changes at a monthly scale and with a ~300 km footprint [9,10,11]. Furthermore, it can be used to estimate GWS changes when other TWS components can be removed or neglected. Many studies have compared the GRACE-derived GWS changes with independent in situ groundwater monitoring observations in many parts of the world. At the early stage, Rodell et al. [12] and Strassberg et al. [13] found good agreement between the GWS estimates from GRACE and the Global Land Data Assimilation System (GLDAS) model and those from in situ measurements in the Mississippi River basin and the High Plains Aquifer, USA, respectively. Subsequently, a number of studies have demonstrated the potential of GRACE to estimate GWS depletion in many regional basins, e.g., northwestern India [2,14,15,16,17], California’s Central Valley, USA [3,18,19], the North China Plain (NCP) [20,21,22], and the Middle East [23].

The West Liaohe River Basin (WLRB), which is located in Northeast China (Figure 1) and is known as the “Granary of Inner Mongolia”, is a major grain-producing region in China, with a rapidly growing population since 1950 [24,25]. The basin is in the transition zone between the Inner Mongolia Plateau and the Songliao Plain. The western part of the basin is the south of the Great Khingan Mountains, while the eastern part is the plain area (Figure 1). The WLRB is also the upstream portion of the Liaohe River and overlies a Quaternary loose rock pore aquifer, where extensive groundwater exploitation occurs. Land use in the WLRB is primarily grassland and cropland, with some tree covered land in the mountain region (Figure S1).

The WLRB is located in a semi-arid region, with an area of 136,000 km². It has a continental climate that experiences dry and windy spring conditions, hot and humid summers, cool autumn conditions, and cold winters with little snow. As a semi-arid region, mean annual precipitation of the WLRB is 392 mm (1961 to 2015) and ranges from 260 to 540 mm. Annual precipitation is spatially variable, with less precipitation in the central plain (see Figure S2). About 80% of the precipitation occurs in summer (June to September) [26,27]; while the potential evaporation ranges from 1190 to 1860 mm, which greatly exceeds precipitation [24,28]. Since 1980, the surface flow of the West Liaohe River began to decline, and it nearly dried up since 2000 [24,29,30]. Water storage in lakes and reservoirs also decreased significantly after 1999 [31]. Groundwater has become a vital resource for agricultural, industrial, and domestic water use in the WLRB [32]. Most areas of the WLRB are equipped for irrigation with groundwater [32], which accounted for ~80% of the total water supply during 2005 to 2015. Thus, it is crucial to evaluate the state of groundwater resources in the WLRB for the development of sustainable agriculture and the preservation of terrestrial ecosystems. Groundwater level declines in the WLRB within the past two decades have been reported by several studies [33,34,35,36], e.g., with a rate of ~−0.18 m/yr from 1999 to 2010 [29]. Average groundwater levels in three main counties of the WLRB in 2006 decreased by 3.8 m, 2.8 m, and 1.9 m, respectively, when compared to those in 1978 [33]. However, quantitative evaluation of the GWS changes in the WLRB for the past decade was still limited.

In this study, we use GRACE, GLDAS land surface models (LSMs), and in situ measurements to evaluate the GWS changes during the GRACE era in the WLRB. Recent studies indicate that there has been an abrupt change in the annual precipitation in the WLRB since 1998, i.e., before the launch of GRACE satellites [37]. Therefore, longer period groundwater storage changes should be evaluated and interpreted. Based on the GRACE data and precipitation observations, we further reconstructed long-term TWS variations in the region prior to the GRACE era. The objective of this study was to address the following questions: (1) How can groundwater monitoring in the WLRB benefit from GRACE satellite data? (2) How do water storages change in response to human demand and climate forcing in the WLRB?

2. Data and Methods

TWS refers to all of the water stored in the land, including surface water, snow and ice, soil moisture, and groundwater. GWS anomalies (GWSA) can be estimated by subtracting surface reservoir storage anomalies (RESSA), snow water equivalent storage anomalies (SWESA), and soil moisture storage anomalies (SMSA) from the GRACE-derived TWS anomalies (TWSA):

$$\mathrm{GWSA}=\mathrm{TWSA}-\mathrm{RESSA}-\mathrm{SWESA}-\mathrm{SMSA}$$

(1)

where RESSA can be estimated from in situ monitoring and is negligible in some cases; and, SWESA and SMSA can be estimated by hydrological models or in situ observations. In this study, in situ groundwater level anomalies (GWLA) were used to validate the GWSA derived from GRACE and ancillary data. In addition, precipitation data and climate indices were used to assess the impact of climate variability on TWS/GWS changes. A summary of the datasets that were used in this study can be found in

Table 1. Summary of the datasets used in this study.

Type	Data Name	Data Version/Source	Spatial and Temporal Resolution, Ttime Period	Data Access
Satellite	TWSA	CSR Mascons	1° × 1°, monthly, 2005–2015	http://www2.csr.utexas.edu/grace/RL05_mascons.html
Models	SMSA+SWESA	GLDAS-1 LSMs	1° × 1°, monthly, 2005–2015	http://mirador.gsfc.nasa.gov/
In situ	GWSA	CIGEM+HBIM	122 wells, monthly, 2005–2015	Groundwater yearbook and Hydrological Bureau of Inner Mongolia
	GWLA		Region mean, annual, 1980–2015	http://116.113.33.52:8989/ and Xu [29]
	RESSA		Region sum, annual, 1997–2015	http://116.113.33.52:8989/
	Precipitation	CMDC	0.5° × 0.5°, monthly, 1961–2015	http://data.cma.cn/
			0.5° × 0.5°, daily, 1961–2013	http://data.cma.cn/
Climate indices	ONI		Time series, monthly, 1980–2015	http://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php

.

2.1. TWSA Derived from GRACE

We used the GRACE Mascons solution that was provided by the Center for Space Research (CSR), the University of Texas at Austin, to derive TWSA [38]. Both mascon and CSR RL05 spherical harmonic products were based on the same GRACE level-1 observations, also with C20 replacement, degree 1 corrections, and glacial isostatic adjustment (GIA) corrections. CSR Mascons (CSR-M) are provided with a 0.5 degree spatial resolution. The regularization constraint on mascon solutions is derived from original GRACE information with no empirical filtering post-processing [38]. Some studies have concluded that, although no empirical post-processing is required, mascon solutions can still provide comparable TWSA estimates as traditional spherical harmonic products [38,39].

2.2. SMSA+SWESA Derived from GLDAS LSMs

We used four versions of LSMs from GLDAS to estimate SMSA and SWESA, i.e., Noah [40], CLM [41], MOSAIC [42], and VIC [43]. The GLDAS LSMs were developed by the National Aeronautics and Space Administration (NASA) and the National Oceanic and Atmospheric Administration (NOAA). Constrained by ground- and space-based observations, the GLDAS LSMs can simulate the optimal fields of land surface states and fluxes [44]. Simulated SMSA from Noah, MOSAIC, VIC, and CLM include changes in soil water contents with depths of 2.0 m, 3.4 m, 3.5 m, and 1.9 m, respectively.

2.3. GWSA from In Situ Observations

The groundwater monitoring well observations were obtained from the groundwater level yearbook that was published by the China Institute of Geological Environment Monitoring (CIGEM) [45] and from the Hydrological Bureau of Inner Mongolia (HBIM) for the period of 2005–2015. The mean groundwater level anomalies (GWLA) from a total of 122 observation wells and the GWLA from each well were shown in Figure S3. Some of the raw groundwater level data are highly inaccurate with variability between consecutive months greater than 2 m in some records, which were eliminated during the preprocessing stage. The month-to-month number of available well records varied from 45 to 103. The GWSA from ground well observations are calculated as follows [7,46]:

$$GWSA={{\displaystyle \sum }}_j^N{S}_j{W}_j∆h_j/{{\displaystyle \sum }}_j^N{W}_j$$

(2)

where S_j is the S_y value of each grid cell j; N represents the total number of grid cells; W_j is the weight of each grid cell (the cosine of the latitude of the corresponding cell); and, ∆h_j refers to the mean value of the groundwater level variations in each grid cell. We divided the study region into a 1° × 1° mesh as Strassberg et al. [13]. In this study, a mean S_y value of 0.042, determined by pumping tests, was used to translate groundwater level changes to water storage changes, as suggested by the Songliao Water Resources Commission (SLWRC) [47].

Besides the monthly groundwater level observations from CIGEM and HBIM, we also collected annual mean groundwater levels in the WLRB for the period of 1980–2010 from Xu [29] and for the period of 1997–2015 from the annual Water Resources Bulletin (WRB) released by the HBIM [32]. Then, we derived a long-term time series of annual GWLA from 1980 to 2015. Note that we used the mean values of annual groundwater levels from the two sources during the overlap period (1997–2010).

2.4. RESSA from Ground Measurements

To estimate GWSA in the WLRB from GRACE-derived TWSA, RESSA, as a potentially important TWS component in some regions, should be taken into account. For the WLRB, annual mean storage data of surface water were collected from the WRB, as published by the HBIM [32]. The RESSA trend was −0.05 km³/yr from 2005 to 2011, and −0.02 km³/yr from 2005 to 2015. Although the contribution of RESSA is rather small, we still removed it from the GRACE-derived TWSA to obtain GWSA.

2.5. Precipitation Data and Climate Indices

Gridded precipitation data over the study region were obtained from the China Meteorological Data Service Center (CMDC). These gridded data were obtained by spatial interpolation of precipitation observations from 2472 weather stations with monthly and daily temporal resolutions and 0.5° × 0.5° spatial resolution over all of China. Additionally, we used NOAA’s Oceanic Niño Index (ONI) to represent the El Niño-Southern Oscillation (ENSO) variability. Positive ONI values represent the warm ENSO phase (i.e., El Niño), while negative ONI values represent the cold ENSO phase (i.e., La Niña).

2.6. Climate-Driven Water Storage Variability

To further assess the impact of climate variability on TWS and better interpret GWS changes in the WLRB, a long-term water storage variation time series from 1980 to 2013 was reconstructed based on long-term precipitation data and eleven-year GRACE CSR-M data. The details on the reconstruction method of climate-driven sub-decadal water storage variability can be found in [48,49]. In brief, first, the effect of one day’s precipitation flux to later days can be expressed by an exponential decay filter:

$$f\left(t\right)=\{\begin{array}{c}0,t<0\\ \frac{e^{-\frac{t}{\tau }}}{\tau },t\ge 0\end{array}$$

(3)

where t is time and $\tau$ is a decay parameter which is calibrated. Then, filtered daily time series was averaged to the monthly resolution using an arithmetic mean. As the exponential decay filtered daily precipitation is highly correlated with interannual and subseasonal TWSA [49], the precipitation-driven TWSA was formulated as following:

$${\mathrm{TWSA}}_{rec}=\beta ·P_{inter+subseas}^{\tau }$$

(4)

where ${\mathrm{TWSA}}_{rec}$ is the reconstructed nonseasonal and detrended TWS variations; $P_{inter+subseas}^{\tau }$ corresponds to the deseasoned and detrended filtered precipitation. The parameters $\tau$ and $\beta$ were calibrated by using daily precipitation from CMDC and ${\mathrm{TWSA}}_{inter+subseas}$ derived from deseasoned and detrended GRACE-derived TWSA. Because the precipitation data from CMDC with daily resolution are available until 2013, we reconstructed the TWSA time series from 1980 to 2013. In the original TWSA reconstruction by Humphrey et al. [49], temperature change was also considered in addition to precipitation. In this study, however, we only used precipitation data because temperature has less effect on TWSA (Vincent Humphrey, personal communication).

2.7. Regression Analysis between Groundwater Flux and Precipitation Anomalies

To investigate the effect of precipitation anomalies on groundwater storage changes, we constructed a linear relationship between the variations in precipitation and GWS, as follows [50]:

$$F=a\ast ∆P+b$$

(5)

where F is annual groundwater flux (GWF, GWS change between two successive years, in km³/yr); ∆P is annual precipitation anomalies relative to the study period, which is multiplied by the area of WLRB to convert it to km³/yr. In addition, a and b are the parameters that are calibrated by regression analysis: a is the ratio of groundwater flux to precipitation anomalies (dimensionless); b is the residual that cannot be explained by precipitation anomalies, which can be attributed to anthropogenic effects.

2.8. Uncertainty Assessment

The TWSA uncertainty was estimated based on the approach suggested by Landerer and Swenson [51] and Scanlon et al. [39]. We first obtained the residuals by removing the linear trend and seasonal components from TWSA time series. The interannual signals were further removed by fitting a 13-month moving average to the residuals because they also contributed significantly to the TWSA in the WLRB. The RMS of the residuals was used as the measurement uncertainty in TWSA. Our uncertainty estimate was relatively conservative, as the residuals would still contain subseasonal signals. The uncertainty in SMSA+SWESA was estimated from the standard deviation (STD) among the time series that were produced by the four GLDAS models.

We computed the trend in GWSA based on GRACE-derived TWSA and mean SMSA+SWESA from the four GLDAS models (Noah, CLM, Mosaic and VIC). The GWS trend estimates were calculated using unweighted least squares to fit a linear trend, plus annual and semiannual sinusoids to the time series of GWSA [15]. The uncertainty in GWS trend was also estimated when considering the error propagation law during the least squares fit and uncertainties in TWSA and SMSA+SWESA. Uncertainty in S_y was estimated to be ±20%, which was used to estimate the uncertainty in the GWS trend from in situ observations.

3. Results

3.1. GWSA from 2005 to 2015

Figure 2a shows the time series of TWSA estimated from GRACE, and mean SMSA+SWESA from four GLDAS LSMs. The shadows show the uncertainties in TWSA (±1.49 cm) and SMSA+SWESA (±0.72 cm). TWSA uncertainty is close to the global mean uncertainty value of TWSA (~2 cm), as suggested by CSR. As shown in Figure 2a, TWSA decreases at a rate of −1.26 ± 0.49 cm/yr from 2005 to 2009 continuously, and then shows relatively large interannual variability from 2010 to 2013, and decreases rapidly again after 2013. The SMSA+SWESA show comparable interannual changes as TWSA from 2005 to 2015. Both TWSA and SMSA+SWESA show two peaks at the end of 2010 and 2012.

Figure 2b shows the time series of GWSA based on GRACE and in situ measurements. During the drought period from 2005 to 2011, the GRACE-derived GWSA in the WLRB declined rapidly at a rate of −0.92 ± 0.49 km³/yr (i.e., −0.68 ± 0.36 cm/yr in equivalent water height (EWH)). For the same period, the GWSA trend from in situ measurements is −0.96 ± 0.19 km³/yr (i.e., −0.70 ± 0.14 cm/yr EWH), which is in good agreement with the GRACE-derived estimate. The time series from GRACE and in situ observations also compare favorably before 2011.

For the whole study period of 2005–2015, GRACE-derived GWSA in the WLRB, as declined at a rate of −0.43 ± 0.26 km³/yr (i.e., −0.32 ± 0.18 cm/yr in EWH), significantly less than the rate from 2005 to 2011. The obvious interannual signals during 2012–2015 result in the difference of GWS rates between the two periods. The GWS rate from well observations is −0.73 ± 0.15 km³/yr from 2005 to 2015.

Figure 2c shows annual precipitation anomalies, annual GWSA results from GRACE, in situ measurements, and WRB. In addition to the obvious long-term GWS depletion for the entire time period, the GRACE-derived GWSA increased in 2012 and 2013 in response to more precipitation in 2012, which is also shown by the GWSA from WRB, but with relatively lower magnitude. The continuous GWS decrease that as observed by in situ measurements also slows down from 2011 to 2013. For the years of 2014 and 2015, the GRACE-derived GWSA decreased rapidly again.

3.2. GWSA from 1980 to 2015

Both GRACE and in situ measurements show significant groundwater depletion in the WLRB from 2005 to 2011. Over a longer time period, Figure 3a shows the mean annual groundwater table changes in the WLRB which declined slowly from 1980 to 1998, and declined more rapidly after 1998. These data indicate that groundwater levels in the WLRB had decreased before the GRACE era. As shown in Figure 4a, the WLRB experienced a long-term drying climate and had the lowest average annual precipitation since 1980 for the period of 1999–2009 [37]. The average annual precipitation during this period decreased by 15% relative to the average annual precipitation from 1961 to 2015. The turning point in precipitation anomalies in 1998 is consistent with the occurrence of one of the largest El Niño events on record in 1997–1998, followed by the long-term drought (Figure 4a).

Although the WLRB did not experience significant droughts before 1998, annual groundwater exploitation increased continuously due to extension of crop cultivation and expanding irrigation (Figure 3b) [52]. Annual groundwater exploitation also shows a jump in 1998, and remains high since then. This feature also agrees well with the occurrence of the 1997–1998 El Niño event and the following long-term precipitation deficit. This jump in the signal is consistent with the acceleration in groundwater level decline after 1998 (Figure 3a). Although the annual water consumption statistics (including agricultural, industrial, and urban use) in the WLRB are only available after 1998, we find that ~80% of water consumption in this region is from groundwater exploitation during 1998–2015 (Figure 3c), consistent with the continuous groundwater level decline (Figure 3a).

Relative to the average of precipitation from 1980 to 2015, precipitation anomalies are positive in most of years prior to 1998, but negative for the continuous eleven years from 1999 to 2009 (Figure 4a). Time series of ONI highlights the strong El Niño/La Niña years from 1980 to 2015 (Figure 4). Precipitation anomalies in the WLRB agree well with ENSO events (Figure 4a). For example, annual precipitation increases significantly during the three El Niño periods (e.g., 1987–1988, 1991–1992, and 1997–1998). In contrast, annual precipitation decreases generally during the La Niña periods (e.g., 1988–1989, 1998–2000, 2007–2008, and 2010–2011).

The reconstructed climate-driven TWSA and detrended mean GWLA in the WLRB show similar interannual changes during 1980–2013 (Figure 4b). Both ONI and reconstructed TWSA were applied with a Vondrak filter to exclude the high frequency signal of less than one year [53,54]. Generally, both show significant increases during El Niño events and decreases during La Niña events. TWSA and GWLA reach their maxima in 1998, consistent with the extreme precipitation increase in this year due to the strong 1997–1998 El Niño event. After the abrupt precipitation increase in 1998, GWLA and TWSA decreased from 1999 to 2010, consistent with the long-term precipitation deficit associated with continuous La Niña periods.

3.3. Groundwater Flux (GWF) from 1980 to 2015

The linear relationship between GWF and precipitation anomalies were constructed using the regression analysis method in Section 2.7. The difference of groundwater levels between two successive years was multiplied with the mean S_y to retrieve the time series of mean GWF over the WLRB. Note that the negative GWF values represent GWS depletion when compared to the previous year; and vice versa. More negative GWF values after 1998 indicate a long-term GWS depletion in the WLRB (Figure 5a), consistent with the results in Figure 3a. Due to distinct features for GWLA before and after 1998, we split the time period into two parts, i.e., 1980–1998 and 1998–2015, and constructed a linear relationship for these two periods.

For the first time period, the GWF is more sensitive to precipitation anomalies. The ratio parameter a in Equation (5) for the period 1980–1998 is 0.088. It indicates that 1 km³ (i.e., ~7.35 mm over the WLRB) of precipitation anomalies can cause a groundwater storage change of 0.088 km³. In other words, ~9% of annual precipitation recharges the aquifers in the WLRB. This ratio is highly consistent with the mean precipitation infiltration coefficient from 1988 to 1997 (i.e., 0.087) in the WLRB estimated by Xu [29]. This good agreement confirms that the linear relationship in this time period is reasonable. Based on Equation (5), we assume that the anthropogenic effect is constant and can be estimated by the parameter b. The anthropogenic contribution to GWF (i.e., discharge of aquifer) is estimated to be −0.501 km³/yr. The coefficient of determination R² between GWF and regression-modeled GWF reaches 0.64; while the RMS between them is 0.601 km³/yr. Good agreement between the GWF and regression-modeled GWF also indicates that it is reasonable to assume a linear relationship between GWF and precipitation anomalies.

For the second time period, however, precipitation anomalies contribute less to GWF, with a ratio of 0.071. The R² between GWF and regression-modeled GWF is only 0.22, which indicates that the linear relationship between the GWF and precipitation is relatively weak for this time period. Figure 5a shows poor agreement between GWF and regression-modeled GWF for the period of 1998–2015.

The discrepancy from 1998 to 2015 indicates that there are still some limitations in this linear regression analysis method. First, we do not consider the time lag of groundwater change in response to precipitation, which could take several months or even more than one year. Because the groundwater table in the WLRB was relatively shallow during 1980–1998 (i.e., 2–3 m, shown in Figure 3a), the effect of the delayed response of groundwater is limited. From 1998 to 2015, however, the groundwater table reaches a depth of 5–6 m. Thus, the time lag might be non-negligible, even on interannual time scales. This could be one of the reasons to explain the discrepancy during 1998–2015. Second, in this method, we assume a constant anthropogenic contribution to groundwater, i.e., a fixed annual groundwater exploitation amount. Under the long-term drying climate in the WLRB since 1998, however, the annual groundwater exploitation amount probably changes over time. It can also explain the poor agreement between GWF and regression-modeled GWF from 1998 to 2015.

Long-term GWSA was compared with the modeled GWSA (Figure 5b). Note that the integral of the regression-modeled GWF with respect to time was computed for the aforementioned two time periods to derive the modeled GWSA from 1980 to 2015. The two time series are highly correlated (R² = 0.99) (Figure 5b), which indicates that the linear regression analysis method can be used to model the GWSA very well. The trends in GWSA and modeled GWSA from 1980 to 1998 are −0.10 km³/yr and −0.06 km³/yr, respectively; while the trends during the second time period are −0.87 km³/yr and −1.04 km³/yr, respectively. Similar to the results in Figure 5a, the time series of GWSA and modeled GWSA agree with each other better for the period of 1980–1998, with consistent interannual fluctuations. From 1998 to 2015, however, some interannual signals in GWSA, e.g., in the 2000s, cannot be modeled well due to the limitations in the linear regression method that are mentioned before.

4. Discussion

4.1. Interannual Variablity of GWS from 2012 to 2015

GRACE-derived GWSA shows a significant increase in GWS from 2012 to 2013 followed by a rapid decline in 2014 and 2015. This interannual variability cannot be found from in situ well observations (Figure 2b). The GWS increase in 2012 and 2013 can be attributed to the abnormal precipitation increase in 2012 (Figure 2c). The statistics from the WRB also show a similar GWS increase in 2012 and 2013 as GRACE, but with lower amplitude (Figure 2c). One possible reason for missing interannual variability in ground observations is the uneven distribution of monitoring wells. As shown in Figure 1, there are few wells in the mountain area and the upstream area of the plain. In addition, in this work, we used a mean S_y (i.e., 0.042) from the SLWRC [47] and failed to derive the spatial information of S_y in the WLRB. It would also cause uncertainties during the estimation of in situ GWSA.

There are uncertainties in GRACE-derived GWSA too. The area of the WLRB is slightly less than the GRACE nominal resolution, i.e, ~200,000 km². But, some studies have concluded that GRACE is able to detect water storage changes in localized regions with areas below GRACE resolution if the signal magnitude is sufficiently large [56,57]. In addition, simulated SMSA+SWESA from global GLDAS LSMs at small spatial scales, such as the WLRB, is likely to be underestimated or overestimated. Nevertheless, interannual fluctuations in GRACE-derived GWSA from 2012 to 2015 are still statistically significant (Figure 2b,c). The time series of GWSA based on other GRACE products and methods also show similar interannual variability from 2012 to 2015 (not shown), which confirm this interannual signal that is shown in Figure 2b,c.

4.2. Uncertainties in GRACE-Derived GWS Rates

Based on the law of error propagation and the uncertainties in GRACE-derived TWSA and model-based SMSA+SWESA, the uncertainty in GRACE-based GWS rate in the WLRB is 0.49 km³/yr and 0.26 km³/yr, respectively, for the time period 2005–2011 and 2005–2015. As the area of the LRB is less than the GRACE nominal spatial resolution, we further compared the results based on various GRACE solutions and methods to validate the GWS depletion signal in the WLRB.

As shown in

Table 2. GWS rates in the WLRB estimated from different GRACE solutions. Center for Space Research (CSR)-M represents the CSR Mascons used in this study [38]; i.e., JPL-M.dsf represents the Jet Propulsion Laboratory (JPL) Global Mascons [58] with grid scaling factors applied; GSFC-M represents the Goddard Space Flight Center (GSFC) Global Mascons [59]; CSRT-GSH represents GRCTellus Land data [51]; while, CSR-DDK3 represents CSR spherical harmonic solutions with DDK3 filtered [60]. The uncertainties were given as the standard deviations, which means that the error bounds represent the 68% confidence interval.

GRACE Solutions	Rate (km3/yr)	Rate (km3/yr)
	January 2005–December 2011	January 2005–December 2015
CSR-M	−0.92 ± 0.49	−0.43 ± 0.26
JPL-M.dsf	−1.54 ± 0.42	−0.45 ± 0.21
GSFC-M	−1.47 ± 0.46	−0.62 ± 0.23
CSRT-GSH	−0.84 ± 0.47	−0.51 ± 0.24
CSR-DDK3	−1.67 ± 0.65	−0.65 ± 0.33

, the GWS rate estimates during 2005–2015 derived from different GRACE solutions varied from −0.43 km³/yr to −0.65 km³/yr, agreeing with each other with respect to the error bounds. Additionally, all of the GRACE solutions show a higher GWS depletion rate during 2005–2011 than that for the entire period 2005–2015. For the period 2005–2011, GWS rates from various GRACE solutions show some variability, ranging from −0.84 km³/yr to −1.67 km³/yr, but still agree with each other within the uncertainty envelope. The GWS rate estimated from GRACE CSR Mascons and used in this study is −0.92 ± 0.49 km³/yr during 2005–2011, most consistent with the estimate from in situ measurements (−0.96 ± 0.19 km³/yr). The comparison of various GRACE solutions confirmed the significant depletion of GWS occurred in the WLRB.

4.3. Comparison with Other Studies

Xu [29] constructed a distributed hydrological model for basin scale water cycle simulation (MODCYCLE) over the WLRB for 1988–1997 and 2001–2010. The input data contained basic geographic information data, meteorological observations (e.g., precipitation, air temperature, wind speed, etc.), soil parameters, groundwater level data, and water supply and use data, etc. Based on the MODCYCLE model over the WLRB, the mean recharge to the unconfined aquifers was 3.94 km³/yr, and the mean discharge from the unconfined aquifers was 5.17 km³/yr from 2001 to 2010. Thus, the GWS depletion rate for 2001–2010 is −1.23 km³/yr, which agrees well with our GRACE estimate (i.e., −0.92 ± 0.49 km³/yr) from 2005 to 2011. Li et al. [24] modeled the groundwater balance in the plain area of the WLRB from July 2003 to June 2004 based on the Visual MODFLOW software, with groundwater table, precipitation infiltration amount, and groundwater exploitation amount data as inputs. Based on their model, GWS in the unconfined aquifers declined at a rate of −1.23 km³/yr from July 2003 to June 2004, while the confined aquifers were in an equilibrium state. Their estimate is also consistent with the results in this study and from Xu [29].

Besides monthly groundwater level observations, we also collected annual groundwater levels from the WRB from 2005 to 2015 [32]. GWS depletion rate in the WLRB from these annual well observations is −0.78 km³/yr for 2005–2011 and −0.75 km³/yr for 2005–2015, consistent with the results from monthly level observations and GRACE.

A recent study by Zhao et al. [36] only used groundwater level data from the CIGEM to estimate GWS changes in the WLRB plain and found a rate of −2.76 ± 0.15 cm/yr from 2005 to 2009, which significantly overestimated the actual GWS depletion rate when compared to results from our study and the aforementioned previous studies. As shown in Figure S4, the distribution of groundwater wells from the HBIM is relatively uniform; while the monitoring wells from the CIGEM are mainly located in the two cones of groundwater depression centered within the WLRB, explaining the likely overestimation of GWS depletion. In this paper, however, combining the well data from two sources provides a more accurate estimate of GWSA.

Several studies have revealed the long-term GWS depletion in the North China Plain (NCP), which is located in the ~600 km southwest of the WLRB [20,21,22]. For example, Feng et al. [20] found a GWS depletion rate of 8.3 ± 1.1 km³/yr from 2003 to 2010 in the NCP. Although the GWS decline signal in the WLRB is less than one-eighth of that in the NCP, GRACE still detects it successfully, which is confirmed by in situ observations.

4.4. ENSO’s Effect on Groundwater from a Long-Term Perspective

It is well known that ENSO is an important factor that is influencing the precipitation in some river basins in China, such as the Yellow River basin, the Yangtze River basin, and the Pearl River basin [61,62]. In this study, the significant influence of ENSO on precipitation is also found in the WLRB. In general, precipitation anomalies are positive during the warm phase of ENSO (i.e., El Niño), such as 1987–1988, 1991–1992, and 1997–1998; while, the WLRB experiences the dry conditions during the La Niña events, such as 1988–1989, 1998–2000, 2007–2008, and 2010–2011. Long-term time series of GWLA show a slow decline prior to 1998, followed by a significant rapid decline after that (Figure 3a). The remarkably different behaviors of GWLA before and after 1998 can be attributed to one of the strongest El Niño events on record in 1998 and the following continuous La Niña periods. Detrended time series of GWLA also shows the interannual variability that is related to ENSO (Figure 4b), i.e., GWS increase with El Niño events and decrease with La Niña events.

4.5. Groundwater Recharge and Discharge

According to the water balance, groundwater storage changes in aquifers represent the difference between groundwater recharge and discharge. For the WLRB, to further identify the driver of groundwater depletion, the recharge and discharge status of the aquifers over the WLRB was analyzed. The unconfined aquifer recharge in the WLRB is mainly from the infiltration of precipitation (precipitation infiltration recharge, PIR), which accounts for ~75% of total recharge during normal years, ~57% during dry years, followed by surface water leakage recharge and seepage of irrigation water [29]. Based on the mean atmospheric precipitation infiltration coefficients during the wet (1988–1997) and dry periods (2001–2010), as estimated from Xu [29], the mean PIR decreased significantly from 5.85 km³/yr in the wet years (1984–1998) to 2.35 km³/yr during the dry years (1999–2015). The RESSA also declined rapidly since 1999 in the WLRB [32]. For example, the RESS in the WLRB plain was only 0.13 km³ at the end of 2004, accounting for ~10% of the total design capacity, and most of the reservoirs approached to their dead storage [24]. Therefore, the mean surface water leakage recharge in the WLRB would also decrease due to the reduced surface water storage.

Groundwater discharge is dominated by groundwater exploitation and evaporation from water-table (phreatic) aquifers in the WLRB. These two mechanisms of groundwater discharge account for more than 90% of the total groundwater discharge [24]. As the surface water was largely utilized in the upstream of the WLRB, the source of agricultural irrigation water in the plain region changed from surface water to groundwater since the 1980s [24,63]. Figure 3c shows the annual surface water use and groundwater exploitation in the WLRB since 1998. In the WLRB, mean surface water use was only 0.90 km³/yr from 1998 to 2015, but mean groundwater exploitation was nearly 3.70 km³/yr. The annual groundwater withdrawal for agriculture in the WLRB plain region increased steadily from 0.89 km³ to 1.87 km³ before 1998, but dramatically rose to 2.67 km³ in 1999, and remained at this level after that (Figure 3c). Note that agricultural groundwater exploitation accounts for ~82% of the total groundwater exploitation in the WLRB [29,47]. The increasing large percentage of groundwater use in total water use over the WLRB is consistent with the accelerated groundwater table decline and relatively less precipitation since 1999 (Figure 3a).

Based on the statistics available over the WLRB plain region (Table S1), although the annual precipitation decreased from 426 mm in the 1990s to 323 mm in the 2000s, groundwater irrigation area expanded from 3406 km² to 5641 km²; meanwhile, crop yield increased from ~3.1 × 10⁶ tons to ~3.8 × 10⁶ tons. Increasing groundwater irrigated area and crop yield, but decreasing precipitation also indicates the serious groundwater overdraft in the WLRB in the 2000s, as consistent with accelerated groundwater level/storage declines after 1998 (Figure 3a and Figure 5b).

5. Conclusions

GWSA in the WLRB was estimated using GRACE and in situ well observations from 2005 to 2015. Significant GWS depletion and interannual variability in the WLRB were detected. GRACE detects a GWS depletion rate of −0.92 ± 0.49 km³/yr, which is equivalent to a total loss of 6.44 km³, during the dry period of 2005–2011. From 2012 to 2015, however, GRACE-derived GWSA shows a prominent inverse V-shape interannual change due to an abnormal precipitation increase in 2012. For the entire 2005–2015 time period, the GWS depletion rate in the WLRB from GRACE is −0.43 ± 0.25 km³/yr. In situ well observations, however, indicate the steady GWS depletion for the entire time period, with a rate of −0.96 ± 0.19 km³/yr from 2005 to 2011, and −0.73 ± 0.15 km³/yr from 2005 to 2015. A remarkable precipitation deficit during the 2000s, associated with La Niña, and continuous groundwater overdraft result in the long-term GWS depletion in the WLRB from 2005 to 2015.

Over a period of almost four decades, long-term interannual TWSA and GWSA in the WLRB are highly correlated with precipitation anomalies, which are driven by ENSO. In general, WLRB receives more precipitation during El Niño and less precipitation during La Niña periods. Annual precipitation is relatively stable during 1980–1997. However, due to the strong 1997–1998 El Niño and the following continuous La Niña, an 11-year drying period after 1998 and the consequent extensive groundwater overexploitation result in significant groundwater depletion in the WLRB. In other words, due to precipitation deficits in the 2000s, groundwater depletion in the WLRB is attributed to both the reduced groundwater recharge from precipitation infiltration and the increasing groundwater discharge due to irrigation.

While GRACE provides a powerful tool to monitor regional GWS variations, there are still large uncertainties in the GRACE-derived GWSA estimates, which mainly result from the GRACE measurement and post-processing errors and the uncertainty in estimated SMSA+SWESA from LSMs. The launch of GRACE Follow-On in early 2018 should extend the TWS observations. The Next Generation Gravity Mission (NGGM) is able to provide higher spatial resolution gravity fields, and will be of great benefit to hydrologic applications [64,65].

In addition, the ongoing development of the new national groundwater monitoring project in China will provide more accurate estimates of GWSA and improved spatial coverage (http.www.mlr.gov.cn/). More space- and ground-based observations with better spatiotemporal resolution and coverage in the near future will enhance our understanding of groundwater changes at the global and regional scales, and guide water resource management for policy makers.