GRACE-Based Terrestrial Water Storage in Northwest China : Changes and Causes

Monitoring variations in terrestrial water storage (TWS) is of great significance for the management of water resources. However, it remains a challenge to continuously monitor TWS variations using in situ observations and hydrological models because of a limited number of gauge stations and the complicated spatial distribution characteristics of TWS. In contrast, the Gravity Recovery and Climate Experiment (GRACE) could overcome the aforementioned restrictions, providing a new reliable means of observing TWS variation. Thus, GRACE was employed to investigate TWS variations in Northwest China (NWC) between April 2002 and March 2016. Unlike previous studies, we focused on the interactions of multiple climatic and vegetational factors, and their combined effects on TWS variation. In addition, we also analyzed the relationship between TWS variations and socioeconomic water consumption. The results indicated that (i) TWS had obvious seasonal variations in NWC, and showed significant decreasing trends in most parts of NWC at the 95% confidence level; (ii) decreasing sunshine duration and wind speed resulted in an increase in TWS in Qinghai province, whereas the increasing air temperature, ameliorative vegetational coverage, and excessive groundwater withdrawal jointly led to a decrease in TWS in the other provinces of NWC; (iii) TWS variations in NWC had a good correlation with water storage variations in cascade reservoirs of the upper Yellow River; and (iv) the overall interactions between multiple climatic and vegetational factors were obvious, and the strong effects of some climatic and vegetational factors could mask the weak influences of other factors in TWS variations in NWC. Hence, it is necessary to focus on the interactions of multiple factors and their combined effects on TWS variations when exploring the causes of TWS variations.


Introduction
Terrestrial water storage (TWS), including surface water, groundwater, snow cover, and soil moisture, not only plays an important role in the global hydrological cycle, but also provides abundant freshwater resources for human society and the terrestrial ecological system [1][2][3].With a changing environment (climate change and human activities), TWS has obvious variations in many river basins and/or regions, which bring big challenges to local water resource management and utilization [3][4][5].It is imperative to continually monitor TWS variations.However, it is impossible to do so using in situ observations and hydrological models due to a limited number of gauge stations and the complicated spatial distribution characteristics of TWS [6][7][8].Accordingly, a more advanced technology is hoped to effectively tackle this tricky issue.
To observe spatial and temporal changes of Earth's gravity field, the National Aeronautics and Space Administration (NASA) and the German Aerospace Center (DLR) jointly developed the Gravity Recovery and Climate Experiment (GRACE) satellite mission [9].Since the launch of the GRACE satellite mission in March 2002, an increasing number of hydrologists realized that it could act as a new effective tool for measuring TWS variations.GRACE could continuously measure TWS variations through conversion of the gravitational field, breaking the limitations of spatial heterogeneity for TWS.Furthermore, many recent research results showed that GRACE is able to provide reliable data for an investigation of TWS in a large area on a monthly basis [10][11][12][13][14].
Until now, most studies were concerned with the spatiotemporal characteristics of TWS variations based on GRACE, paying less attention to the causes of TWS variations [3].To determine the causes of TWS variations, some scholars investigated the relationships between TWS and climatic and vegetational factors.Li et al. [15] found a good correlation between TWS and vegetational coverage in western and central Europe during droughts, on the basis of a severe deficiency in TWS harming vegetational growth to a large extent.Zhou et al. [16] suggested that TWS was dominated by precipitation in the Poyang Lake basin.Deng and Chen [3] showed the effects of precipitation and air temperature on TWS in Central Asia over the past decade.These studies substantially promoted the causal analysis of TWS variations.However, most previous studies only showed a correlation between GRACE-based TWS and one factor on a monthly basis without considering the combined effects and the interactions of multiple factors.This probably generated some unreliable correlations, maybe even misleading the causal analysis of TWS variations.For example, Deng and Chen [3] found that monthly TWS and monthly average air temperature (AT) had an obviously positive correlation in Xinjiang; however, they failed to clearly explain the effect of AT on TWS due to the combined effects of air temperature, precipitation, and wind speed on TWS.In contrast, this study paid more attention to the interactions of multiple factors and their combined effects on TWS variations.
Precipitation (P) and evapotranspiration (Et) directly influence TWS variation; however, it is difficult to access actual Et in a large region.According to the Penman-Monteith equation [17,18], Et is influenced by many climatic factors, such as sunshine duration (SD), air temperature (AT), and wind speed (WS).In addition, vegetational growth also has important effects on Et [19].The normalized difference vegetation index (NDVI), an index for the greenness of vegetation, could characterize vegetation growth [20].Thus, multiple climatic and vegetational factors including P, SD, AT, WS, and NDVI were adopted to account for TWS variations in this study.
In addition to climatic and vegetational factors, some studies also analyzed the relationship between TWS variations and socioeconomic water consumption.Rodell et al. [21] suggested that groundwater depletion was mainly attributable to irrigation in northwest India.Feng et al. [22] believed that a decrease in TWS was caused by agricultural water consumption in North China.Moreover, Deng and Chen [3] deemed that agricultural diversion and groundwater withdrawal were primarily attributable to the decrease in TWS in Central Asia.Artificial reservoirs play a significant role in utilizing and managing surface water resources [23]; however, few studies explore the relationship between TWS changes and reservoir operation.In addition to discussing the effects of groundwater withdrawals, this study also briefly analyzed the correlation between water storage variations in reservoirs and TWS variations.
Northwest China (NWC) is a typical arid and semi-arid region playing an important strategic role in the social and economic development of China [24][25][26].Since the implementation of the Belt and Road initiative of China, this region faces an unprecedented new opportunity for development [27].With rapid socioeconomic development, the water shortage problem became increasingly severe in NWC.To our knowledge, this study is the first study employing GRACE data to characterize TWS variations in NWC.The results of this study are proposed to hopefully help the management of local water resources in this region.
Two targets were established for this study.The first was to examine TWS variations in NWC, and to identify the relationships between TWS and climate change, vegetational coverage, and socioeconomic water consumption.The second was to show the interactions of multiple climatic and vegetational factors, and their combined effects on TWS variations.The paper is organized as follows: Section 2 describes the study area, data collection and processing, and two correlation analysis methods for hydrological and meteorological time series; Section 3 describes the detailed TWS variations in NWC, and an analysis of relationships between TWS and climate change, vegetational coverage, and socioeconomic water consumption; Section 4 presents the major conclusions of this study.

