Divergent Trends of Water Storage Observed via Gravity Satellite across Distinct Areas in China

: Knowledge of the spatiotemporal variations of terrestrial water storage (TWS) is critical for the sustainable management of water resources in China. However, this knowledge has not been quantiﬁed and compared for the di ﬀ erent climate types and underlying surface characteristics. Here, we present observational evidence for the spatiotemporal dynamics of water storage based on the products from the Gravity Recovery and Climate Experiment (GRACE) and the Global Land Data Assimilation System (GLDAS) in China over 2003–2016. Our results were the following: (1) gravity satellite dataset showed divergent trends of TWS across distinct areas due to human factors and climate factors. The overall changing trend of water storage is that the north experiences a loss of water and the south gains in water, which aggravates the uneven spatial distribution of water resources in China. (2) In the eastern monsoon area, the depletion of water storage in North China (NC) was found to be mostly due to anthropogenic disturbance through groundwater pumping in plain areas. However, precipitation was shown to be a key driver for the increase of water storage in South China (SC). Increasing precipitation in SC was linked to atmospheric circulation enhancement and Paciﬁc Ocean warming, meaning an unrecognized teleconnection between circulation anomalies and water storage. (3) At high altitudes in the west, the change of water storage was a ﬀ ected by the melting of ice and snow due to the rising temperatures, yet the topography determines the trend of water storage. We found that the mountainous terrain led to the loss of water storage in Tianshan Mountain (TSM), while the closed basin topography gathered the melted water in the interior of the Tibetan Plateau (ITP). This study highlights the impacts of the local climate and topography on terrestrial water storage, and has reference value for the government and the public to address the crisis of water resources in China.


