Human-Induced and Climate-Driven Contributions to Water Storage Variations in the Haihe River Basin, China

Abstract

Terrestrial water storage (TWS) can be influenced by both climate change and anthropogenic activities. While the Gravity Recovery and Climate Experiment (GRACE) satellites have provided a global view on long-term trends in TWS, our ability to disentangle human impacts from natural climate variability remains limited. Here we present a quantitative method to isolate these two contributions with reconstructed climate-driven TWS anomalies (TWSA) based on long-term precipitation data. Using the Haihe River Basin (HRB) as a case study, we find a higher human-induced water depletion rate (−12.87 ± 1.07 mm/yr) compared to the original negative trend observed by GRACE alone for the period of 2003–2013, accounting for a positive climate-driven TWSA trend (+4.31 ± 0.72 mm/yr). We show that previous approaches (e.g., relying on land surface models) provide lower estimates of the climate-driven trend, and thus likely underestimated the human-induced trend. The isolation method presented in this study will help to interpret observed long-term TWS changes and assess regional anthropogenic impacts on water resources.

1. Introduction

Terrestrial water storage (TWS) is the sum of all forms of water storage over land surfaces, which plays a key role in the global and regional hydrological cycle [1,2,3]. TWS varies on a wide range of temporal and spatial scales, and exerts important effects upon the Earth’s climate system [4,5,6]. After the launch of the Gravity Recovery and Climate Experiment (GRACE) satellites in 2002, many studies have quantitatively estimated the TWS or groundwater storage (GWS) changes using GRACE data in global and regional basins [7,8,9,10]. GRACE satellites can derive global monthly temporal gravity fields with a spatial resolution of ~300 km, which can be further used to estimate water storage changes over land when other mass change effects are removed from the gravity fields [11,12]. However, the lack of in situ observations limits our knowledge of the relationship between natural climate variability, and human activity impacts upon terrestrial water systems [13,14,15,16].

Water resources in basins are mainly influenced by climate changes and anthropogenic activities [14]. Assessing changes in water storage on the global scale under the impact of climate variability and human activity has drawn great attention to the growing water scarcity around the world [15,17,18,19,20,21,22,23].

On one hand, the climatic fluctuations in ocean and atmosphere, such as El Niño–Southern Oscillation (ENSO), resulting from the large-scale ocean–atmosphere interactions over the equatorial Pacific [24], can strongly affect global water storage by influencing precipitation patterns [25,26,27]. During the warm phase of ENSO (i.e., El Niño), heavy rains generally occur in South America, while droughts occur in Southeast Asia and Australia. On the contrary, during the cold phase of ENSO (i.e., La Niña), heavy rains occur in Southeast Asia and Australia, while droughts occur in Southwestern USA and the West Coast of South America. Other climatic fluctuations, e.g., Pacific Decadal Oscillation (PDO), Indian Ocean Dipole (IOD) and North Atlantic Oscillation (NAO) also impact TWS anomalies (TWSA) at different temporal and spatial scales [28,29]. On the other hand, due to the increasing water resource management and consumption all around the world, impacts of anthropogenic activities on water storage become more important, and should be taken into account in the hydrologic cycle and modeling [30,31]. Impacts of direct human activities could be through water withdrawal [32,33,34], the building of reservoirs [35,36], or land use and land cover change [37].