Study Area
NWC, shown in Figure 1, includes the Shaanxi, Ningxia, Gansu, Qinghai, and Xinjiang provinces.NWC is in the interior of China and is located far from oceans, having a total area of 3.11 million km 2 .Its terrain is very complex, consisting of mountains, plateaus, basins, and plains.Among the five provinces, the area of Xinjiang (1.66 million km 2 ) is the largest, and the average elevation of Qinghai (3000 m above sea level) is the highest.Most parts of NWC are mainly characterized by an arid or semi-arid climate, and the annual average precipitation in Shaanxi (634 mm) and Xinjiang (144 mm) are the maximum and minimum values, respectively.Overall, the ecosystem of NWC is rather vulnerable, particularly in Xinjiang where 64% of the area is desertified land as of the end of 2014 [28].The upstream of the Yellow River flows through the Qinghai, Gansu, Ningxia, and Shaanxi provinces.In the upper Yellow River, Longyangxia (LYX) and Liujiaxia (LJX) in the Qinghai and Gansu provinces, respectively, are two large cascade reservoirs with a total storage capability of 30.4 billion m 3 [29].
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 21 variations in NWC.The results of this study are proposed to hopefully help the management of local water resources in this region.Two targets were established for this study.The first was to examine TWS variations in NWC, and to identify the relationships between TWS and climate change, vegetational coverage, and socioeconomic water consumption.The second was to show the interactions of multiple climatic and vegetational factors, and their combined effects on TWS variations.The paper is organized as follows: Section 2 describes the study area, data collection and processing, and two correlation analysis methods for hydrological and meteorological time series; Section 3 describes the detailed TWS variations in NWC, and an analysis of relationships between TWS and climate change, vegetational coverage, and socioeconomic water consumption; Section 4 presents the major conclusions of this study.

Study Area
NWC, shown in Figure 1, includes the Shaanxi, Ningxia, Gansu, Qinghai, and Xinjiang provinces.NWC is in the interior of China and is located far from oceans, having a total area of 3.11 million km 2 .Its terrain is very complex, consisting of mountains, plateaus, basins, and plains.Among the five provinces, the area of Xinjiang (1.66 million km 2 ) is the largest, and the average elevation of Qinghai (3000 m above sea level) is the highest.Most parts of NWC are mainly characterized by an arid or semi-arid climate, and the annual average precipitation in Shaanxi (634 mm) and Xinjiang (144 mm) are the maximum and minimum values, respectively.Overall, the ecosystem of NWC is rather vulnerable, particularly in Xinjiang where 64% of the area is desertified land as of the end of 2014 [28].The upstream of the Yellow River flows through the Qinghai, Gansu, Ningxia, and Shaanxi provinces.In the upper Yellow River, Longyangxia (LYX) and Liujiaxia (LJX) in the Qinghai and Gansu provinces, respectively, are two large cascade reservoirs with a total storage capability of 30.4 billion m 3 [29].