Introduction
Water shortage is undoubtedly a global crisis [1]. Terrestrial water storage (TWS) is an important index to evaluate water deficit, which includes surface water, groundwater, soil water, ice and snow, canopy water, and organism water [2]. Under the background of global warming, the intensification of water circulation will lead to the imbalance of TWS among regions [3], which threaten human production and domestic water supply. Urbanization and population explosion lead to the rapid increase of irrigation water, domestic water, and manufacturing water demand, which also aggravates the shortage of water resources. Since the 21st century, TWS dynamics have changed greatly due to the interference of human factors and climate factors [4].
The traditional way to estimate the TWS is the water balance equation (that is precipitation-evapotranspiration-runoff) at a small basin scale. However, this is extremely challenging on a national scale [5]. First, the spatial distribution of the observation stations is sparse. Secondly, the data collected by the stations have observation deviation, which is easily affected by the surrounding environment of the stations. Earth system and conceptual model (such as Global Land Data Assimilation System (GLDAS) [6], WaterGap [7], and PCRaster GLOBal Water Balance model [8]) estimates TWS by simulating water cycle process, which is another important way to obtain the spatiotemporal dynamics of TWS. However, the TWS estimated by the models are limited by the physical mechanism and the input parameters [9,10], and the TWS output by the model is often inaccurate in reflecting the real water storage change. Therefore, there is an urgent need for an efficient and accurate pathway to describe the space-time dynamics of TWS.
The Gravity Recovery and Climate Experiment (GRACE) satellite, successfully launched in March 2002, provides an unprecedented opportunity for monitoring TWS at the global, national, and basin scales [11]. The equivalent water height product driven by GRACE overcomes the shortcomings of traditional ground observation data and avoids the operation error of the hydrological model, and can directly monitor the change of TWS at a large-scale via measuring the change of temporal Earth's gravity field [12,13]. At present, GRACE has many applications in hydrology and climate variation. For example,  used GRACE satellite data to find that the TWS of global basin topography was lost in a large area, and that the rate exceeded 100 Gt year −1 [14]. They believed that the evaporation enhancement of the basin topography, frequent drought, and human intervention jointly lead to the loss of water [14]. Asoka et al. (2017) found that the strength of the Indian monsoon affected the intensity and frequency of precipitation, and directly or indirectly drove the change of TWS in most areas of India [15]. Recent reports indicated that 40% of the global TWS are affected by vegetation greenness, and vegetation is an important reason for affecting water storage [16]. Different climate and topographies will have a great impact on TWS [17]. However, the current research lacks in the exploration the effects of different climatic conditions and underlying surface attributes on the spatial differentiation of TWS.
China has abundant water resources. However, due to the large population, diverse climate, and complex landform, the contradiction between the supply and demand of water resources in local areas is very prominent [18,19]. The water resources bulletin shows that China's available water per capita is approximately 2200 m 3 , which is only a quarter of the world average [20]. This is also unevenly distributed in space. The southern area accounts for 81% of the total water resources in China, while the northern area only accounts for 19% [21]. To better understand the current situation of China's water resources, previous studies used GRACE gravity satellite data to analyze the spatial change information of China's water storage. However, the majority of the current research is limited to individual river basins and areas, such as the Yangtze River Basin [22], the North China Plain [23], and the Qinghai-Tibet Plateau [24]. However, the relationship between the space-time dynamics of TWS and climate change across distinct areas in China is still unclear. The current research is mainly limited to the impact of climatic factors on TWS (such as precipitation and temperature). Therefore, it is necessary to analyze the relationship between TWS and teleconnection factors, which can deepen the understanding of the causes of TWS spatiotemporal changes.
Recent reports indicated that the change of TWS in China is mainly due to monsoon and teleconnection, such as El Nino Southern Oscillation (ENSO) and East Asian monsoon [25,26]. He et al. (2020) found that ENSO has a more significant impact on precipitation and TWS in East China, where China's social economy is developed [27]. Zhang et al. (2015) found that there was a time delay effect between TWS and ENSO in the Yangtze River Basin [28]. Different time delays may be attributed to different underlying surface properties, including the water recharge process and topography, which indicates that the response time of TWS in China to ENSO is divergent in different areas [28]. Han et al. (2019) investigated that Arctic Oscillation (AO) and ENSO have a great influence on the TWS in Yunnan Province [29], in particular when ENSO is in the positive (negative) phase, the TWS are in a deficit (abundant) phase. Zhang et al. (2018) showed that drought caused by circulation anomaly made the TWS significantly lower than in normal years [30], which indicates that a climate anomaly will have a far-reaching impact on TWS, and the impact will be more severe in the future [31]. The above studies have discussed the possible effects of climate anomalies on the temporal and spatial pattern of TWS in China. However, the analysis of the impacts of multiple teleconnection factors on water storage in different regions of China has not been discussed.
Large-scale and long-time series dynamic monitoring of TWS is helpful to improve the accuracy level of water resources management. In this paper, the TWS driven by GRACE equivalent water height (2003-2016, monthly) products were taken as indicators to analyze the spatiotemporal variation of water storage comprehensively and quantitatively over China. Specifically, we attempted to investigate three areas: (1) whether the TWS trends show consistent patterns between model simulation and satellite observation; (2) the spatial evolution characteristics of TWS at different scales; (3) to explore the main drivers affecting the variations of TWS across hotspot areas.

