Terrestrial Water Storage in China: Spatiotemporal Pattern and Driving Factors

China is the largest agricultural country with the largest population and booming socio-economy, and hence, remarkably increasing water demand. In this sense, it is practically critical to obtain knowledge about spatiotemporal variations of the territorial water storage (TWS) and relevant driving factors. In this study, we attempted to investigate TWS changes in both space and time using the monthly GRACE (Gravity Recovery and Climate Experiment) data during 2003–2015. Impacts of four climate indices on TWS were explored, and these four climate indices are, respectively, El Niño Southern Oscillation (ENSO), Indian Ocean Dipole (IOD), North Atlantic Oscillation (NAO), and Pacific decadal oscillation (PDO). In addition, we also considered the impacts of precipitation changes on TWS. We found significant correlations between climatic variations and TWS changes across China. Meanwhile, the impacts of climate indices on TWS changes were shifting from one region to another across China with different time lags ranging from 0 to 12 months. ENSO, IOD and PDO exerted significant impacts on TWS over 80% of the regions across China, while NAO affected TWS changes over around 40% of the regions across China. Moreover, we also detected significant relations between TWS and precipitation changes within 9 out of the 10 largest river basins across China. These results highlight the management of TWS across China in a changing environment and also provide a theoretical ground for TWS management in other regions of the globe.


Introduction
The warming climate, and its impacts on the hydrological cycle at regional and global scales, has aroused remarkable human concerns in recent decades [1][2][3]. Intensifying hydrologic cycle, due to the warming climate has the potential to trigger occurrences of hydrometeorological hazards, such as floods and droughts [4][5][6]. A range of researches indicated intensified drought severity and expanded drought-affected regions, due to warming climate [7][8][9] although what are the drought behaviors in warming climate is still a scientific issue open for debate [10,11]. In addition, Hu et al. [12] found that, at a global scale, the occurrence of floods is increasing; therefore, flood-affected population, flood-related mortalities are increasing as well. Meanwhile, amplification of floods is featured by  Table 1.

GRACE Data
In this study, TWS data from release-5 GRACE were retrieved from NASA's GRACE Tellus website [42]. The data have been pre-processed to remove the noise signal from the atmosphere and ocean [43,44]. The spatial resolutions of the monthly TWS data are 1° by 1°. There are currently three publicly released GRACE TWS data by CSR (Center for Space Research), GFZ (Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum) and JPL (Jet Propulsion Laboratory), respectively. To diminish the uncertainties associated with data processing, the overall average TWS data were calculated based on the three research centers.
In regional studies, the application of appropriate scaling factor methods is essential for accurate quantification of the TWS monitored by GRACE satellites [45], which can correct and restore the GRACE signal loss in low-pass filtering processes, such as strip removal, truncation and filtering. Since Community Land Model 4.0 (CLM4.0) from National Center for Atmospheric Research (NCAR) explains the interaction between surface and groundwater, as well as irrigation and river diversion, and provides reliable flow and hydrological status information in areas where human activity is dense [45], we utilized the scaling factor from NCAR's CLM4.0 to correct and restore the GRACE signal loss during low-pass filtering [43]. Furthermore, the terrestrial water storage anomaly (TWSA)  Table 1.

GRACE Data
In this study, TWS data from release-5 GRACE were retrieved from NASA's GRACE Tellus website [42]. The data have been pre-processed to remove the noise signal from the atmosphere and ocean [43,44]. The spatial resolutions of the monthly TWS data are 1 • by 1 • . There are currently three publicly released GRACE TWS data by CSR (Center for Space Research), GFZ (Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum) and JPL (Jet Propulsion Laboratory), respectively. To diminish the uncertainties associated with data processing, the overall average TWS data were calculated based on the three research centers.
In regional studies, the application of appropriate scaling factor methods is essential for accurate quantification of the TWS monitored by GRACE satellites [45], which can correct and restore the GRACE signal loss in low-pass filtering processes, such as strip removal, truncation and filtering. Since Community Land Model 4.0 (CLM4.0) from National Center for Atmospheric Research (NCAR) explains the interaction between surface and groundwater, as well as irrigation and river diversion, and provides reliable flow and hydrological status information in areas where human activity is dense [45], we utilized the scaling factor from NCAR's CLM4.0 to correct and restore the GRACE signal loss during low-pass filtering [43]. Furthermore, the terrestrial water storage anomaly (TWSA) data were obtained by subtracting the monthly TWS data by a historical mean from 2004 to 2009 according to the anomaly baseline reported in GRACE product highlights (http://www2.csr.utexas.edu/grace/RL05_ mascons.html). The uncertainties of TWSA time series were evaluated at 95% confidence level.