Some previous studies have focused on analyzing the relative contribution of TWSA from natural and anthropogenic effects [13,14,17,38,39]. For instance, Yi et al. [39] estimated an anthropogenic TWS depletion estimate of −187 ± 38 km³/yr in Asia by using a linear relationship between variations in precipitation anomalies and water storage fluxes. Huang et al. [14] isolated the human-induced TWSA by subtracting climate-driven TWSA simulated by models from GRACE-observed TWSA. As GRACE can detect the total TWS changes affected by both climate variability and human activities, whereas land surface model (LSM) simulations generally represent the climate-related TWS changes in most cases. Felfelani et al. [13] quantified the human-induced TWSA by calculating the difference between the two hydrological models, in which human impacts were considered in one model, but not included in the other one. However, there are some limitations and uncertainties in terms of modeling and forcing data. For example, Yi et al. [39] assumed that the anthropogenic component is constant over time, and the climate-driven effect is zero from 1979 to 2015, and they analyzed both impacts on the annual scale. A systematic difference between various precipitation-forcing data was also highlighted [39]. Felfelani et al. [13] also found that the TWSA trend uncertainty caused by the errors in forcing datasets was as high as the TWSA trend differences resulting from different hydrological models. Moreover, global hydrological models and land surface models may exhibit large biases in TWSA amplitude [8,20].

Here we quantify the contribution of climate-driven TWSA and its trend based on the reconstruction of long-term climate-driven water storage variations using one decade of GRACE observations and a more than four-decade precipitation dataset, and further evaluate the human-induced TWSA by removing the climate-driven TWSA from the GRACE-derived TWSA. We employ the information from the current GRACE record and avoid the difficulties of using hydrological models by using a statistical approach. The upper and lower bounds of the climate-driven TWSA trend (Trend_CD) and human-induced TWSA trend (Trend_HI) are estimated based on a relatively long time period of reconstructed, climate-driven TWSA time series.

2. Materials and Methods

2.1. Study Area

The Haihe River Basin (HRB) is selected as a case study, which is located in North China (~320,000 km²; Figure 1a). The HRB has a total population of ~124 million, including Beijing, Tianjin, and most of the area of Hebei, which is the political, economic and cultural center of China (Figure A1). As a major grain-producing region in China, it produces 10% of China’s total grains [40]. The major crops in this area are winter wheat and summer maize.

In this basin, the mean annual precipitation is ~530 mm/yr from 1961 to 2013. The runoff is mostly generated in the mountainous regions; its mean runoff in the hydrologic control station of Haihezha is only 0.82 km³/yr from 1960 to 2010 [41]. The water consumption is concentrated on the plain region, and it contains the main agricultural areas, industries and cities.

Agriculture depends largely on groundwater exploitation in the HRB, which is similar to the High Plains in the USA and northwest India, and it has become the major groundwater depletion areas around the world [42]. The GWS suffers a prolonged declining trend of −17.8 ± 0.1 mm/yr during 1971–2015 in the North China Plain [43], which covers most of the region of the HRB.

In recent years, with the launch of the GRACE satellite, the water depletion in the HRB has been a research hotspot, given its significant depletion rate and important role in China [43,44,45,46,47,48,49]. Most of the water depletion is caused by agricultural irrigation in the HRB [17,44,50]. Previous study has tried to quantify the climate-driven and human-induced TWSA in the HRB. For example, Yi et al. [39] found a human-induced TWSA trend of −8.4 ± 0.9 km³/yr in the HRB from 2003 to 2014 based on a linear regression of annual precipitation anomalies against GRACE-derived TWSA. In this study, we aim to isolate climate-driven and human-induced contributions to GRACE-derived TWSA in the HRB on the basis of the statistical reconstruction method.

2.2. GRACE Data