Data Collection
Three Release-05 gravity-field solutions were obtained separately from three different data processing centers: the University of Texas Center for Space Research (CSR; http://www.csr.utexas.edu/grace/), the German Research Center for Geosciences (GFZ; http://isdc.gfz-potsdam.de/grace),and the Jet Propulsion Laboratory (JPL; http://grace.jpl.nasa.gov/).The differences between the solutions obtained from the JPL, CSR, and GFZ were used to infer the uncertainty in the Level-2 GRACE fields [30].Three sets of monthly gravity-field time series of NWC were collected from April 2002 to March 2016, a total of 168 months (13 months with missing data).TWS includes surface water, groundwater, snow cover, and soil moisture.Thus, to evaluate the products is challenging because of the unavailability of a reference dataset (field measurements) for each component of TWS in the study region.Nevertheless, numerous studies demonstrated the reliability of these datasets for monitoring the variation in TWS worldwide on monthly, seasonal, and annual bases [21,22,31,32].
Monthly P, SD, AT, and WS data (from January 1982 to March 2016) were collected from 192 meteorological stations (Figure 1) distributed across NWC, which could be downloaded from the National Meteorological Information Center of China (http://data.cma.cn/).Remote sensing data (MOD13C2 product) based on the moderate resolution imaging spectroradiometer (MODIS) sensors aboard the Terra satellite were downloaded from the Land Processes Distributed Active Archive Center (LP DAAC) of the National Aeronautics and Space Administration (https://lpdaac.usgs.gov/).Data on the monthly variation in storage of the LYX and LJX reservoirs (from April 2002 to March 2016) were obtained from the water year books and water-resources bulletins for Yellow River, which were released by the Yellow Conservancy Commission of the Ministry of Water Resources of China (http://www.yellowriver.gov.cn/).In addition, the annual groundwater withdrawal data (from 2002 to 2015) were collected from the water-resources bulletins released by the respective departments (http://www.sxmwr.gov.cn/;http://www.nxsl.gov.cn/;http://www.gssl.gov.cn/;http://www.qhsl.gov.cn/; http://www.xjslt.gov.cn/) of the five provinces in NWC.

Data Processing
(1) TWS Data Processing In this study, the average monthly gravity field (between January 2004 and December 2009) was taken as the baseline of the monthly gravity-field time series (from April 2002 to March 2016), i.e., a monthly gravity field anomaly was generated following the subtraction of the baseline from the monthly gravity field.The missing monthly GRACE data between April 2002 and March 2016 were interpolated using the linear interpolation method.To eliminate the influence of noise, a Gauss smoothing kernel [33] was introduced into the calculation of gravity-field anomalies, which was expressed as an equivalent water height by the following equation [34]: where ∆h is the equivalent water height; θ is the colatitude; λ is the longitude; a is the equatorial radius (6378 km); ρ ave is the mean density of the earth (5517 kg/m 3 ); ρ water is the density of water (1000 kg/m 3 ); n is the degree of decomposition; m is the order; k n is the loading love number of the nth degree; W n denotes Gauss smoothing kernel related to the nth degree, calculated by the recurrence formula )); r is the filter radius; ∆C nm and ∆S nm are the gravity spherical harmonic coefficient and normalized Stokes coefficient residuals relative to the baseline, respectively; and P nm (cos(θ)) is the nth degree and the mth-order fully normalized Legendre function.
In Equation ( 1), the maximum value of the degree (n max ) and the filter radius (r) were two significant parameters, which were neither too large nor too small.If n max was greater than 60, a larger error would result in the GRACE data [35].If n max was too small, the resolution of the gravity field would decrease.A large r could reduce the noise of the GRACE data; however, it probably resulted in some additional errors when it exceeded a study area [36].In this study, n max = 60 and r = 300 km.In addition, the C 20 (degree 2, order 0) obtained from GRACE was replaced by the C 20 from satellite laser ranging (SLR), because the GRACE-C 20 values had a larger uncertainty than those of the SLR-C 20 values [37,38].
Monthly variations in land gravity fields are mainly caused by monthly TWS changes [39,40].Hence, a monthly equivalent water height represents a monthly anomaly in TWS, herein named TWSA.Based on the GRACE data from three centers (CSR, GFZ, and JPL), three sets of monthly TWSA time series (from April 2002 to March 2016) for NWC were available to support this study.Because the accuracy differences between three sets of time series were not clear, the cross-correlation coefficients between them were adopted to calculate the weights of the data from the three centers as follows: (2) where r(X, Y) is the Pearson correlation coefficient between two time series, X and Y [41]; X k and Y k denote the kth sample values of X and Y, respectively; X a and Y a represent the averages of X and Y, respectively; n is the length of the two time series; w i is the data weight of the ith center (CSR, GFZ, and JPL were separately numbered 1, 2, and 3); r ij denotes the correlation coefficient between two monthly TWSA time series from the ith and jth centers; TWSA is the weighted TWSA time series; and TWSA i is the monthly TWSA time series from the ith center.
When obtaining monthly TWSA time series, the monthly TWS variation could be calculated as shown below [42].
where ∆TWS(i) is the TWS variation of the ith month, and TWSA(i + 1) is the TWSA of the (i + 1)th month.
(2) Climatic Data Processing For a province in NWC, its monthly climatic time series (from January 1982 to March 2016) was calculated using the Thiessen polygon method as follows [43]: where M is the climatic factor (P, SD, AT, and WS) value of a province, l is the total number of meteorological stations in the province, f i denotes the Thiessen polygon area controlled by the ith station, and F is the total area of the province.Specifically, the procedures of the Thiessen polygon method are described below.
(i) Connect the meteorological stations in NWC with straight lines, creating a Delaunay triangulation.
(ii) Determine three perpendicular bisectors and a circumcenter for each triangle in the Delaunay triangulation, thus forming many polygons which take perpendicular bisectors and/or the outline of NWC as boundaries.(iii) Each polygon is controlled by one meteorological station, where the measured climatic factor values represent those over the whole polygon.(iv) Calculate the climatic factor values over a province according to Equation (6).
The Thiessen polygon method assumes that the climatic factors change linearly between two adjacent meteorological stations.The method is widely used in areas where meteorological stations are not uniformly distributed [43].
(3) NDVI Data Processing The monthly remote sensing data (from January 1982 to March 2016) were adopted from a MOD13C2 product with a spatial resolution of 0.05 degrees.The remote sensing data were translated into NDVI data by means of the ENVI software [44].The NDVI data of a province in NWC was the weighted mean value of all gridded NDVI data in the province.Specifically, the weight of each grid was equal to the ratio of overlapping area between the grid and the province to the area of province.