Study Area
China's terrain from west to east has a ladder-like distribution according to altitude (Figure 1). Affected by the topography, many rivers originate from the mountains and plateaus in the west, and flow eastward and into the ocean. Overall, the spatial distribution of water resources in China was extremely uneven. The south of China in the monsoon zone was rich in rainfall and runoff, and had relatively large water storage. The north of China was also a monsoon region, but there was less precipitation and a large population, and the problem of water shortage was particularly prominent. In Northwest China, precipitation was scarce, but evaporation was strong, and social-economic development was restricted by water resources. The Tibet Plateau was an Asian water tower [32], where many large rivers in Asia develop. The water storage was abundant, but most of it is frozen on the plateau in the form of glaciers and snow.
In this study, based on geographical location, climate type, and underlying surface, we selected South China (SC), North China (NC), the Interior of the Tibetan Plateau (ITP), and Tianshan Mountains (TSM) as hotspot areas ( Figure 1). The detailed rules for selecting hotspots were as follows. First, whether it is monsoon has an enormous impact on water storage. SC and NC were chosen for famous monsoon regions, while the ITP and TSM were inland far from the sea. The climatic factors (heat and humidity) on water storage should be considered. NC and SC were representative of warm-wet and warm-dry monsoon climate, respectively, but ITP and TSM are typical mountain climate with temperature limitation. Secondly, the potential impact of the underlying surface on water storage was considered. NC was mainly farmland ecosystems on plains, SC was forest ecosystems on hills, and TSM and ITP were alpine ecosystems. In particular, ITP had a closed basin topography, while the TSM had a mountain topography (water can flow out). Thirdly, compared to NC and SC, there was little contribution of human activity to the region's water storage in ITP and TSM.

TWS Dataset from Gravity Satellite
We used newly released GRACE-based (Gravity Recovery and Climate Experiment) TWS (terrestrial water storage) products from the NASA Jet Propulsion Laboratory at 0.5 • spatial resolution (monthly, 2003-2016) [33]. The unit of TWS is the equivalent water thickness, which represent the deviations of mass in terms of the vertical extent of water [34]. TWS has been widely applied in hydrological research to represent terrestrial water storage changes [12,33], including groundwater, soil water, canopy water, and channel water, in terms of anomalies relative to the average equivalent water height value during 2003-2009. More detailed algorithms, theories, and acquisition of the TWS can be found in previous studies [35,36]. In this paper, TWS datasets were comprehensively evaluated for water storage changes in the space-time level covering the whole of China.

Climate and Hydrological Dataset
The dataset used in this paper included three parts (1) Soil moisture storage (SMS, monthly, 0.25 • spatial resolution) was derived from the Global Land Data Assimilation System (GLDAS) soil moisture product [6]. The GLADS Noah hydrological model outputs four layers of soil moisture dataset for 0-10, 10-40, 40-100, and 100-200 cm, the sum of different thicknesses was regarded as the soil moisture storage at the pixel level for the same period.
In addition, the time series of TWS, SMS, PRE, TEM, and ET were then resampled to generate a 0.25 • dataset in this study for further analysis.

Linear Regression Models
The trends of water storage from 2003 to 2016 at the pixel scale were calculated by linear regression methods, and the t-test at the 5% was defined as significant. The TWS trend of 34 provinces and 10 basins (Supplementary Table S1 and S2) in China were also calculated.
In this formula, the slope is the one-variable linear regression equation fitted by the variable to indicate the interannual trend of the indicator; i is the time variable, and n is the year of annual time series. When the slope < 0, it means decreases with time; when the slope > 0, this indicates the variable increases with time; the larger the absolute value of slope, the larger the trend of the variable.
The multiple linear regression also was applied to calculate the most important environmental factors between TWS and other variables in the time series.
where TWS is the terrestrial water storage; X1 is the precipitation; X2 is the average temperature; and X3 is the evapotranspiration. The most important factor is the environmental factor corresponding to the maximum value in the determination coefficient (a, b, and, c).

Breakpoint Algorithm
To detect and characterize changes in water storage, the BFAST (Breaks For Additive Season and Trend) algorithm was used on TWS monthly time series. BFAST decomposed the original time series into three parts: trend, seasonal, and residual components by using a piecewise linear trend and seasonal components [37]. In this paper, assuming 30 months as a moving data window, the bandwidth parameter is set to 0.14, and the significance level was set to 0.05 in the BFAST analysis. The 'bfast01 function in the "bfast" package in the R software was used to find the most important breakpoint. The advantage of the bfast01 algorithm is not only the ability to detect breakpoint in time series, but also the ability to detecting the characteristics [38]. In addition, the breakpoint characteristics of the TWS time series were divided into four categories according to the previous studies [39].