The GRACE mascon product is provided by the Center for Space Research (CSR); the University of Texas at Austin is used (http://www2.csr.utexas.edu/grace/). CSR mascons (CSR-M) are provided with a 0.5-degree spatial resolution. The C20 replacement, degree 1 corrections and the glacial isostatic adjustment (GIA) correction have been applied to CSR-M, similar as standard processing for spherical harmonic products. More details about the mascon product can be found in Save et al. [51]. This dataset has been used in several studies to derive global and regional TWSA estimates [10,20,33,49].

We also use two other mascon solutions for intercomparison of reconstruction and trend estimation: one is the Jet Propulsion Laboratory mascons (JPL-M) [52], and the other is the Goddard Space Flight Center mascons (GSFC-M) [53]. The JPL-M data is also processed based on the Level-1 GRACE observations, and the same corrections are applied as CSR-M. The JPL-M uses a priori constraints to estimate global monthly gravity fields in terms of 4551 equal-area, surface spherical-cap mass-concentration functions to minimize the effect of measurement errors [52]. The gain factors are not applied in our study, since they are not suitable to quantify trends from the JPL-M [52]. GSFC-M uses different least squares constraint strategies to develop its mascon solutions relative to JPL-M [53].

The GSFC solution uses many more and much smaller mascons than JPL-M to represent global mass changes. We extract the mascon grids in our study area. The GRACE gridded data used in this study covered a total period of 124 months (January 2003 to December 2013).

2.3. Precipitation Data

Gridded daily precipitation data over the study region are obtained from the China Meteorological Administration (CMA) for the period of 1961–2013. The reconstruction period is from 1967 for discarding six spin-up years as suggested by Humphrey et al. [54]. The gridded data were obtained by partial thin plate smoothing splines interpolation of precipitation observations from 2472 weather stations with daily temporal resolution and 0.5° × 0.5° spatial resolution over the whole of China, and 258 weather stations were used in the HRB (Figure A1). These datasets were validated with gauge observations by generalized cross-validation, Root Mean Square Error (RMSE), absolute error and relative error [55]. Humphrey and Gudmundsson [56] highlighted that the quality of the reconstruction is heavily dependent upon the quality of the input precipitation forcing. Two global precipitation forcing datasets applied by Humphrey and Gudmundsson [56] are also used to compare with the CMA precipitation dataset, including the European Centre for Medium-Range Weather Forecast (ECMWF) re-analysis (ERA-Interim) dataset and the multi-source weighted-ensemble precipitation dataset (MSWEP). We validate these three precipitation datasets with annual mean precipitation statistics from the Haihe River Water Resources Bulletin (HRWRB) (http://www.hwcc.gov.cn/hwcc/wwgj/xxgb/szygb/). The statistics from HRWRB are more reliable due to more precipitation data from local meteorological stations included, but only annual mean precipitation statistics over the whole HRB are available. The comparison with statistics from HRWRB indicates that the CMA precipitation data used in this study are more reliable than those from ECMWF and MSWEP over the HRB (Text A1, Figure A2,

Table A1. Comparison of different precipitation-forcing datasets with precipitation statistics from the Haihe River Water Resources Bulletin on the annual scale.

Method	R	MBE (mm)	RMSE (mm)
CMA	0.99	6.67	9.24
ERA-Interim	0.37	21.85	76.86
MSWEP	0.97	−69.63	70.92

), consistent with previous findings from Pan et al. [57].

2.4. Land Surface Models

The Global Land Data Assimilation System (GLDAS) Noah is also used to estimate climate-related TWSA. In GLDAS LSMs, the TWS is the sum of soil moisture content and snow water equivalent. 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 optimal fields of land surface states and fluxes [58]. The data has a 1-degree spatial resolution over global land, with a monthly temporal resolution.

Climate Prediction Center (CPC) global monthly soil moisture data is also used to estimate climate-related TWSA [59]. The dataset has a 0.5-degree spatial resolution from 1948 to the present. The CPC dataset contains monthly averaged soil moisture water height equivalents.

2.5. Climate-Driven Water Storage Variability

Humphrey et al. [54] proposed a data-driven statistical model where precipitation is filtered using an exponential decay function that is optimized to yield the best correlation with the GRACE-observed TWSA [60]. Here we conduct our study based on the reconstructed climate-driven TWSA. However, we did not consider interannual temperature anomalies [33], since there are non-significant correlation coefficients between interannual temperature and TWSA in the HRB [60]. Different methods with and without considering temperature anomalies [33,54], as well as a new method proposed by Humphrey and Gudmundsson [56], were tested, but results indicate that differences among these methods are negligible in the HRB (Text A2, Figure A3,

Table A2. Comparison of different methods with or without considering the temperature, and a new method from Humphrey and Gudmundsson [56] for the GRACE period. All results are compared with deseasoned and detrended GRACE observations from 2003 to 2013.

Method	NSE	RMSE (mm)	R
This study	0.41	16.34	0.79
Humphrey et al. [54]	0.40	16.38	0.79
Humphrey and Gudmundsson [56]	0.45	16.09	0.80

,

Table A3. Intercomparison of among reconstruction terrestrial water storage anomalies results based on the three methods from 1985 to 2013.

Method	NSE	RMSE (mm)	R
This study vs Humphrey et al. [54]	1.00	1.24	1.00
This study vs Humphrey and Gudmundsson [56]	0.96	6.50	0.98
Humphrey et al. [54] vs Humphrey and Gudmundsson [56]	0.95	7.03	0.98

). The climate-driven TWSA is formulated as follows:

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

(1)

where TWSA_rec is the reconstructed, deseasoned and detrended TWS variations; $P_{inter+subseas}^{\tau }$ corresponds to the deseasoned and detrended monthly precipitation. The parameters τ and β are calibrated against the monthly TWSA derived from deseasoned and detrended GRACE-observed TWSA. We apply this method on the basin scale, which is different from that of the grid scale considered in Humphrey et al. [54].

2.6. Quantifying the Contributions of Climate-Driven and Human-Induced TWSA Trends

As shown in Figure 1c, we reconstruct the long-term, climate-driven TWSA from January 1967 to December 2013. The linear trend of TWSA for the reference period or the reconstruction period (i.e., 1967.01–2013.12) is zero, according to the definition of reconstruction [54], but the TWSA time series for the study period (i.e., 2003.01–2013.12) have a trend of 4.56 ± 0.63 mm/yr (expressed as: ${\mathrm{Trend}}_{196701–201312}^{2003–2013}$). It should be noted that the trend of reconstructed climate-driven TWSA for the study period will vary corresponding to different reference periods (see

Table 1. Climate-driven TWSA trends for the period of 2003–2013 corresponding to different reference periods (unit: mm/yr).

Number i	Reference Period	Trend
1	1967.01–2013.12	4.56 ± 0.63 (${\mathrm{Trend}}_{196701–201312}^{2003–2013}$)
2	1967.02–2013.12	4.56 ± 0.63 (${\mathrm{Trend}}_{196702–201312}^{2003–2013}$)
…	…	…
157	1980.01–2013.12	3.82 ± 0.63 (${\mathrm{Trend}}_{198001–201312}^{2003–2013}$)
…	…	…
324	1993.12–2013.12	3.97 ± 0.63 (${\mathrm{Trend}}_{199312–201312}^{2003–2013}$)
Mean		4.31 ± 0.72 (${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$)

and Movie S1). We further estimate an upper and lower bound of the climate-driven trend by considering as many possible reference periods. The mean climate-driven trend is estimated as Equation (2).

$${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\mathrm{Trend}}_{ith–201312}^{2003–2013}}$$

(2)

where ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$ is the mean climate-driven trend, ith is the start year/month of the reference period (i.e., 1967.01, 1967.02, 1967.03, ⋯, 1993.12), and n represents the number of considered reference periods, equivalent to the number of months from January 1967 to a suitable time (i.e., December 1993 and n = 324, see Table 1). We select December 1993 to December 2013 as the last reference period, so that the reconstruction period can be at least 20 years or longer. The superscript of ${\mathrm{Trend}}_{ith–201312}^{2003–2013}$, i.e., 2003–2013, represents the study period to estimate the climate-driven trend. The subscript of ${\mathrm{Trend}}_{ith–201312}^{2003–2013}$ represents the reference period used for reconstruction. The climate-driven trend estimated under a longer reconstruction period tends to be more stable, as it is less influenced by extreme precipitation anomalies in some years.

Then, the human-induced TWSA trend is estimated as follows:

$${\mathrm{Trend}}_{\mathrm{HI}}{=\mathrm{Trend}}_{\mathrm{GRACE}}{-\mathrm{Trend}}_{\mathrm{CD}}$$

(3)

where Trend_HI is the trend of human-induced TWSA, Trend_GRACE is the total trend of GRACE-derived TWSA, and Trend_CD is the climate-driven trend estimate.

2.7. Error Estimation

The uncertainty of TWSA is estimated following the approach used in Landerer and Swenson [8] and Scanlon et al. [10]. Details can be found in the supporting information in Scanlon et al. [10]. The uncertainty of ${\mathrm{Trend}}_{ith–201312}^{2003–2013}$ is estimated as the formal error of trend estimate by applying least squares to fit a linear trend to the constructed climate-driven TWSA time series from 2003 to 2013 using the reference period of ith-201312 [33,61]. The error of ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$ is estimated as the combination of the mean formal error of different ${\mathrm{Trend}}_{ith–201312}^{2003–2013}$ and the standard deviation (STD) of different trend estimates considering all the different reference periods. The equation is as follows:

$${\sigma }_{\mathrm{TCD}}=\sqrt{{\sigma }_T^2+\mathrm{STD}{(Tts)}^2}$$

(4)

where ${\sigma }_{\mathrm{TCD}}$ is the error of ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$, ${\sigma }_T$ is the mean formal error of different ${\mathrm{Trend}}_{ith–201312}^{2003–2013}$, and $\mathrm{STD}(Tts)$ represents the STD of different trend estimates.

The error of the ${\mathrm{Trend}}_{\mathrm{HI}}^{2003–2013}$ is estimated using conventional error propagation as Equation (5), with ${\sigma }_{\mathrm{GRACE}}$ representing the error of the GRACE TWSA trend.

$${\sigma }_{\mathrm{THI}}=\sqrt{{\sigma }_{\mathrm{TCD}}^2+{\sigma }_{\mathrm{GRACE}}^2}$$

(5)

The uncertainties of annual climate-driven and human-induced TWSA are also estimated. We first estimate the uncertainties of the monthly TWSA time series, then the uncertainty of reconstructed TWSA, which is estimated by the root mean square error (RMSE_rec-obs) between reconstructed and GRACE observations [54], and then variations of monthly TWSA (var) for every month which are estimated from the STDs of 324 climate-driven TWSA time series (see Section 2.6). Thus, the uncertainties of annual climate-driven TWSA can be estimated as follows:

$${\sigma }_{\mathrm{AnCD}}=\sqrt{\frac{{\mathrm{RMSE}}_{rec–obs}^2+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$12$}\right.{\displaystyle \sum _{i=1}^{12}{\mathrm{var}}^2}}{12}}$$