Cross-Wavelet Transformation
The cross-wavelet transformation (CWT) is able to present the correlation between two time series in both time and frequency domains, combined with the wavelet transformation and cross-spectrum analysis [45].The CWT is often used to explore the correlations between two annual hydrological and climatic time series [46].In this study, a CWT based on a Morlet wavelet [47] was mainly employed to investigate resonant periods between the monthly TWSA time series and the monthly climatic and vegetational factor time series.
During the CWT process, W n X (L) and W n Y (L) separately denoted the continuous wavelet transformations of two time series, X and Y, respectively, where n is the length of two time series, and L is the time lag.The CWT between X and Y is expressed as is the complex conjugation of W n Y (L).The power spectrum of the wavelet is defined as |W n XY (L)|, which reflects the correlation degree between X and Y.The Fourier transformation power spectrum (P k ) of red noise was calculated using the following equation [48]: Here, b is the autoregressive coefficient of order 1 for a red-noise time series, and N denotes the length of the red-noise time series.
Supposing that the expectation spectra of the two time series, X and Y, are P k X and P k Y , the distribution of the cross-wavelet power spectrum [49] is expressed as follows: where σ X and σ Y denote the standard deviations of the two time series, X and Y, respectively; Z ν (p) is the confidence level corresponding to the probability p for a probability distribution function by the square root of two χ 2 distributions; and ν is the freedom degree.
In Equation ( 8), if the left-hand value is larger than the upper limit of the power spectrum of red noise under the significance level, α, it is thought that the correlation between X and Y is significant.

Pearson Correlation Coefficient Test
The Pearson correlation coefficient (PCC) denotes a certain dependency relationship between two time series, X and Y, including the direction and degree of correlation.The PCC in Equation ( 2) is a traditional and effective statistic for correlation analysis.Its plus/minus sign and absolute value reflect the direction and degree of correlation between two time series, respectively.If X and Y are stationary, a statistic is built as follows [50]: where t is the statistic following a Student's t-distribution with a degree of freedom, n − 2; r is the PCC between two time series; and n is the length of the two time series.When the significance of r is tested, the null hypothesis is that r is equal to zero, i.e., the two time series have no correlation.At the significance level, α, the null hypothesis is rejected if the absolute value of t is larger than t α/2 .
In reality, a time series may be non-stationary, which could reduce the effectiveness of the PCC.In addition, the periodic components in both time series probably lead to a false correlation.Thus, the resonant periodic components and their respective trend components are removed from X and Y before analyzing their correlation.Specifically, it is assumed that X and Y have significant resonant periods, and their separate periodic components are X P and Y P , respectively.In addition, both X-X P and Y-Y P have separate, overall linear trend components, X T and Y T , respectively.The PPC between X-X P -X T and Y-Y P -Y T is r' with the statistic t'.If |t'| is larger than t α/2 , X and Y have a significant correlation depending on r'; otherwise, it is insufficient to determine the correlation between X and Y based on r'.
For brevity, r and r' were named original PCC and filtered PCC, respectively.It is assumed that m factors (X 1 , X 2 , . . ., X m ) jointly influence a factor (Y). Furthermore, r i (i = 1, 2, . . ., m) and r' i are the original and filtered PCCs between the X i and Y time series.According to the (positive or negative) sign of r*•r' i and the value of |r' i |, the effects of different factors (X 1 , X 2 , . . ., X m ) on Y could be ranked.If r*r' j > 0 and r*r' k < 0 (j, k = 1, 2, . . ., m), the effect of X j on Y is larger than that of X k .If r*r' j and r*r' k have the same sign and |r' j | > |r' k |, the effect of X j on Y is larger than that of X k ; otherwise, the effect of X j on Y is smaller than that of X k .

Data Weights of the Three Centers
Figure 2 shows the correlations of monthly TWSA time series measured by the three centers (CSR, GFZ, and JPL).All correlation coefficients were greater than 0.7, indicating that the TWSA data of NWC measured by the three centers were highly consistent.The correlation between CSR and GFZ was the highest, while the correlation between GFZ and JPL was the lowest.In addition, the correlations of the three centers decreased from east to west in NWC, which probably resulted from the increase in gravity-field inversion error related to the increases in regional size and topographic complexity in this direction [51].Based on the correlations, the data weights of the three centers were calculated, and are presented in

TWS Variations in the NWC
Figure 3 shows the intra-annual variations in TWS in the five provinces of NWC.It can be seen that the TWS showed distinct seasonal variations, and TWS variations in Xinjiang province were quite different from those in the other four provinces.In the Shaanxi, Ningxia, Gansu, and Qinghai provinces, TWS increased from April to September, and basically decreased from September to April.In the Xinjiang province, TWS decreased from April to October, and increased from October to April.Yang and Chen [12] found that TWS reached maximum and minimum values in April and October, respectively, in the Central Xinjiang province, which is consistent with the changing characteristics of TWS in the Xinjiang province in this study.
Using a Mann-Kendall test [52], the trends in Figure 4a show that TWS in Qinghai had a significant trend, with an increasing rate of 1.47 mm/a in Figure 4b, while TWS in the other four provinces (Shaanxi, Ningxia, Gansu, and Xinjiang) obviously decreased.Figure 4b     Using a Mann-Kendall test [52], the trends in Figure 4a show that TWS in Qinghai had a significant trend, with an increasing rate of 1.47 mm/a in Figure 4b, while TWS in the other four provinces (Shaanxi, Ningxia, Gansu, and Xinjiang) obviously decreased.Figure 4b

Correlations between TWS and Climatic and Vegetational Factors
CWT was utilized to analyze the correlations between ΔTWS and climatic and vegetational factors in NWC.For brevity, Figure 5 only shows the CWTs between the monthly ΔTWS and P series in the Xinjiang province of NWC.The main information of the cross-wavelet power spectrum, like that in Figure 5a, could be described as follows: (1) the influencing cone of the wavelet, namely the area surrounded by the fine arc, is set to avoid the boundary effect and the false information outside