Cross-Wavelet Transformation
In order to judge the teleconnection between TWS and the atmospheric circulation index as well as sunspots, the Morlet wavelet was used to analyze the correlation between the monthly time series x (water storage) and y (PDO, ENSO, AO, and Sunspots) in the time domain and frequency domain. The cross-wavelet energy spectrum can show the same energy spectrum area of two time series after the wavelet transform, which reflects the structural characteristics of two time series after wavelet transform in the time-frequency domain [40].

Empirical Orthogonal Function Analysis
Empirical orthogonal function (EOF) is a method to analyze the structural characteristics and obtain information of the hydrological data on time series [41]. The eigenvector corresponds to the spatial sample, and also known as spatial eigenvectors or spatial modes (EOF), which reflects the spatial distribution characteristics of the element field to a certain extent; the principal component corresponds to the time variation (PC), reflecting the weight change of the corresponding spatial mode with time. In this paper, empirical orthogonal function (EOF) analysis was applied to assess the role of large scale climate anomaly on GRACE-based TWS during 2003-2016. Considering that the dominant mode can explain the largest variance and can best reflect the spatial characteristics of TWS, the dominant mode (EOF-1) and the corresponding principal component (PC-1) were calculated.

Spatiotemporal Characteristics of Terrestrial Water Storage
Linear regression was performed for 2003-2016 using the equivalent water height (terrestrial water storage, TWS) from the GRACE satellite and soil moisture storage (SMS) from the GLDAS dataset to assess and compare the spatial dynamics of water storage ( Figure 2). The trend of GRACE-based TWS and GLDAS SMS was −0.93 mm per year and 0.04 mm per year across the entire country, respectively. Two-fifths of the area was markedly declining and 30% was markedly increasing from the GRACE product (Figure 2a), and the soil moisture storage showed less decreasing and increasing (19% and 27%, respectively) over China (Figure 2b). The annual changes in the water storage anomaly, including total and soil, showed strong spatial heterogeneity. The decrease was prominently clustered in one area-most notably in Northern China, whereas increase was South China. More than 40% of the trend between TWS and SMS showed a positive correlation, and only 11% showed a negative correlation (Figure 2c). The correlation of the TWS trend and SMS trend of each province and basin were found to be 0.77 and 0.81, respectively, which indicated that the trends of TWS and SMS in most provinces and basins were consistent (Figure 3a,b). Linear regression and correlation analysis jointly proved that the trends of TWS and SMS were consistent in space at the pixel, province, and basin scale.
The spatial patterns of TWS trends at the province and basin levels are shown in Figure 4. Overall, TWS had clear spatial differentiation across distinct areas in China. At the provincial level, the majority of the provinces in northern China had a negative trend of TWS, and the loss rate of TWS in Hebei, Shandong, Henan, Anhui, Tianjin, and Beijing was more than 10 mm per year. On the contrary, in most of the southern provinces, the TWS was increasing, and Jiangxi, Zhejiang, Fujian, and Guizhou exceeded 10 mm per year (Figure 4a). Compared to the TWS changes at the basin level, in which the loss rate of water storage in the Haihe River and Huai River exceeded 10 mm per year, the gain rate was mainly in the Pearl and Yangtze River (Figure 4b).
Abrupt change analysis revealed that the variation of water storage in China was diverse ( Figure 5). The dominant types of breakpoints in northern China were 'increase to decrease', whereas in southern China were 'increase and increase'. The major breakpoint types in western China were 'increase to decrease', and 'decrease to decrease' (Figure 5a). The average year of breakpoints on the TWS time series was 2011, and the pixels with abrupt changes from 2010 to 2014 accounted for 48% of the total pixels in China (Figure 5b). The majority of the pixels had significant abrupt changes, and non-significant (p < 0.05) changes were only 2% in China (Figure 5c). In conclusion, the results of the breakpoint algorithm also show that TWS in the north was decreasing, and in the south was increasing. Figure 6 illustrates that the GRACE-based TWS showed remarkably divergent trends in four key areas over the years. In the monsoon areas, the SC water storage decreased steadily from 2000 to 2015 and recovered slightly in 2016 (Figure 6a). However, the NC water storage showed a certain increase over the past decade and reached the highest in 2016 (Figure 6b). For mountain areas, the average trend for water storage decreased significantly in TSM, even in 2016 with weak growth (Figure 6c). In contrast, the ITP water storage increased significantly from 2000 to 2012 and reached a relatively high level (Figure 6c).