(6)

where ${\sigma }_{\mathrm{AnCD}}$ is the uncertainty of annual climate-driven TWSA.

The uncertainty of annual human-induced TWSA is estimated with conventional error propagation considering both annual uncertainties of GRACE and climate-driven TWSA, similar to Equation (5). It should be noted that the uncertainties provided in this study are provided as 68.3% confidence intervals.

3. Results

3.1. Trend of Climate-driven and Human-Induced TWSA

As shown in Figure 1b, our reconstructed time series of TWSA agree well with the GRACE results, with a high correlation coefficient (R) of 0.79 and Nash-Sutcliffe Efficiency (NSE) of 0.41. The RMSE of two time series of TWSA is 1.6 cm (Table A2), in the same level of GRACE errors, which also indicates the reliability of our reconstructed results.

As shown in Figure 1c, although the trend of reconstructed time series of TWSA from 1967 to 2013 is zero according to the definition of reconstruction method, there are significant fluctuations on interannual to decadal time scales. We compute the trend of reconstructed climate-driven TWSA for the period of 2003–2013 using different reference periods as explained by Equation (2) in Section 2.6 (i.e., ${\mathrm{Trend}}_{196701–201312}^{2003–2013}$, ${\mathrm{Trend}}_{196702–201312}^{2003–2013}$, ⋯,${\mathrm{Trend}}_{199312–201312}^{2003–2013}$). Further, the time series of human-induced trends of TWSA for the period of 2003–2013 are estimated by removing climate-driven TWSA trends from GRACE-derived TWSA trends following Equation (3).

