Terrestrial Water Storage in China: Spatiotemporal Pattern and Driving Factors

Abstract

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.

1. 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 elevating frequency of floods [12]. Furthermore, the warming climate drives higher risks of floods and droughts, and evidence from model simulation results documents that the dry regions of the world will be drier, while the wet regions will become wetter [13,14]. This phenomenon was known as the rich-get-richer mechanism [15,16]. As one of the major water resources for human activities, such as agricultural irrigation, the storage of groundwater is decreasing globally under climate warming [17,18,19,20]. In this case, we attempted to use GRACE data to analyze terrestrial water storage (TWS) across China.

El Niño Southern Oscillation (ENSO) is an important subsystem of the global climate system. ENSO usually reaches its peak values in the northern hemisphere during winter and continues in the spring of the subsequent year and until the summer recession [21]. Although ENSO occurs in the equatorial Pacific Ocean, its impacts extend more than 75% of the Earth’s area [22,23]. A bunch of studies have shown that ENSO has a greater impact on climate anomalies in East Asia [24,25,26,27,28,29,30]. Meanwhile, most of the studies focused on the effects of spring and summer ENSO on precipitation in China, whereas Zhou et al. [31] indicated that ENSO has a significant positive impact on precipitation in winter (November-February) in southern China. Jiang et al. [32] also showed that ENSO has a good teleconnection with droughts and floods in the Yangtze River basin, China. Luo and Lau [33] suggested that ENSO exerts an amplifying effect on the heatwave activity in China. Meanwhile, the hydro-climatic conditions in East Asia are also affected by the Indian Ocean Dipole (IOD), North Atlantic Oscillation (NAO) and Pacific Decadal Oscillation (PDO) [34,35,36,37,38]. Besides, climate indices, such as ENSO were also proved to have impacts on terrestrial total water storage. Xie et al. [39] investigated the spatial patterns and temporal dynamics of the terrestrial total water storage across Australia under impacts of ENSO, IOD and Southern Annular Mode (SAM), and found multi-climate modes interactions and related impacts on Australia’s water storage dynamics. Vissa et al. [40] assessed the impacts of climate perturbations, such as ENSO for the groundwater storage changes and found that the decline of groundwater storage changes over northern India enhanced during the period of ENSO. However, no report is available now addressing the variability of total water resources in all the large basins across China and related driving factors.

Since 2002, Gravity Recovery and Climate Experiment (GRACE) satellite monitors global water storage, including groundwater and terrestrial water storage, at monthly to decadal time scales and has provided comprehensive data information on water cycle [3]. GRACE satellite observations have recorded monthly changes in earth surface gravity field that reflects the mass changes with unparalleled accuracy [41]. The major factor driving the variations of measured mass at monthly time scales is the water redistribution, and it allows GRACE to detect the changes of freshwater on the land surface. With the monitoring of GRACE satellite, the total amount of water, including snow, surface water, groundwater and soil moisture that enters or leaves an area each month can be weighted at an accuracy of 1.5 cm equivalent water height [3].

Therefore, the major objectives of this current study are: (1) To characterize variability and availability of TWS in both space and time; (2) to investigate attribution of TWS changes in both space and time from a perspective of climate indices, such as ENSO and so on; and (3) to quantify relationships between TWS changes and precipitation variations across China. This study will be of great importance in water resources management in a changing environment.

2. Materials

2.1. Study Region

China has complicated topography across its vast territory with various climates shifting from dry seasons and wet monsoons. In winter, cold and dry winds are from northern high-latitude areas; while in summer, warm and moist winds are from southern coastal areas at lower latitudes. The diverse climates in different regions of China are the results of highly complex topography and also the location of China, i.e., between the largest Eurasia and the largest Pacific Ocean. As a performance of diverse climates, there is a varied amount of precipitation falling in different regions (basins) of mainland China. Also, water resources are unevenly distributed in each regions (basins). Therefore, the analyze of TWS changes based on different regions (basins) will be done in this study. There are 10 large drainage basins in China (Figure 1). In this study, we took each river basin as the study unit to characterize relationships between TWS changes, climate indices and precipitation variations as well (

Table 1. Information of Basin Division.

Basin Label	Basin Name	Abbreviation	Area (×105 km2)	Pixel Numbers
a	Yangtze River	YRB	17.7	221
b	Southeast River	SRB	2.5	35
c	Hai River	HRB	2.7	51
d	Yellow River	YB	8.4	129
e	Huai River	HuaiRB	3.3	52
f	Liao River	LRB	3.7	57
g	Songhua River	ShRB	9.3	146
h	Xibei River	XbRB	29.1	444
i	Xinan River	XnRB	8.6	131
j	Pearl River	PRB	5.4	78

).

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

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

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

3. Methods

3.1. 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, $\widehat{g}\left(x\right)$, 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.

$\widehat{g}\left(x\right)$ 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 ${\lambda }_q\left(x\right)$ is the distance of the qth farthest $x_i$ from x, and W is the tricube weight function:

$$W\left(u\right)\{\begin{array}{c}{\left(1-u^3\right)}^3 for 0\le u<1\\ 0 for u\ge 1.\end{array}.$$

(1)

The neighborhood weight for any $x_i$ is

$$v_i\left(x\right)=W\left(\frac{\left|x_i-x\right|}{{\lambda }_q\left(x\right)}\right).$$

(2)

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\left(x\right)$ at ($x_i, y_i$). The value of the locally-fitted polynomial at x is $\widehat{g}\left(x\right)$. 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. ${\lambda }_n\left(x\right)$ is the distance from x to the farthest $x_i$. For q > n we define ${\lambda }_q\left(x\right)$ by:

$${\lambda }_q\left(x\right)={\lambda }_n\left(x\right)\frac{q}{n}.$$

(3)

Then we proceed as before in the definition of the neighborhood weights using this value of ${\lambda }_q\left(x\right).$

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 ${\sigma }^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\left(x\right)$ 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:

$$O_t=T_t+S_T+e_t,\mathrm{for}t=1,\dots ,\mathrm{n},$$

(4)

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

3.2. Statistical Methods

3.2.1. 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]:

$$\frac{dS}{dt}=P-ET-R,$$

(5)

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:

$$m=a\left(p-p_0\right)+b,$$

(6)

where m is the annul TWS change between two successive years (in mm), p is the annual precipitation (in mm/year), p₀ is the long-term annual average precipitation (in mm/year), a is the slope of the regression in (year⁻¹), 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.

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

3.2.3. 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:

$$r_j=\frac{{{\displaystyle \sum }}_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}},$$

(7)

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, $\overline{x}$ and $\overline{y}$ is the sample average of the climate indices and TWSA, respectively.

4. Results and Discussion

4.1. 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 increased continental average TWS during the period from 2003 to 2004 across China. Besides, we observed relatively high TWS during 2004–2005 when compared to the TWS for the whole period considered in this study. Afterwards, TWS decreased gradually until 2008, and the terrestrial water storage of China was relatively low in the late 2007 and early 2008. After 2008, TWS was in relatively unstable changes with its values fluctuating in large amplitude. Two peak values and two trough values in the fluctuation with peak values in 2010 and 2013, respectively, while trough values in 2011 and 2015, respectively.

4.2. Cross-Correlations of Climate Modes and TWS

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.

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 4 displays the maximum Pearson correlation coefficients between climatic drivers and TWS of China in 2003–2015, and the corresponding time lag in months between the study variables was also depicted in Figure 5. Based on the results in Figure 4 and Figure 5, we found that with varying strength of correlations and time lags at different regions in China, the four climatic drivers have close coupling relationship with TWS in the period of 2003–2015. In particular, MEI, PDO and DMI have even more widespread impacts on TWS than NAO (Figure 4).

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. Area percentage of negative and positive Pearson correlation coefficient © in basins of China.

Basin		YRB	SRB	HRB	YB	HuaiRB	LRB	ShRB	XbRB	XnRB	PRB	CHN
MEI-TWSr	<0	38.9%	65.7%	37.3%	33.3%	71.2%	45.6%	4.8%	30.5%	39.7%	65.4%	33.3%
	>0	47.5%	34.3%	62.7%	49.6%	28.8%	54.4%	86.3%	54.1%	43.5%	5.1%	53.2%
DMI-TWSr	<0	27.1%	2.9%	98.0%	76.0%	17.3%	84.2%	87.7%	45.2%	68.7%	0.0%	49.7%
	>0	59.3%	97.1%	0.0%	10.1%	34.6%	12.3%	1.4%	24.2%	20.6%	97.4%	31.6%
NAO-TWSr	<0	17.6%	40.0%	0.0%	6.2%	3.8%	54.4%	21.9%	38.5%	5.3%	34.6%	25.7%
	>0	31.2%	42.9%	0.0%	0.0%	0.0%	0.0%	0.0%	1.8%	26.7%	64.1%	14.0%
PDO-TWSr	<0	47.1%	82.9%	56.9%	39.5%	36.5%	64.9%	10.3%	35.1%	16.8%	55.1%	36.0%
	>0	30.8%	0.0%	41.2%	50.4%	55.8%	14.0%	71.2%	47.5%	75.6%	14.1%	47.1%

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

4.3. 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 (Figure 6 and Figure 7). The results indicated that YRB (Figure 6a), SRB (Figure 6b), XnRB (Figure 7d) and PRB (Figure 7e) are at the top of the river basin list considered in this study in terms of annual precipitation amount across China with annual precipitation amount of, respectively, 81 mm, 127.9 mm, 86.1 mm and 118.7 mm. However, among the top 4 river basins in precipitation, XnRB (Figure 7d) is the only river basin which is not the top 4 river basins in the annual amount of TWS in that XnRB (Figure 7d) is one of the arid regions in China. In addition, it is also the only river basin where TWS shows a decreasing trend from 2003 to 2015 by the declining rate of −6.2 ± 3.1 mm/a, while TWS within YRB (Figure 6a), SRB (Figure 6b), and PRB (Figure 7e) are persistently increasing. Apart from YRB (2.64 ± 1.92 mm/a), SRB (4.92 ± 1.56 mm/a) and PRB (4.56 ± 3 mm/a), ShRB (2.28 ± 1.44 mm/a) is the only river basin being characterized by increasing TWS trend (Figure 7b).

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. The highest correlation coefficient between changes in TWS and precipitation and the corresponding shifted a month.

Basin	YRB	SRB	HRB	YB	HuaiRB	LRB	ShRB	XbRB	XnRB	PRB
Coefficient	0.8	0.83	0.84	0.66	0.56	0.7	0.7	0.43	0.86	0.76
Shifted month	0	3	4	5	9	2	3	6	3	4

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

4.4. 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 (Figure 6, Figure 7 and Figure 11).

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
Spring	L	L	G	G	G	G	L	G	G	L	G
Summer	L	L	G	G	G	G	L	G	G	L	G
Autumn	L	L	G	G	G	G	L	G	G	L	G
Winter	L	L	G	G	G	G	L	G	G	L	G
Year	L	L	G	G	G	G	L	G	G	L	G

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.

4.5. Influences of Climate Modes and Precipitation to Variations of TWS

The results from Figure 4 and Figure 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 Figure 8 and Figure 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].

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