Correlations between TWS and Climatic and Vegetational Factors
CWT was utilized to analyze the correlations between ∆TWS and climatic and vegetational factors in NWC.For brevity, Figure 5 only shows the CWTs between the monthly ∆TWS and P series in the Xinjiang province of NWC.The main information of the cross-wavelet power spectrum, like that in Figure 5a, could be described as follows: (1) the influencing cone of the wavelet, namely the area surrounded by the fine arc, is set to avoid the boundary effect and the false information outside the cone; (2) the area surrounded by the thick real line denotes the cross-wavelet power spectrum passing the test of standard spectrum of red noise under the significance level, α = 0.05; (3) the numbers on the right of the color bar represent the relative power spectrum values (dimensionless value); (4) darker red corresponding to a bigger spectrum value denotes a more significant resonant period (RP), while darker blue corresponding to a smaller spectrum value indicates a less significant RP; and (5) the arrow represents the phase relationship between factor 1 (such as precipitation) and factor 2 (such as ∆TWS).Specifically, "→" denotes that the variations of factors 2 and 1 are synchronous, "↓" indicates that the variation of factor 2 lags behind that of factor 1 with one-fourth of an RP, "←" implies that the variation of factor 2 lags behind that of factor 1 with half of an RP, and "↑" shows that the variation of factor 2 lags behind that of factor 1 with three-fourths of an RP.Additional information about wavelet cross spectra was given by Herrera et al. [56].
Without eliminating periodic and trend components, as shown in Figure 5a, monthly ∆TWS and P had a significant resonant period of 12 months.The ∆TWS value was negatively correlated with P as a whole, but monthly ∆TWS and P had a positive relationship according to the terrestrial water budget equation.In the Xinjiang province, the multi-year average ratio of P to potential Et was less than 0.1 [57]; thus, the negative correlation between ∆TWS and Et was far stronger than the positive correlation between ∆TWS and P. As shown in Figure 5b, the negative correlation between ∆TWS and P was greatly weakened, while the positive correlation was enhanced after removing their periodic and trend components.In NWC, monthly ∆TWS and climatic and vegetational factors had a significant resonant period of 12 months, based on the CWT analysis.To obtain reliable correlations, it was necessary to reduce the periodic and trend components before analyzing a correlation between the two factors on a monthly basis.
Overall, correlations between the climatic and vegetational factors were noteworthy, as shown in Table 2.In NWC, the negative correlations were significant between monthly P and monthly SD, AT, and NDVI.Monthly SD had significant positive correlations with monthly AT and NDVI.Monthly AT showed significant positive correlations with monthly NDVI.Monthly WS had no stable correlations with the other four factors.In the Shaanxi, Ningxia, Gansu, and Qinghai provinces, monthly WS had negative correlations with monthly P, and showed positive correlations with monthly SD, AT, and NDVI.However, monthly WS showed a significant positive correlation with monthly P and negative correlations with monthly SD, AT, and NDVI in the Xinjiang province.The significant correlations showed strong interactions between these factors, which could disturb the independent effects of some factors on ∆TWS to a certain extent.
correlation between ΔTWS and P. As shown in Figure 5b, the negative correlation between ΔTWS and P was greatly weakened, while the positive correlation was enhanced after removing their periodic and trend components.In NWC, monthly ΔTWS and climatic and vegetational factors had a significant resonant period of 12 months, based on the CWT analysis.To obtain reliable correlations, it was necessary to reduce the periodic and trend components before analyzing a correlation between the two factors on a monthly basis.Vegetational coverage was influenced by multiple climatic factors, particularly P, SD, and AT [58,59].On a monthly scale, the increase of P led to the decreases of SD and AT, which hindered the normal growth process of vegetation, although soil moisture was replenished by P. Thus, the negative correlations between monthly P and NDVI were reasonable in the five northwest provinces.Nevertheless, the vegetational coverage conditions generally had a good positive correlation with antecedent P in an area; for example, Gessner et al. [60] found vegetational coverage was sensitive to P anomalies of the previous one to three months in Central Asia.
Wind speed has multiple effects on the transport of water vapor [61,62].On one hand, wind can transport water vapor away from a region, improving regional Et.On the other hand, wind can also transport abundant water vapor to a region, promoting regional P. Thus, WS usually had complicated relationships with other climatic and vegetational variables because of its multiple effects on the transport of water vapor.In particular, Xinjiang is far from oceans, and its P is mainly supplied by water vapor from the North Atlantic Ocean and the Indian Ocean [63,64].Thus, a higher WS Note: r' is the filtered PPC between the two factor time series, achieved by removing the periodic and trend components; bold data are significant values under the significance level, α = 0.05.
Strong effects of one factor were also able to disturb the influences of the weaker factors on ∆TWS in NWC.In the Shaanxi, Ningxia, Gansu, and Qinghai provinces, the effects of P completely overshadowed the influences of SD, AT, WS, and NDVI.Conversely, the effect of P was thoroughly overshadowed in the Xinjiang province, where the annual average P (144 mm) was far less than that of the other four provinces Shaanxi (634 mm), Ningxia (275 mm), Gansu (295 mm), and Qinghai (356 mm) from 1982 to 2015.Deng and Chen [3] found that monthly TWS and P had a significant positive correlation in Southern Xinjiang, and showed a negative correlation in Northern Xinjiang.It was noted that the P in Northern Xinjiang obviously exceeded that in Southern Xinjiang [65].In addition, Deng and Chen [3] found that monthly TWS and AT had significant positive correlations in most parts of Xinjiang, which was probably a result of disturbances in their periodic and trend components.