The GRACE-derived TWSA trend from January 2003 to December 2013 is −8.56 ± 0.79 mm/yr (i.e., −2.74 ± 0.25 km³/yr), while the Trend_CD ranges from 3.63 to 4.77 mm/yr depending on the various reference periods (Figure 2, Text A3). The mean Trend_CD is estimated to be 4.31 ± 0.72 mm/yr (i.e., 1.38 ± 0.23 km³/yr). The Trend_CD time series show some fluctuations over the different reference periods, but their standard deviation is only 0.34 mm/yr. As shown in Figure 2a, when the reconstruction period becomes longer than 37 years (i.e., the starting year ranges from 1967 to 1977), the trend estimates tend to become more stable.

As a direct result of Equation (3), variations of Trend_HI are symmetric to the Trend_CD time series (Figure 2a). The Trend_HI is estimated to be −12.87 ± 1.07 mm/yr (i.e., −4.12 ± 0.34 km³/yr, Figure 2b), which is stronger than the GRACE-derived TWSA trend.

3.2. Interannual Variations of TWSA

To better illustrate the year-to-year variability of TWSA, monthly GRACE-derived, climate-driven and human-induced TWSA are averaged to annual means (Figure 3a). The climate-driven TWSA increased with fluctuations, while the human-induced TWSA continuously decreased, indicating a continuous human-induced water mass loss. The climate-driven TWSA decreased from 2004 to 2007 and increased continuously after 2007, which is consistent with year-to-year precipitation anomalies during the same time period (Figure A4). However, the human-induced TWSA shows a long-term continuous decrease without significant interannual variability, indicating a stable negative effect on water storage from anthropogenic activities (Figure 3a). As the leading factor of human activities, the long-term extensive water consumption with a rate of ~26 km³/yr based on the statistics from HRWRB causes the continuous decrease of human-induced TWSA (Figure A4).