Precipitation Data
TWS changes were believed to be associated with precipitation variations [46]. We adopted three monthly precipitation datasets, i.e., Global Precipitation Climatology Project (GPCP), Global Precipitation Climatology Centre (GPCC), and NOAA's precipitation reconstruction over Land (NOAA) in this study (http://www.esrl.noaa.gov/psd/data/gridded/tables/precipitation.html). The spatial resolutions are 2.5 • by 2.5 • for GPCP and 1 • by 1 • for GPCC and NOAA, and GPCP dataset was interpolated to a spatial resolution of 1 • by 1 • to match the other datasets [47].
The three datasets were averaged to quantify relations between precipitation changes and TWS variations in both space and time. Among the three datasets, they are in similar dynamic patterns, although there would be systemic variations. Hence, the uncertainty can be estimated in two steps [47]: (1) Remove the mean value in each dataset, the systemic variations can also be removed in this step for the reason that the mean value in each dataset would be zero; (2) for each grid point, it has three data value from the three datasets, then the variation among the three values is set as the uncertainty. By this way, the systemic variation would not be included in the uncertainties.

Climate Index Data
To investigate impacts of climatic systems on TWS changes in China, four climatic indices representing four climatic drivers were used in this study, i.e., the Multivariate ENSO Index (MEI), Indian Ocean Dipole Mode Index (DMI), North Atlantic Oscillation Index (NAO) and Pacific Decadal Oscillation Index (PDO). MEI is a method for describing the intensity of ENSO events and was considered to be the most comprehensive index for monitoring ENSO events [48], since it combines the analysis of multiple meteorological and oceanic components. The MEI datasets used in this study were extracted from NOAA's Physical Sciences Division (http://www.esrl.noaa.gov/psd/enso/mei/). DMI can measure the zonal sea surface temperature gradient in the equatorial Indian Ocean [49,50] and indicate the dynamic of IOD. IOD is an aspect of the overall global climate cycle, interacting with similar phenomena, such as ENSO in the Pacific Ocean [35]. The DMI datasets were obtained from the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) and can be downloaded from https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/DMI/.
NAO is a weather phenomenon in the North Atlantic that controls the strength and direction of westerlies and the storm tracks across North Atlantic by the fluctuations in the pressure difference at sea level between the Icelandic low and the Azores high [51]. The NAO index is derived from NOAA's Climate Prediction Center (http://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/nao.shtml).
PDO is an active and recurring oceanic, atmospheric climate change pattern across the mid-latitude Pacific basin. The PDO index is the main empirical orthogonal function of the monthly sea surface temperature anomaly (SST-A) after the surface sea surface temperature minus the global mean sea level temperature in the North Pacific Ocean (north of 20 • N), and is the normalized principal component time sequence [52]. The PDO index adopted in this study was taken from the Earth System Research Laboratory (ESRL) of NASA (http://www.esrl.noaa.gov/psd/data/correlation/pdo.data). The period of the values for MEI, DMI, NAO and PDO covers from 2002 to 2015 here in this study.

Seasonal Decomposition of Time Series by LOESS
To better assess the annual variations of TWSA, a standard time series decomposition method based on Locally weighted regression smoother (LOESS), namely, Seasonal-Trend decomposition by Loess model (STL), was used in this paper. Loess Regression is a non-parametric method used to smoothen a volatile time series, and least squares regression is performed in localized subsets, which makes it a suitable choice for smoothing any numerical vector [53].
We take x i and y i , for i = 1 to n, as the samples of independent and dependent variables, respectively. The loess regression curve,ĝ(x), is a smoothing of y given that x can be computed for any value along the scale of the independent variable. It means that loess can deal with missing values and detrend the seasonal component in a straightforward way in STL. For STL, only the case of one independent variable is needed, although loess can smooth y as a function of any number of independent variables.
g(x) is estimated as follows: A positive integer, q, is chosen. The q values of the x i that are closest to x are selected for the moment that q ≤ n and each is given a neighborhood weight based on its distance from x. We assume that λ q (x) is the distance of the qth farthest x i from x, and W is the tricube weight function: The neighborhood weight for any Thus, the x i close to x have the largest weights; the weights decrease as the x i increase in distance from x and become zero at the qth farthest point. The next step is to fit a polynomial of degree d to the data with weight v i (x) at (x i , y i ). The value of the locally-fitted polynomial at x isĝ(x). In this paper, we will use d = 1 or 2; that is, the fitting is locally-linear or locally-quadratic. Now suppose that q > n. λ n (x) is the distance from x to the farthest x i . For q > n we define λ q (x) by: Then we proceed as before in the definition of the neighborhood weights using this value of λ q (x). Suppose each observation (x i , y i ) has a weight p i that expresses the reliability of the observation relative to the others. For example, if the y i had variances σ 2 k i where the k i were known, then p i might be 1/k i . We can incorporate these weights into the loess smoothing in a straightforward way by using p i v i (x) as the weights in the local least-squares fitting. As we will see, this provides a mechanism by which we can easily build robustness into STL [53].
The STL processor enables conveniently decomposing time series into trend, seasonal and noise components, which can be defined in the form of: where O t is the observed data at time t, T t denotes the trend component, S t refers to the seasonal component, and e t means the noise component. The noise component e t clarifies the remaining variation in the time series after trend and seasonal component [53]. Following the above STL procedure, the TWSA time series data during 2003-2015 across China were de-seasonalized and smoothed in this paper. The coding environment and calculating process were done in R [54].

Linear Regression
The total amount of water is equal during the process of the water cycle, which can be expressed as the water budget equation [55]: where dS/dt is total water storage change, P is precipitation, ET is evapotranspiration and R is runoff.
To explore how changes in annual precipitation force the variations in TWS, we assumed their linear relationship [46] as: where m is the annul TWS change between two successive years (in mm), p is the annual precipitation (in mm/year), p 0 is the long-term annual average precipitation (in mm/year), a is the slope of the regression in (year −1 ), i.e., the sensitivity factor of TWS to precipitation, and b is the average annual depletion water storage without the influence of precipitation (in mm). The linear regression with the ordinary least squares method (OLS) was applied to weigh the model parameters a and b.

Spatial and Temporal Variability
For the period from January 2003 to December 2015, mean values of monthly TWSA and de-seasonalized TWSA for China were computed, and the annual average, together with the average seasonal rate of change in TWS images were also processed to show the spatial and temporal variability of the TWS through time.

Cross-correlation Analysis
To concern TWS variations during study periods to the large-scale climate modes, such as ENSO, IOD, NAO and PDO across China, spatial cross-correlation analysis was applied to quantify the relationships between TWS and the major climate drivers. As the influence of climate modes on TWS is expected to show finite time lags (0 to 12 months), we only discussed positive time lags in which climate modes changes are leading the variations of TWS. We investigated the cross-correlations between TWSA and the four climate indices on a per-pixel level at different time lags by shifting the climate indices time series forward one month per time. Following this procedure, the maximum Pearson's correlation coefficient (r) and the corresponding time lag (in months) were then gained to revel the strength of the relationships between TWS and four climate modes with the optimal lag at any given pixel across China. The Pearson's coefficient r was calculated as: where j = 0, 1, . . . , 12, and r j is the Pearson's correlation coefficient in 0 to 12 months lag, respectively, and n is the sample size, i.e., the length of the time series, x refers to the climate indices, y is the TWSA, x and y is the sample average of the climate indices and TWSA, respectively.

Temporal Variation Characteristics of the TWS across China
The GRACE gravity datasets derived monthly TWSA with its smoothed and de-seasonalized values during 2003-2015 were shown in Figure 2. We found from Figure 2 strong TWS fluctuations at both seasonal and inter-annual scales with basically seasonal peak values in summer (June, July and August, JJA), while trough values in winter (December, January and February, DJF). We also found  Figure 3 illustrates the fluctuation ranges in MEI, DMI, NAO and PDO climate indices from 2002 to 2015. Among the four indices, MEI is an index used to characterize the intensity of ENSO event based on six climatic variables, i.e., sea-level pressure, zonal and meridional surface winds, sea surface temperature, surface air temperature and total cloudiness fraction of the sky. Positive MEI index indicates starting of El Niño phenomenon and negative MEI index the La Niña phenomenon. DMI index characterizes the dynamic changes of IOD events. Positive DMI index indicates the positive phase of the IOD event, and the sea surface temperature in the western Indian Ocean is higher than the average level, the precipitation is greater than the average level, the sea temperature in the eastern Indian Ocean is correspondingly decreased, and vice versa. The NAO exhibits considerable inter-seasonal and interannual variability, and prolonged periods (several months) of both positive and negative phases of the pattern are common. Positive NAO index indicates that the pressure ratio between the Iceland low and the Azores high is greater than the average, west wind dominates and blows warm air; while NAO index is negative, the pressure difference between the two places is below average, the winds with cold air from east and northeast are more frequent. For PDO index, when PDO is in a positive phase, the Pacific would become cooler, and part of the eastern ocean become warmer; while, when PDO is in a negative phase, the eastern Pacific becomes cooler and part of the western ocean become warm.  Among the four indices, MEI is an index used to characterize the intensity of ENSO event based on six climatic variables, i.e., sea-level pressure, zonal and meridional surface winds, sea surface temperature, surface air temperature and total cloudiness fraction of the sky. Positive MEI index indicates starting of El Niño phenomenon and negative MEI index the La Niña phenomenon. DMI index characterizes the dynamic changes of IOD events. Positive DMI index indicates the positive phase of the IOD event, and the sea surface temperature in the western Indian Ocean is higher than the average level, the precipitation is greater than the average level, the sea temperature in the eastern Indian Ocean is correspondingly decreased, and vice versa. The NAO exhibits considerable inter-seasonal and interannual variability, and prolonged periods (several months) of both positive and negative phases of the pattern are common. Positive NAO index indicates that the pressure ratio between the Iceland low and the Azores high is greater than the average, west wind dominates and blows warm air; while NAO index is negative, the pressure difference between the two places is below average, the winds with cold air from east and northeast are more frequent. For PDO index, when PDO is in a positive phase, the Pacific would become cooler, and part of the eastern ocean become warmer; while, when PDO is in a negative phase, the eastern Pacific becomes cooler and part of the western ocean become warm.

Cross-Correlations of Climate Modes and TWS
pressure ratio between the Iceland low and the Azores high is greater than the average, west wind dominates and blows warm air; while NAO index is negative, the pressure difference between the two places is below average, the winds with cold air from east and northeast are more frequent. For PDO index, when PDO is in a positive phase, the Pacific would become cooler, and part of the eastern ocean become warmer; while, when PDO is in a negative phase, the eastern Pacific becomes cooler and part of the western ocean become warm. In this paper, cross-correlation analysis based on the Pearson correlation method was used to calculate the correlation between the climatic driver variables and TWS variations.  In this paper, cross-correlation analysis based on the Pearson correlation method was used to calculate the correlation between the climatic driver variables and TWS variations.   Figure 4a shows that MEI was remarkably correlated to TWS over China with great regional variations across China. Around 86.5% areas of China were dominated by significant correlations between TWS and MEI, of which 33.3% areas are characterized by negative correlations 53.2% of the total territory of China by positive correlations (Table 2). Moreover, the regions with positive   DMI was also found to be highly related to TWS, while the time lags are relatively shorter in comparison with the time lags between MEI and TWS. In addition, in most regions, the correlation coefficients between DMI and TWS are exactly opposite to the ones between MEI TWS (Figure 4a and  4b). The results demonstrated that the influence of ENSO events on TWS variations is different from the influence of IOD on TWS. However, unlike the powerful impact of ENSO and IOD on TWS changes in China, NAO was related to TWS in part of the regions of China (Figure 4c). Only 39.7% of the areas in China have significant relationships between NAO and TWS ( Table 2). The primary portions of the regions with positive correlation were mainly located in southern China, while the major negative correlation areas were found mainly in western China (southern part of XbRB). Compared with MEI and DMI, PDO was highly correlated with TWS over China as well, but much more similarity with the correlation between MEI and TWS (Figure 4a, 4b and 4d).

Relations between Precipitation and TWS
To investigate the relations of TWS vs. precipitation, we quantified relations between monthly precipitation and TWS within each river basin across China (Figures 6 and 7). The results indicated that YRB (Figure 6a Figure 4a shows that MEI was remarkably correlated to TWS over China with great regional variations across China. Around 86.5% areas of China were dominated by significant correlations between TWS and MEI, of which 33.3% areas are characterized by negative correlations 53.2% of the total territory of China by positive correlations (Table 2). Moreover, the regions with positive correlation are mainly in northern China with 86.3% areas of ShRB, 54.4% areas of LRB, 62.7% areas of HRB, 49.6% areas of YB and 54.1% areas of XbRB comprised of 71% areas of the total regions with positive correlation. These abbreviated names of the river basins can be found in Table 1. The other regions with positive correlations are mainly located in southwestern China and the lower YRB. However, the regions with active correlations are mainly concentrated within three areas, i.e., southern XbRB, southeastern China and northern HuaiRB. As for the time lag, as shown in Figure 5a, the regions with the highest correlation coefficients and the regions with the longest time lags in months are inconsistent in the spatial sense. More specifically, for the regions in southern XbRB, the correlation coefficients between MEI and TWS range from −0.6 to −0.1, while the time lag in months ranges between 0-12. DMI was also found to be highly related to TWS, while the time lags are relatively shorter in comparison with the time lags between MEI and TWS. In addition, in most regions, the correlation coefficients between DMI and TWS are exactly opposite to the ones between MEI TWS (Figure 4a,b). The results demonstrated that the influence of ENSO events on TWS variations is different from the influence of IOD on TWS. However, unlike the powerful impact of ENSO and IOD on TWS changes in China, NAO was related to TWS in part of the regions of China (Figure 4c). Only 39.7% of the areas in China have significant relationships between NAO and TWS ( Table 2). The primary portions of the regions with positive correlation were mainly located in southern China, while the major negative correlation areas were found mainly in western China (southern part of XbRB). Compared with MEI and DMI, PDO was highly correlated with TWS over China as well, but much more similarity with the correlation between MEI and TWS (Figure 4a,b,d).

Relations between Precipitation and TWS
To investigate the relations of TWS vs. precipitation, we quantified relations between monthly precipitation and TWS within each river basin across China (Figures 6 and 7). The results indicated that YRB (Figure 6a   Specifically, in comparison with other river basins considered in this study, we found larger decreasing magnitude of TWS in the HRB (Figure 6c) with a decreasing rate of −9.36 ± 1.44 mm/a. It should be clarified that the HRB is located in the North China Plain, one of the largest irrigation areas in China. Therefore, massive exploitation of groundwater for irrigation leads to remarkable TWS changes in the HRB [56]. Furthermore, the annual average TWS of HRB (Figure 6c) is also the lowest within all the river basins considered in this study, which is −27.68 mm. For the YB (Figure 6d) and the HuaiRB (Figure 6e), the annual average precipitation of the HuaiRB (71.39 mm) is larger than YRB (37.92 mm), whereas the annual average TWS of the HuaiRB (−23.77 mm) is less than the YRB (−17.25 mm). The decreasing magnitude of TWS in the HuaiRB is −7.08 ± 1.92 mm/a, while the decreasing rate of TWS in the YB is −5.04 ± 1.2 mm/a. Besides, we also found decreasing TWS in the LRB ( Figure  7a) and the XbRB (Figure 7c) and the decreasing rates of TWS in the LRB and XbRB are, respectively, −3.0 ± 1.08 mm/a and −1.32 ± 0.48 mm/a. Typically, it takes time for the process of precipitation turn into runoff and then transforms into TWS. So, we shifted the time series of TWS in 0 to 12 months later than precipitation time series to search the maximum correlation and found different degrees of sensitive responses of TWS to precipitation changes. It can be seen from Figure 8 and Table 3 that most basins' shifted months are less than half years expect HuaiRB (9 months). The most sensitive basin with 0 month shifted is YRB. We also found that precipitation variations are highly related to the annual TWS changes with a correlation coefficient of 0.8 (Figure 8a). In particular, the relationship between precipitation and TWS Specifically, in comparison with other river basins considered in this study, we found larger decreasing magnitude of TWS in the HRB (Figure 6c) with a decreasing rate of −9.36 ± 1.44 mm/a. It should be clarified that the HRB is located in the North China Plain, one of the largest irrigation areas in China. Therefore, massive exploitation of groundwater for irrigation leads to remarkable TWS changes in the HRB [56]. Furthermore, the annual average TWS of HRB (Figure 6c) is also the lowest within all the river basins considered in this study, which is −27.68 mm. For the YB (Figure 6d) and the HuaiRB (Figure 6e), the annual average precipitation of the HuaiRB (71.39 mm) is larger than YRB (37.92 mm), whereas the annual average TWS of the HuaiRB (−23.77 mm) is less than the YRB (−17.25 mm). The decreasing magnitude of TWS in the HuaiRB is −7.08 ± 1.92 mm/a, while the decreasing rate of TWS in the YB is −5.04 ± 1.2 mm/a. Besides, we also found decreasing TWS in the LRB (Figure 7a) and the XbRB (Figure 7c) and the decreasing rates of TWS in the LRB and XbRB are, respectively, −3.0 ± 1.08 mm/a and −1.32 ± 0.48 mm/a. Typically, it takes time for the process of precipitation turn into runoff and then transforms into TWS. So, we shifted the time series of TWS in 0 to 12 months later than precipitation time series to search the maximum correlation and found different degrees of sensitive responses of TWS to precipitation changes. It can be seen from Figure 8 and Table 3 that most basins' shifted months are less than half years expect HuaiRB (9 months). The most sensitive basin with 0 month shifted is YRB. We also found that precipitation variations are highly related to the annual TWS changes with a correlation coefficient of 0.8 (Figure 8a). In particular, the relationship between precipitation and TWS inYRB during 2009-2013 is the highest when compared to that during the entire study period (Figure 9a). The precipitation variation results in annual TWS change of 2.82 ± 1.54 mm/mm (Figure 10a), while other factors result in annual TWS change of 11.28 ± 113.96 mm/a. Among all the basins in China, XnRB has the highest correlation coefficient of 0.86 (Table 3, Figure 9i). It can be seen from Figure 9i that precipitation and TWS changes are consistent with a time lag of three months. In addition, the linear relationship between precipitation and TWS is highly significant with the slope of the linear fitting model of 4.29 ± 1.76 mm/mm (Figure 10i).
inYRB during 2009-2013 is the highest when compared to that during the entire study period ( Figure  9a). The precipitation variation results in annual TWS change of 2.82 ± 1.54 mm/mm (Figure 10a), while other factors result in annual TWS change of 11.28 ± 113.96 mm/a. Among all the basins in China, XnRB has the highest correlation coefficient of 0.86 (Table 3, Figure 9i). It can be seen from Figure 9i that precipitation and TWS changes are consistent with a time lag of three months. In addition, the linear relationship between precipitation and TWS is highly significant with the slope of the linear fitting model of 4.29 ± 1.76 mm/mm (Figure 10i).

Spatial Distribution of the TWS Trends at Seasonal and Annual Scale
We observed distinctly different TWS changes across China ( Figure 11). The spatial pattern of the seasonal TWS trends indicated that TWS trends are similar with different changing magnitude (Figure 11a-d). In addition, we also observed increasing TWS trends in southeastern China and ShRB (Figures 6, 7 and 11). Table 4 lists the statistical information about the percentage of areas with decreasing TWS trends in each river basin at the seasonal scale and annual scale. Figure 11 and Table 4 show that the proportion of regions with decreasing and/or increasing TWS trends is similar in all the river basins in China. Specifically, the percentage of regions in YRB, SRB, ShRB, and PRB with decreasing TWS trends is less than 50% seasonally, and it is exactly the opposite for the rest of the river basins. Similar findings were also achieved by [57].  As a result, 56.7% of the regions across China are dominated by decreasing TWS trends in spring with a declining rate of from −30 mm/a to 0 mm/a (Figure 11a), while less than half of the regions in China are dominated by increasing TWS trends. However, we found spatial heterogeneity of TWS changes at finer spatial resolutions. The east part of the YRB is experiencing increasing TWS trend which accounts for 57.5% of the total area of the YRB, while the major areas of the SRB, the PRB and the ShRB are characterized by increasing TWS trends as well, and the percentages of the areas are 97.1%, 74.4% and 85.6%, respectively ( Figure 11a). In other words, more than half of the areas are dominated by decreasing TWS trends with varying declining rates in the HRB (100%), the YB (85.3%), the HuaiRB (77.8%), the LRB (94.7%), the XbRB (62.8%) and the XnRB (88.6%).
It can be seen from Figure 11b that the TWS trend is ranging from −30 mm/a to 30 mm/a in summer during 2003-2015, which is the same as the range of TWS trend in spring. However, the declining rate of the TWS in the HRB in summer is a little greater than that in spring (Figure 11b). The exploitation of the groundwater is the major cause behind the larger decreasing magnitude of the TWS since the agricultural irrigation excessively withdrawal groundwater from deep wells every year [58,59]. As shown in Figure 11c,d, the areas with increasing TWS trends in northern China shrink in autumn and winter, while the ones in southeast China expand in autumn and winter. More specifically, the areas with increasing TWS trends in the ShRB are 85.6% and 80.8% in spring and summer, respectively, and then it drops to 69.9% in autumn and 69.2% in winter. In addition, the areas with increasing TWS trend in the PRB are 74.4% and 55.1% in spring and summer, respectively, and then it rises to 97.4% in autumn and 100% in winter.

Spatial Distribution of the TWS Trends at Seasonal and Annual Scale
We observed distinctly different TWS changes across China ( Figure 11). The spatial pattern of the seasonal TWS trends indicated that TWS trends are similar with different changing magnitude (Figure 11a-d). In addition, we also observed increasing TWS trends in southeastern China and ShRB (Figures 6, 7 and 11). Table 4 lists the statistical information about the percentage of areas with decreasing TWS trends in each river basin at the seasonal scale and annual scale. Figure 11 and Table  4 show that the proportion of regions with decreasing and/or increasing TWS trends is similar in all the river basins in China. Specifically, the percentage of regions in YRB, SRB, ShRB, and PRB with decreasing TWS trends is less than 50% seasonally, and it is exactly the opposite for the rest of the river basins. Similar findings were also achieved by [57].  Table 4. The percentage areas in decreasing trend of each basin in seasons and year with L denotes less than 50% and G denotes greater than 50%.   Table 4. The percentage areas in decreasing trend of each basin in seasons and year with L denotes less than 50% and G denotes greater than 50%.

Basin YRB SRB HRB YB HuaiRB LRB ShRB XbRB XnRB PRB CHN
As a result, 56.7% of the regions across China are dominated by decreasing TWS trends in spring with a declining rate of from −30 mm/a to 0 mm/a (Figure 11a), while less than half of the regions in China are dominated by increasing TWS trends. However, we found spatial heterogeneity of TWS changes at finer spatial resolutions. The east part of the YRB is experiencing increasing TWS trend which accounts for 57.5% of the total area of the YRB, while the major areas of the SRB, the PRB and the ShRB are characterized by increasing TWS trends as well, and the percentages of the areas are 97.1%, 74.4% and 85.6%, respectively (Figure 11a). In other words, more than half of the areas are

Influences of Climate Modes and Precipitation to Variations of TWS
The results from Figures 4 and 5 suggest that all the climate modes investigated in this study, namely, ENSO, IOD, NAO and PDO have influenced to variations of TWS. Meanwhile, linkages between precipitation and TWS can be seen from Figures 8 and 9. Previous studies have also testified about the relation between precipitation and TWS [60,61]. Voss et al. [19], Forootan et al. [62] and Forootan et al. [63], for example, demonstrated that a recurring decrease of seasonal precipitation brings about a decline of TWS. Hirschi et al. [64] found similar influence for 37 river basins in Europe, Asia, North America, and Australia. Other studies have indicated that large-scale climate modes, such as ENSO, IOD, NAO, Atlantic Multidecadal Oscillation (AMO) and PDO, had considerable effects on the fluctuation of precipitation over monsoon areas of China in a warming climate, which in turn would indirectly change the water storage [65][66][67][68].

Conclusions
In this study, we investigated the relationship between TWS and precipitation for 10 basins in China during 2003-2015 based on GRACE data in combination with precipitation data from GPCP, GPCC, and NOAA. The impacts of the climate drivers, such as ENSO, IOD, NAO and PDO on TWS variation in China were also identified and discussed.
The results showed that the variations in TWS are correlated with the large-scale climate variability since there are varying significant correlations between the climate indices and TWS across China. However, the influence of each climate model on water storage dynamics varies in regions of China with different time lags. Basically, the sphere of influence on TWS is relatively larger under the climate events of ENSO, IOD and PDO; and the discussion on the linear relationship between TWS change and precipitation variation provides the evidence that there is normally a good consistency between the two variables in most basins except XbRB. Then we further analyzed the TWS trend at the seasonal scale, and found that the characteristics of TWS trend are similar in spring and summer, while the results in autumn and winter have common magnitude and spatial pattern as well.
The current results present evidential information of the influence of climatic variability on China's water storage variations, and also exhibit the value of GRACE observations in the hydroclimatic study as a vital indicator of water storage dynamics, and therefore, can be helpful to the management of water resources.