Potential Factors Affecting Total Water Storage
According to the slopes of linear regression methods, the environment is an important factor contributing to the variations of TWS (Figure 7). For precipitation, approximately 14% of the pixels were positively increased (slope > 0 and p < 0.05) for the period 2003-2016, and most of the significant trends of water storage were in South China (Figure 7a). About 53% were changed for temperature, and the proportions of significant positive and negative correlation was 46% and 7%, respectively. The positive increasing pixels were mainly concentrated in the Tibetan Plateau and Northwest China (Figure 7b). The linear trend of evapotranspiration in South China increased significantly (Figure 7c). On the whole, precipitation and evapotranspiration were the main environmental factors affecting TWS over China, accounting for 34% and 47%, respectively. Precipitation was mainly controlled in the south and Tibetan Plateau, while evapotranspiration was mainly controlled in the north. The temperature was also an important contributing factor affecting the Tibetan Plateau (Figure 7d). Telecorrelation factors can further affect regional and global precipitation anomalies through global atmospheric circulation anomalies, resulting in the redistribution of water storage in time and space [42]. The results of TWS telecorrelation analysis showed that atmospheric circulation anomalies had an impact on the NC, SC, TSM, and ITP areas, and the impact intensity was significantly different (Figure 8). Except for the TSM, the effects of various teleconnection factors on TWS had a cycle of 8-16 months, indicating the annual cycle change of TWS across regions. Among them, the TWS of NC had significant negative and positive correlations with PDO and ENSO in the entire sequence, indicating that the change in the positive and negative phases of PDO and ENSO was the main factor driving water storage (Figure 8a). The SC and each teleconnection factor demonstrated a significant correlation, as southern China is a typical monsoon climate area, and is close to the ocean, which is seriously affected by Western Pacific subtropical highs (Figure 8b). As the TSM is located inland, the water storage was less affected by the precipitation factor, and thus the telecorrelation factor showed insignificant influence (Figure 8c). In addition, due to its high altitude and low latitude, ITP was less affected by AO (Figure 8d), while the significant negatively correlated changes shown by PDO and ENSO may be affected by the warm-humid airflow in the Bay of Bengal [29]. In general, the monsoon regions (SC and NC) were more sensitive to factors such as precipitation than the non-monsoon regions (TM and ITP), therefore, they showed more significant changes with the teleconnection factors. To further appreciate the importance of atmospheric circulation caused by sea surface temperature (SST) anomalies on TWS, we then compared the teleconnection between water storage and SST in two monsoon areas (Figure 9). The spatial structure of the dominant mode (EOF-1) obtained from EOF analysis was very similar to water storage trends for 2003-2016, which explained approximately 89% and 85% of the total squared covariance for NC and SC, respectively (Figure 9a,b). The principal component of EOF-1 corresponding to PC-1 showed that water storage in NC and SC decreased and increased continuously, respectively (Figure 9c,d), which indicated that the dominant mode represented the change rate in the water storage anomaly. Correlation analysis indicated the PC-1 and SST in the Pacific Ocean had a negative relation in NC, but had a positive relation in SC (Figure 9e,f). In particular, the water storage gain in SC was highly correlated with the rise of SST (R > 0.5).