Effects of Climate and Vegetation Changes on TWS Variations
As shown in Figure 6, climatic and vegetational factors had different trends in the five provinces of NWC from 1982 to 2015.Annual AT and NDVI had significant increasing trends in the Shaanxi, Ningxia, and Gansu provinces.Wang et al. [66] also found that vegetation cover had a notable improvement in NWC from 1981 to 2013.Annual SD and WS obviously decreased, while annual AT notably increased in the Qinghai province.Annual P and AT showed significant increasing tendencies in the Xinjiang province, and Deng and Chen [3] detected the same trends of AT and P in the Xinjiang province during recent decades.
Table 4 shows the effects of climatic (P, SD, AT, and WS) and vegetational (NDVI) factors on TWS in the five provinces of NWC.In the Shaanxi, Ningxia, and Gansu provinces, increasing AT and NDVI jointly resulted in decreases in TWS, because the other factors had no significant trends.In the Qinghai province, the decreases in SD and WS jointly caused an increase in TWS, because the increasing AT had a weaker negative influence on TWS than SD and WS.In the Xinjiang province, the increase in AT resulted in a decrease in TWS, because the negative effect of increasing AT was stronger than the positive effect of increasing P.
of NWC from 1982 to 2015.Annual AT and NDVI had significant increasing trends in the Shaanxi, Ningxia, and Gansu provinces.Wang et al. [66] also found that vegetation cover had a notable improvement in NWC from 1981 to 2013.Annual SD and WS obviously decreased, while annual AT notably increased in the Qinghai province.Annual P and AT showed significant increasing tendencies in the Xinjiang province, and Deng and Chen [3] detected the same trends of AT and P in the Xinjiang province during recent decades.3; and *, ↑, and ↓ separately denote non-significant, significantly positive, and significantly negative trends, respectively, based on the results in Figure 6.
The significance trends of climatic and vegetational factors changed the hydrological pattern in NWC.In particular, potential evapotranspiration significantly decreased in the western area of NWC because of the decrease in annual WS, while it increased in the eastern area of NWC because of the increase in AT [67].The increasing AT accelerated the melting of glaciers and snow cover in Xinjiang, leading to a decrease in TWS in the Tianshan Mountains and an increase in TWS in Southern Xinjiang [3,68].Long-term vegetation restoration reduced soil water storage on the Loess Plateau [69,70].In addition, vegetation transpiration apparently increased groundwater evapotranspiration in the arid inland river basins of NWC [71].

Connections between TWS Variations and Socioeconomic Water Consumption
Notably, socioeconomic water consumption (surface water and groundwater consumption) always has direct links with changes in TWS [72].Based on abundant data on well levels and GRACE-based TWS data, many scholars investigated spatiotemporal variations of regional groundwater [73][74][75][76].Nevertheless, it is beyond the scope of this study to examine the linkage between the temporal evolution of groundwater storage and groundwater withdrawal, given the challenge of deriving groundwater storage constrained by in situ observations in the study region.This study investigated the relationships between monthly TWS variations and reservoir operation and groundwater withdrawal in NWC.
As shown in Table 5, ∆TWS in the Shaanxi, Ningxia, Gansu, and Qinghai provinces had significant positive correlations with monthly water storage variation in the cascade reservoirs (LYX and LJX) of the upper Yellow River.Thus, the water storage variation of cascade reservoirs could influence ∆TWS in the four provinces.In addition, the correlation coefficients gradually decreased from east to west in NWC.In Figure 7, in addition to Shaanxi, annual water withdrawals from the Yellow River also decreased in this geographical direction.A reasonable explanation was that the downstream water demands influenced the water storage variation of both cascade reservoirs, and further influenced ∆TWS in the three provinces.Though the water withdrawal in Shaanxi was not the largest, the water withdrawals of Gansu and Ningxia probably supply Shaanxi in the form of groundwater, indirectly causing a good correlation between water storage variation in the two cascade reservoirs and ∆TWS in Shaanxi [77].As shown in Figure 8, the annual groundwater withdrawal of Xinjiang province had an obvious increasing trend, while the annual groundwater withdrawals of the other four northwest provinces of NWC had no noticeable trends.Groundwater exploitations are often conducted in the Shaanxi, Ningxia, Gansu, and Xinjiang provinces because of large socioeconomic water demands in these regions [78,79].As shown in Figure 8, the ranking of the annual average groundwater withdrawals from largest to smallest was Shaanxi (16.3 mm), Ningxia (8.2 mm), Gansu (6.0 mm), Xinjiang (5.0 mm), and Qinghai (0.8 mm), which was consistent with the ranking of decreasing TWS shown in Figure 4b.However, this is not a mere coincidence, as considerable amounts of groundwater are converted into surface water during groundwater withdrawal activities, thereby promoting actual Et and reducing groundwater storage in a province, even under the same climatic and vegetational conditions [80].As shown in Figure 8, the annual groundwater withdrawal of Xinjiang province had an obvious increasing trend, while the annual groundwater withdrawals of the other four northwest provinces of NWC had no noticeable trends.Groundwater exploitations are often conducted in the Shaanxi, Ningxia, Gansu, and Xinjiang provinces because of large socioeconomic water demands in these regions [78,79].As shown in Figure 8, the ranking of the annual average groundwater withdrawals from largest to smallest was Shaanxi (16.3 mm), Ningxia (8.2 mm), Gansu (6.0 mm), Xinjiang (5.0 mm), and Qinghai (0.8 mm), which was consistent with the ranking of decreasing TWS shown in Figure 4b.However, this is not a mere coincidence, as considerable amounts of groundwater are converted into surface water during groundwater withdrawal activities, thereby promoting actual Et and reducing groundwater storage in a province, even under the same climatic and vegetational conditions [80].
In NWC, agricultural irrigation remains the largest source of water usage, accounting for over 70% of total water consumption [79,81].Irrigation consumes substantial amounts of groundwater every year, leading to severe decreases in groundwater tables in many irrigated farming areas worldwide [82,83].The groundwater depth gradually increased with the growing number of irrigation wells on the Guanzhong Plain of the Shaanxi province from 1977 to 2010 [84].The groundwater table significantly declined in the Dunhuang oasis of the Gansu province from 1987 to 2007 because of increased pumping for irrigation [85].The groundwater levels obviously decreased in the Tarim River basin of the Xinjiang province with the expansion of irrigation between 2004 and 2010 [3].With increasing AT, irrigation probably consumed more groundwater in some irrigated areas, thus resulting in a further reduction of TWS in NWC.
from largest to smallest was Shaanxi (16.3 mm), Ningxia (8.2 mm), Gansu (6.0 mm), Xinjiang (5.0 mm), and Qinghai (0.8 mm), which was consistent with the ranking of decreasing TWS shown in Figure 4b.However, this is not a mere coincidence, as considerable amounts of groundwater are converted into surface water during groundwater withdrawal activities, thereby promoting actual Et and reducing groundwater storage in a province, even under the same climatic and vegetational conditions [80].In NWC, agricultural irrigation remains the largest source of water usage, accounting for over 70% of total water consumption [79,81].Irrigation consumes substantial amounts of groundwater every year, leading to severe decreases in groundwater tables in many irrigated farming areas worldwide [82,83].The groundwater depth gradually increased with the growing number of irrigation wells on the Guanzhong Plain of the Shaanxi province from 1977 to 2010 [84].The groundwater table significantly declined in the Dunhuang oasis of the Gansu province from 1987 to 2007 because of increased pumping for irrigation [85].The groundwater levels obviously decreased

