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

: Terrestrial water storage (TWS) can be inﬂuenced 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 ﬁnd 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.


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 3 /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.

Study Area
The Haihe River Basin (HRB) is selected as a case study, which is located in North China (~320,000 km 2 ; 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. 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 3 /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.

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 sphericalcap 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- 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 3 /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 3 /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.

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).

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), consistent with previous findings from Pan et al. [57].

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.

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 Remote Sens. 2019, 11, 3050 5 of 19 differences among these methods are negligible in the HRB (Text A2, Figure A3, Table A2, Table A3). The climate-driven TWSA is formulated as follows: where TWSA rec is the reconstructed, deseasoned and detrended TWS variations; P τ inter+subseas 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].

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: Trend 2003-2013 196701-201312 ). 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 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).
Trend 2003-2013 ith-201312 (2) where Trend 2003-2013 CD 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 Trend [2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013] ith-201312 , i.e., 2003-2013, represents the study period to estimate the climate-driven trend. The subscript of Trend [2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013] ith-201312 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: 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.

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 Trend [2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013] ith-201312 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]. HI is estimated using conventional error propagation as Equation (5), with σ GRACE representing the error of the GRACE TWSA trend.
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: var 2 12 (6) where σ 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.

Trend of Climate-driven and Human-Induced TWSA
As shown in Figure 1b 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 3 /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 3 /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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 19 reference periods (Figure 2, Text A3). The mean TrendCD is estimated to be 4.31 ± 0.72 mm/yr (i.e., 1.38 ± 0.23 km 3 /yr). The TrendCD 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.

Interannual Variations of TWSA
To better illustrate the year-to-year variability of TWSA, monthly GRACE-derived, climatedriven 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 3 /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  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 3 /yr, Figure 2b), which is stronger than the GRACE-derived TWSA trend.

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 3 /yr based on the statistics from HRWRB causes the continuous decrease of human-induced TWSA ( Figure A4). Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 19

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, 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 2 ), 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 TrendHI estimated by this method is −11.05 ± 2.46 mm/yr, agreeing reasonably well with our TrendHI 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 climaterelated TWSA. The two climate-driven TWSA time series are relatively similar ( Figure A6) in the HRB. The corresponding TrendCD 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 The human-induced TWF are all in negative values, ranging from −6.28 ± 3.06 to −1.09 ± 3.06 km 3 /yr.

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, 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 2 ), 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.

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  (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 TrendCD and TrendHI 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 2 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 TrendCD and TrendHI 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 TrendCD and TrendHI 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 TrendCD estimates proposed in this study.

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 3 /yr) due to the wet period in the past decade relative to the climate of 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 2 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.

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 3 /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 3 /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. Error (RMSE) (Table A1). Results indicate that precipitation data from CMA are more reliable in the HRB, which is in line with Pan's conclusion [57].

Text A2. Comparison of Different Reconstruction Methods
Results with or without considering the interannual temperature in the reconstruction are shown in Figure A3, which are applied in Humphrey et al. [54] and Zhong et al. [33], respectively. The interannual temperature data are obtained from the ERA-Interim atmospheric reanalysis [64]. A new reconstruction method proposed by Humphrey and Gudmundsson [56] is also taken into comparison. The results are evaluated with three indices: the Nash-Sutcliffe efficiency (NSE) [66], RMSE and R (Table A2, Table A3). The comparison results indicate that the impact of temperature is limited in the HRB. Considering that there is non-significant R between the interannual temperature and TWSA, we choose the method without considering the temperature.

Text A3. A Relatively Conservative Estimate of Trend CD and Trend HI
If we assume the climate-driven trend is zero for the longest reference period (i.e., Trend [2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013] 196701-201312 = 0), then Trend CD is estimated to be 4.56 mm/yr. The subsequent trend time series are used to estimate the upper and lower bounds of Trend CD (Figure 2a) new reconstruction method proposed by Humphrey and Gudmundsson [56] is also taken into comparison. The results are evaluated with three indices: the Nash-Sutcliffe efficiency (NSE) [66], RMSE and R (Table A2, A3). The comparison results indicate that the impact of temperature is limited in the HRB. Considering that there is non-significant R between the interannual temperature and TWSA, we choose the method without considering the temperature.

Text A3. A Relatively Conservative Estimate of TrendCD and TrendHI
If we assume the climate-driven trend is zero for the longest reference period (i.e.,