Further, the year-to-year terrestrial water fluxes (TWF) are calculated as the difference between annual mean TWSA in two consecutive years. As shown in Figure 3b, the interannual GRACE TWF are dominated by the climate-driven TWF. The climate-driven TWF between 2003 and 2004 is 16.80 ± 2.15 km³/yr, corresponding to the positive precipitation anomaly in 2003 and 2004. The abnormal negative precipitation anomaly in 2006 (i.e., −79.14 mm relative to average precipitation of 1961 to 2013, Figure A4) corresponds to minimum TWF (−6.92 ± 2.15 km³/yr) between 2005 and 2006. The human-induced TWF are all in negative values, ranging from −6.28 ± 3.06 to −1.09 ± 3.06 km³/yr.

4. Discussion

4.1. Comparison of Different GRACE Solutions and Methods

The method to separate contributions of climate-driven and human-induced activities in this study is also applied by using JPL-M and GSFC-M. As shown in

Table 2. Estimates of GRACE-derived, climate-driven, and human-induced TWSA trends in the HRB from 2003 to 2013 based on different GRACE solutions and methods (unit: mm/yr).

GRACE-Derived TWSA Trend	Climate-Driven TWSA Trend	Human-Induced TWSA Trend	Description (GRACE Solution, Method)
−8.56 ± 0.79	4.31 ± 0.72	−12.87 ± 1.07	CSR-M, reconstruction
−9.86 ± 0.69	3.35 ± 0.66	−13.21 ± 0.96	GSFC-M, reconstruction
−14.20 ± 0.99	5.66 ± 0.95	−19.86 ± 1.37	JPL-M, reconstruction
−8.56 ± 0.79	1.09 ± 0.70	−9.65 ± 1.06	CSR-M, climate-driven TWSA from Noah
−8.56 ± 0.79	2.48 ± 0.69	−11.04 ± 1.05	CSR-M, climate-driven TWSA from CPC
−8.56 ± 0.79	2.49 ± 2.58	−11.05 ± 2.46	CSR-M, linear regression (Yi et al. 2016)
−8.56 ± 0.79	4.33 ± 0.63	−12.89 ± 1.01	CSR-M, reconstruction and linear regression
−8.56 ± 0.79	4.31 ± 0.71	−12.87 ± 1.07	CSR-M, reconstruction (detrended annual precipitation used)