Conclusions
This study investigated GRACE-based TWS variations in NWC between April 2002 and March 2016, including seasonal and trend variations.Furthermore, the relationships between TWS variations and climatic and vegetational factors, as well as socioeconomic water consumption, were analyzed.The major conclusions are as follows: (1) TWS showed distinct seasonal variations and a significant decreasing tendency in NWC as a whole.In particular, TWS obviously decreased in the Shaanxi, Ningxia, Gansu, and Xinjiang provinces, while TWS notably increased in the Qinghai province.Increases in AT and NDVI were the main causes of the decreases in TWS in the Shaanxi, Ningxia, and Gansu provinces.
The decreases in SD and WS resulted in an increase in TWS in the Qinghai province, while the decrease in TWS was caused by the obvious increase in AT in the Xinjiang province.(2) The interactions of climatic and vegetational factors were significant, and strong effects of some factors could weaken the influences of other factors on TWS variations in NWC.In particular, the negative effects of SD, AT, and NDVI jointly masked the positive effects of P and WS on TWS in the Xinjiang province, whereas the positive effect of P masked the negative effects of the other factors on TWS in the other provinces in NWC.Accordingly, we should emphasize the analysis of the interactions and combined effects of multiple effects on TWS variations in a region.(3) TWS in the Shaanxi, Ningxia, Gansu, and Qinghai provinces had good correlations with the variation in water storage in the cascade reservoirs of the upstream of the Yellow River, and the correlation coefficients gradually decreased from east to west in NWC.In addition, increasing AT could promote actual Et when more groundwater is converted into surface water in irrigated areas, thus resulting in a further reduction in TWS in NWC.

Figure 1 .
Figure 1.Location of Northwest China (NWC), including meteorological stations and cascade reservoirs in the upstream of the Yellow River.

Figure 1 .
Figure 1.Location of Northwest China (NWC), including meteorological stations and cascade reservoirs in the upstream of the Yellow River.

Figure 2 .
Figure 2. Correlation coefficients between monthly time series of terrestrial water storage anomalies (TWSA) of NWC measured by three data centers from April 2002 to March 2016.CSR-the University of Texas Center for Space Research; GFZ-the German Research Center for Geosciences; JPL-the Jet Propulsion Laboratory.

Figure 3 .
Figure 3. Multi-year averages of monthly TWS in the five provinces in NWC between April 2002 and March 2016.