Water Storage Changes Aggravate the Spatial Heterogeneity of Water Resources
Here, we provide a comprehensive national-scale analysis of spatiotemporal variation in GRACE-based TWS during 2003-2016 in China. We demonstrate that the spatial pattern of TWS changes in China is extremely unbalanced and that geographical differentiation is prominent. Trend analysis and breakpoint analysis jointly revealed that the most general shift characterization of the space-time pattern is that TWS rapidly and continuously decreased in the North, while the TWS widely increased in the South. These results are consistent with previous studies [19,43].
In general, the water resource pattern of China is that the north and northwest are relatively arid due to the limitation of the precipitation supply, while the South has very abundant precipitation. We found that the shifting trend of TWS was consistent with the spatial distribution of water resources in China. Concerning the relation between the TWS trend and the annual precipitation (Figure 10a,b), we noted that the water storage gain at the provinces and basins level was more obvious with high-precipitation (in the South), while the water storage loss was very serious with low-precipitation (in the North, Northwest, and Tibet). That is to say, the water storage in the arid areas decreased, but in the humid areas increased. Not only was the precipitation less for the northern and northwestern provinces, but the evapotranspiration was strong. This asymmetry of precipitation and evapotranspiration leads to the water balance being easily disrupted. In the past 14 years, the water storage loss directly led to the breaking of the water circulation system, which had less water supply and further aggravated the crisis of water resources in semihumid and arid areas. On the contrary, the annual precipitation in southern provinces was mostly above 800 mm, and water resources were not the main factor restricting sustainable development. The continuous rise of water storage in South China only reflects the enhanced the probability of heavy rainfall and flood. Consequently, we propose that divergent change rates in water storage aggravate the contradiction of uneven spatial distribution of water resources. By comparing SMS, we further evaluated the accuracy and regional applicability of the TWS dataset from the GRACE satellite. Our results indicate that TWS changes were generally consistent with SMS simulated from the GLDAS model, and those were highly correlated at the provincial and watershed levels. The cross-test proved that the water storage observed by gravity satellite was an efficient and reliable data supplement, which provided a way to monitor and manage water storage on a large scale. Although some studies indicated that SMS is a traditional and widely used indicator of water storage [44], it was obvious that SMS was not sufficient to reflect the changing pattern of water storage on a large scale. For example, Rodell et al. (2009) reported that the Indian subcontinent experienced rapid depletion of water storage due to excessive pumping over the past decade [1], GRACE-based TWS sensitively observed a decline in water storage, while SMS did not detect this. Since the water storage provided by the GLDAS hydrological model is only the soil moisture content of the 0-200 cm soil layer, it cannot capture the groundwater variations [5]. However, the GRACE water storage signal comes from gravity changes of all matter, which means that any form of water (groundwater, soil water, surface water, and ice) on terrestrial can be monitored. We found that the spatial matching between TWS and SMS was good in eastern China, but not in the western region. This is due to the fact that, compared with the western region, the eastern region had abundant precipitation and a higher soil moisture, which accounted for a larger proportion of terrestrial water storage [29]; thus TWS and SMS demonstrated a similar trends. However, in the western part of the terrestrial water storage, groundwater accounted for a larger proportion [45], and precipitation and runoff were lower. Considering that GLDAS can only represent the soil water composition, this may be the reason for the TWS mismatch in the eastern and western regions.