, the GRACE-derived TWSA trend from JPL-M is significantly larger than those from CSR-M and GSFC-M. Scanlon et al. [10] also found larger trends in JPL-M after a global evaluation of different GRACE solutions, especially in medium and small basins (≤500,000 km²), which can be caused by the different processing strategy and a priori constraints used in JPL-M [52,62]. The climate-driven trend estimates from CSR-M, GSFC-M and JPL-M are 4.31 ± 0.72, 3.35 ± 0.66, and 5.66 ± 0.95 mm/yr, respectively. The mean value of the three estimates is 4.44 mm/yr, and it is closer to the CSR-M estimate, which we used finally in this study. Due to the positive trends of climate-driven TWSA during 2003–2013, the negative trends of the human-induced TWSA are all larger than those of GRACE-derived total TWSA.

We reproduce the method of Yi et al. [39] by assuming a simple linear relationship between annual precipitation anomalies and annual TWF (Figure A5). The Trend_HI estimated by this method is −11.05 ± 2.46 mm/yr, agreeing reasonably well with our Trend_HI estimate within uncertainties (Table 2). Note that only precipitation anomalies and TWF on an annual scale can be used to establish a reasonable linear relationship between them in this method. However, in our data-driven reconstruction method, a more complicated decayed response of TWSA to daily precipitation is considered [54].

We also follow the method of Huang et al. [14], who take the LSM simulations as the climate-related TWSA. The two climate-driven TWSA time series are relatively similar (Figure A6) in the HRB. The corresponding Trend_CD from the Noah and CPC models are 1.09 ± 0.70 and 2.48 ± 0.69 mm/yr, respectively, during 2003–2013 (Table 2). The increased TWSA from LSMs are consistent with the positive climate-driven TWSA trend reconstructed in this study. However, climate-driven TWSA trends estimated by LSMs appear smaller than our estimates. This finding is consistent with a recent study indicating that the amplitude of TWSA trends is systematically underestimated in LSMs (Scanlon et al. 2018), partly because of missing storage compartments such as lakes and aquifers which, on the other hand, are implicitly included in our data-driven reconstruction method. Further, human-induced TWSA time series are estimated by subtracting model-simulated, climate-driven TWSA from GRACE-derived TWSA. The corresponding Trend_HI estimates based on the Noah and CPC models are −9.65 ± 1.06 and −11.04 ± 1.05 mm/yr, respectively (Table 2), which are smaller compared to our estimates using the data-driven reconstruction method due to the underestimation of climate-driven TWSA in LSMs.

4.2. Impact of Annual Precipitation Trend

The intrinsic trend in annual precipitation could also contribute to the changes of TWS, as over the Tibetan Plateau in China [17,63], northern North America, Tropical western Africa and Northern Congo [17]. In this study, we assumed that the long-term climatological precipitation remains stable over the reference period, so the trends of the time series of precipitation and climate-driven TWSA are zero over the reference period. In fact, long-term annual precipitation trend may cause an increase or decrease of TWS. In other words, the long-term climate-driven TWSA trend might not be zero.