Figure 2 .
Figure 2. Correlation coefficients between monthly time series of terrestrial water storage anomalies (TWSA) of NWC measured by three data centers from April 2002 to March 2016.CSR-the University of Texas Center for Space Research; GFZ-the German Research Center for Geosciences; JPL-the Jet Propulsion Laboratory.
Figure3shows the intra-annual variations in TWS in the five provinces of NWC.It can be seen that the TWS showed distinct seasonal variations, and TWS variations in Xinjiang province were quite different from those in the other four provinces.In the Shaanxi, Ningxia, Gansu, and Qinghai provinces, TWS increased from April to September, and basically decreased from September to April.In the Xinjiang province, TWS decreased from April to October, and increased from October to April.Yang and Chen[12] found that TWS reached maximum and minimum values in April and October, respectively, in the Central Xinjiang province, which is consistent with the changing characteristics of TWS in the Xinjiang province in this study.Using a Mann-Kendall test[52], the trends in Figure4ashow that TWS in Qinghai had a significant trend, with an increasing rate of 1.47 mm/a in Figure4b, while TWS in the other four provinces (Shaanxi, Ningxia, Gansu, and Xinjiang) obviously decreased.Figure4bshows the linear rates of TWS in the five provinces.Specifically, ranking the rates of TWS decrease from fastest to slowest yields Shaanxi (−2.88 mm/a), Ningxia (−2.78 mm/a), Gansu (−1.7 mm/a), and Xinjiang (−1.24 mm/a).According to the available literature[3,12,[53][54][55], TWS had a significant decline in the Xinjiang province, particularly in the Tianshan Mountains and surrounding areas.TWS had an obvious increase at the intersection between Tibet and the Xinjiang and Qinghai provinces.TWS notably decreased in the Shaanxi and Ningxia provinces between 2003 and 2013.In addition, TWS also obviously decreased in the North Gansu province from 2003 to 2012.

Figure 2 .
Figure 2. Correlation coefficients between monthly time series of terrestrial water storage anomalies (TWSA) of NWC measured by three data centers from April 2002 to March 2016.CSR-the University of Texas Center for Space Research; GFZ-the German Research Center for Geosciences; JPL-the Jet Propulsion Laboratory.

Figure 3 .
Figure 3. Multi-year averages of monthly TWS in the five provinces in NWC between April 2002 and March 2016.

Figure 3 .
Figure 3. Multi-year averages of monthly TWS in the five provinces in NWC between April 2002 and March 2016.
shows the linear rates of TWS in the five provinces.Specifically, ranking the rates of TWS decrease from fastest to slowest yields Shaanxi (−2.88 mm/a), Ningxia (−2.78 mm/a), Gansu (−1.7 mm/a), and Xinjiang (−1.24 mm/a).According to the available literature [3,12,53-55], TWS had a significant decline in the Xinjiang province, particularly in the Tianshan Mountains and surrounding areas.TWS had an obvious increase at the intersection between Tibet and the Xinjiang and Qinghai provinces.TWS notably decreased in the Shaanxi and Ningxia provinces between 2003 and 2013.In addition, TWS also obviously decreased in the North Gansu province from 2003 to 2012.

Figure 4 .
Figure 4. Trends of monthly TWSA time series in NWC from April 2002 to March 2016: (a) significance test results of the overall trends; Z is the statistic of the trend of a time series, and Z_ub and Z_lb separately denote the upper and lower bounds of non-significant trends under the significance level, α = 0.05; (b) linear trend rates of TWS in the five provinces of NWC.

Figure 4 .
Figure 4. Trends of monthly TWSA time series in NWC from April 2002 to March 2016: (a) significance test results of the overall trends; Z is the statistic of the trend of a time series, and Z_ub and Z_lb separately denote the upper and lower bounds of non-significant trends under the significance level, α = 0.05; (b) linear trend rates of TWS in the five provinces of NWC.

Figure 5 .
Figure 5. Cross-wavelet transformations (CWTs) between time series of the variations in TWS (∆TWS) and precipitation (P) in Xinjiang from May 2002 to February 2016: (a) Original ∆TWS and P time series; (b) ∆TWS and P time series with periodic components and linear trend components removed.

Figure 6 .
Figure 6.Significance tests of the trends of annual climatic (precipitation (P), sunshine duration (SD), air temperature (AT), and wind speed (WS)) and vegetational (normalized difference vegetation index (NDVI)) factor time series in NWC from 1982 to 2015.(Z is the statistic of the trend of a time series, and Z_ub and Z_lb separately denote the upper and lower bounds of non-significant trends under the significance level, α = 0.05.).

Figure 7 .
Figure 7. Annual water withdrawals of the four provinces from the Yellow River between 2002 and 2015.

Figure 7 .
Figure 7. Annual water withdrawals of the four provinces from the Yellow River between 2002 and 2015.

21 Figure 8 .
Figure 8. Annual groundwater withdrawals of the five provinces in NWC.

Figure 8 .
Figure 8. Annual groundwater withdrawals of the five provinces in NWC.

Table 1 ,
thus obtaining the monthly weighted TWSA time series of NWC.
Note: CSR denotes the University of Texas Center for Space Research; GFZ denotes the German Research Center for Geosciences; and JPL denotes the Jet Propulsion Laboratory.

Table 1 .
Weights of monthly the Gravity Recovery and Climate Experiment (GRACE) data released by the three centers from April 2002 to March 2016.
Note: CSR denotes the University of Texas Center for Space Research; GFZ denotes the German Research Center for Geosciences; and JPL denotes the Jet Propulsion Laboratory.

Table 3 .
PCCs between monthly variations in terrestrial water storage (∆TWS) and monthly climatic (P, SD, AT, and WS) and monthly vegetational (NDVI) factor time series in NWC from May 2002 to February 2016.

Table 4 .
Analysis of climatic (P, SD, AT, and WS) and vegetational (NDVI) factor effects on TWS in the five provinces in NWC.: + and − denote positive and negative effects on TWS, respectively; numbers represent the ranking of effects according to r and r' in Table Note

Table 5 .
Correlations between monthly ∆TWS in the four northwest provinces and the monthly water storage variation of cascade reservoirs at the upper Yellow river between May 2002 and February 2016.
Note: bold data are significant values under the significance level, α = 0.05.2016.Note: bold data are significant values under the significance level, α = 0.05.