The Dominant Driver is Different across East Asian Monsoon Areas
We next investigated what leading control factors explained the divergent trends of TWS from 2000 to 2016 over the monsoon regions. As NC and SC are located in the core area of human activities, climate anomalies and human disturbances affect the water storage of NC and SC. However, we found that the dominant factors controlling the change of TWS were different.
On one hand, atmospheric circulation controls the global water vapor transport and distribution [42], and are generally considered to be the key factors affecting the spatiotemporal dynamics of precipitation [5,29]. Han et al. (2019) proposed that a single circulation index is not sufficient to explain changes in TWS [29]. We used multiple circulation indices and sunspots to jointly reveal that the TWS in eastern China was indeed affected and controlled by climate anomalies. However, Han's research was limited to the research in Yunnan Province [29], ignoring that telecorrelation factors such as ENSO and AO are atmospheric oscillations on a global scale, and ignoring that the impact of TWS was different in distinct regions. Our results showed that the influence of atmospheric circulation and sunspots in NC was weak, but was very strong in SC (mainly PDO and ENSO). On the other hand, an important characteristic of precipitation in the East Asian monsoon area is to be easily affected by the Pacific surface temperature [17,46]. Recently it was reported that the change of TWS in eastern China was directly affected by the precipitation brought by the East Asian summer monsoon [27]. According to previous steps [15], here we continued to explore the impact of ocean warming on TWS in the monsoon region by applying the empirical orthogonal function and teleconnection analysis. We found that the decline of NC water storage was weakly negatively correlated with the rise of SST, while the increase of SC water storage had a very significant positive correlation the rise of SST. Therefore, we believe that the sea temperature anomaly may be another important teleconnection factor that affects TWS. Recently studies claimed that TWS distribution pattern in China was potentially linked to regional changes of the Hadley-type meridional circulation which aggravated the unevenness between NC and SC water storage distributions over the monsoon regions [27].
In the past few decades, China's grain production and population demonstrated a rapid upward trend, and NC is one of the most concentrated areas of farmland and population in China. Due to the needs of socioeconomic development, human beings have overexploited the groundwater resources over the North China Plain [47,48]. Crop planting and urban expansion require more water resources; however, with less precipitation it is difficult to meet the supply. Pumping groundwater has become the most efficient pathway to obtain water with the lowest cost. The previous study analyzed the change rate of the groundwater level for NC using the observation well data, and found that the groundwater of NC was being lost at a rate of −6 mm year −1 [23]. This confirms that the depletion of TWS observed by the GRACE satellite may be caused by excessive groundwater extraction in this area. Moreover, vegetation can also lead to TWS anomalies [16]. For the North China Plain, groundwater is the main source of water for crop growth [23], because higher vegetation greenness in farmland means that more groundwater resources need to be pumped. The global surface water storage accounts for 36.08% of the TWS, which is an important part of the water cycle [14]. The surface water resources of NC also showed a decreasing trend over a long time [19], we propose that the change of open water body area may be another reason for the decline in TWS.
The TWS of SC showed a significant increase, which may be more closely related to the precipitation. The atmospheric circulation was previously shown to alter the large-scale water supply and distribution [42]. We indicated that the enhancement of atmospheric circulation and the rise of SST made SC precipitation increase continuously over the past 10 years. This is consistent with Zhang's [28] results in which they claimed that an ENSO and water storage were significantly correlated in SC, and that a La Niña (El Niño) event was associated with low (high) TWS. Ocean warming and increased water vapor transport change regional hydrological factors, which makes water resources redistribute in time-space, and then leads to an upward trend of TWS of SC. A recent study reported an increasing trend in water storage of SC can be confirmed by the vegetation greening via ecological engineering [49]. On the one hand, the increasing precipitation means that the Southern Hills collect more water and promote growth. On the other hand, the forest ecosystem has a strong ability to absorb and conserve water, and vegetation restoration can also maintain more water in the soil and groundwater [50]. Therefore, the increasing water supply and the underlying surface of the forest ecosystem provided a clear SC water storage increase during 2003-2016. In particular, some large-scale water conservancy projects may also lead to the surplus of regional water storage [44]; therefore, the construction of reservoirs and wetlands is also a reason for the obvious raise of TWS. To sum up, the decreases in TWS in NC were mainly driven by human activities, while the decreases in SC were mainly driven by climate factors.