We examined the trend of annual precipitation for different reference periods (i.e., 1967–2013, 1968–2013, ⋯, 1993–2013) using the Mann-Kendall (M-K) trend test (Figure 4a). There are two significant annual precipitation trends during 1980–2013 and 1992–2013 corresponding to the two low values of ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$ estimates (Figure 4a), which indicates the impact of annual precipitation trend on the ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$ estimation. The precipitation trend remains stable for the reference periods from 1967–2013 to 1977–2013 (Figure 4a), so the ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$ estimates are also stable during these reference periods (Figure 2a). A simple linear regression was further performed between annual precipitation trends (predictor) and ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$ estimates (dependent variable) during different reference periods (Figure 4b). The regression coefficient shows that the ${\mathrm{Trend}}_{\mathrm{CD}}^{2003–2013}$ is 4.33 ± 0.63 mm/yr when the annual precipitation trend is zero (Table 2), which is consistent with our estimate. Here we assume that the long-term, climate-driven TWSA trend is zero when the annual precipitation trend is zero.

In addition, we also check the Trend_CD and Trend_HI estimates by removing the annual precipitation trend before reconstruction. The trend of the annual precipitation during 1961–2013 is found to be −0.87 mm/yr² based on the M-K trend test, but not significant (p > 0.1). Nevertheless, we still consider its contribution as follows. We first removed the annual precipitation trend by adding the corresponding amount of precipitation in the daily precipitation time series, so the annual precipitation trend for the period of 1961–2013 is close to zero after this processing. The detrended daily precipitation time series are further used to reconstruct the climate-driven TWSA. Finally, the Trend_CD and Trend_HI are estimated to be 4.31 ± 0.71 and −12.87 ± 1.07 mm/yr, respectively (Table 2). There is almost no difference between these estimates and the foregoing trend estimates using the original precipitation. Thus, we can conclude that the trend of annual precipitation has no obvious influence on our Trend_CD and Trend_HI estimates in the HRB.

However, it is worth noting that special attention might be required when applying the method to a region with a significant annual precipitation trend during the reconstruction period. Technically, the potential impact of a precipitation trend can be tested by applying regression analysis between annual precipitation trends and Trend_CD estimates proposed in this study.

5. Conclusions

In this study, we use a data-driven reconstruction method to quantify the contributions of climate and anthropogenic activities to water storage changes. This method takes advantage of the information from GRACE and in situ precipitation data. The HRB, suffering a long-term water depletion, is taken as a case study. Although the apparent rate of TWSA in the HRB measured by GRACE is −8.56 ± 0.79 mm/yr over 2003–2013, we estimate a climate-driven TWSA trend of 4.31 ± 0.72 mm/yr (i.e., 1.38 ± 0.23 km³/yr) due to the wet period in the past decade relative to the climate of the past half century. Therefore, the human-induced TWSA trend is −12.87 ± 1.07 mm/yr (i.e., -4.12 ± 0.34 km³/yr) for the period of 2003–2013, which indicates that the human-induced water depletion in the HRB is more serious (i.e., ~50% larger) than that estimated by GRACE observations alone.

While annual climate-driven TWSA shows some fluctuations with a positive trend, prolonged anthropogenic water abstractions lead to a continuous decrease in TWSA. Although the climate-driven TWSA increase partially compensates the human-induced water depletion in the case of the HRB, it should be noted that there might also be opposite cases in which a negative climate-driven TWSA trend along with an existing anthropogenic TWSA depletion will exacerbate the total water loss in some regions of the world.

The potential effect of the long-term annual precipitation trend in the HRB is also investigated, with the results showing that it has almost no effect in the HRB. However, we recommend that caution should be exercised in regions with the significant precipitation trend when applying the reconstruction method.

The method presented in this study provides a way to isolate the individual contribution of climate variability or human activities to decadal water storage changes observed by GRACE, which would definitely benefit the interpretation of observed net changes in water storage and water resources management. Our results indicate that although the climate-driven water storage varies corresponding to the precipitation changes, the water storage variations due to anthropogenic activities show a continuous decreasing rate in the HRB, which highlights the possible underestimation or overestimation of actual human-induced water depletion in other regions of the world when GRACE data alone was used to quantify regional water storage changes.