The Topography is the Key in High Elevation Areas
According to the teleconnection analysis, the influence of atmospheric circulation and sunspot on TWS of TSM and ITP was relatively weak, which indicates that precipitation may not be the dominant factor to impact TWS in high elevation areas. TSM and ITP are barren land, and the impact of human activities on TWS in no man's land is very weak [24,51]. The rise of global temperature is bound to lead to glacier retreat and snow melting [32], which will affect the TWS dynamics on the TSM and ITP. Therefore, the TSM and ITP may be more strongly affected by temperature. However, our study found that the change rate of TWS on TSM and ITP was completely the opposite.
The loss of TWS for TSM may be mainly related to the rapidly rising temperature. Recent studies have reported that the average temperature of TSM was increasing at the rate of 0.38 • C 10 −1 year −1 [52], which makes the mountain glaciers and snow melt when the temperature exceeds 0 • C. The decrease of TWS of TSM indicated that the solid ice and snow in the high mountain area was converted into liquid or gaseous water [51]. Liquid water can easily enter the runoff and soil, and a part of the soil water seeps into the groundwater. Driven by gravity, this water flows from high-altitude mountainous areas to low-altitude farmland and grassland ecosystems, and so the TWS of TSM in mountainous terrain decreased rapidly in the past 10 years.
However, compared with TSM, ITP has a very unique topography, surrounded by high mountains, which is a completely closed basin. The gradually rising global temperature leads to the melting of ice-snow on the Tibetan Plateau [32]. First, as a depression basin on the Tibetan Plateau, the glacier and snow melted by the nearby high mountains can be gathered [53], which makes the TWS of the ITP increase greatly. Secondly, the increase of the number of lakes and the lake area on the ITP represents an obvious increase of surface water storage [54], which is a direct reflection of the gain of TWS. Finally, the altitude of the ITP area is above 4000 m, and the evaporation is very weak. Water resources are not easily lost through evaporation [55], which may be the reason why TWS can gather and maintain growth for a long time. Here we presented that rising temperature was the main reason for the change of TWS for TSM and ITP; however, the opposite trend of TWS was caused by different topographies.
In summary, the melting ice-snow flows out of the TSM through the runoff and soil due to the mountain terrain, while the water resources of glacier retreat and snow melting were gathered in the closed basin.

Conclusions
As an important indicator of the earth system, TWS is the summation of surface water, groundwater, soil water, glacier, and snow cover. Based on the GRACE-based terrestrial water storage, we successfully applied multiple methods to evaluate the spatiotemporal trends and associated drivers of TWS across distinct areas in China. Compared with GLDAS soil water storage, we propose that the equivalent water height from gravity satellite is a better indicator to describe large-scale water storage. During our study period of 2003 to 2016, one-third of the water storage variations presented a considerable surplus, that is, gaining more water resource, while two-fifths of the water storage gradually depleted, that is a loss. We found that the gain and loss of water storage are generally associated with the south and the north, respectively, which likely aggravated the contradiction of water resource allocation over China. The most prominent characteristic was that the water storage of the southern provinces and basins increased with more precipitation, while the water storage of the northern provinces decreased with less precipitation. We concluded that there were divergent trends of water storage due to human factors and climatic factors across distinct areas over the entirety of China. We combined the changing trend of water storage, climate type, and topography to identify four hotspot areas, and then explored the dominant driver of water storage variation. In the monsoon area, the trends of water storage in NC and SC are the opposite. The loss of water storage in the north was mainly caused by excessive pumping of groundwater, while the surplus of water storage in the south was mainly due to the increase of rainfall, which was caused by the enhancement of large-scale atmospheric circulation and the rise of ocean surface temperature. In high mountain areas, the trend of water storage was found to be controlled by the melting of ice-snow caused by rising temperatures. However, different topographies led to inconsistent trends of water storage. The TSM is a mountainous topography, and the melted water under gravity flows out of the mountain, and thereby the water storage is reduced. However, the ITP has a closed basin topography, high around the exterior and low in the interior, which absorbs the melting water from the surrounding mountains and leads to increased water reserves. We suggest that more measures should be taken to adjust the spatial imbalance of water storage in the future. Inter-basin water transfer and water-saving irrigation may be effective measures to alleviate the water shortage in China.