Regional Characteristics and Impact Factors of Change in Terrestrial Water Storage in Northwestern China From 2002 to 2020

This article characterized the linear trends and interannual signals of terrestrial water storage (TWS) and meteorological variables including precipitation (P) and evapotranspiration (ET) over the arid Northwestern China (NWC). The relative impaction of P, ET, and human water utilization (HU) on TWS variation among the 10 watersheds of NWC (i.e., 5 watersheds in Xinjiang, 3 watersheds in Hexi Corridor, and 2 watersheds in Qinghai) were then investigated. The result indicated that groundwater storage (GWS) was the main contributor to the TWS variation and matched well with TWS in spatial features or watershed-scale variations. The entire NWC presented growth trends for P (0.05 cm/year) and ET (0.22 cm/year) and decline trends for TWS (−0.19 cm/year) and GWS (−0.20 cm/year). The watersheds in Qinghai province where mainly affected by natural factors showed the increasing TWS/GWS trend. The watersheds in Xinjiang and Hexi Corridor, which had strong impact from human activities generally showed the declining TWS/GWS trends, but Xinjiang showed more intensive declining trend than Hexi Corridor. The analysis of HU indicated that water sustainable management and water-saving technologies had effectively kept down the tendency of TWS/GWS declining in the watersheds in Hexi Corridor, however, they were not sufficient to address the water shortage caused by farmland expansion, slight P growth, and high ET growth in Xinjiang. Groundwater use, as the main source to compensate for the increase in HU (especially agricultural water use), exacerbated TWS/GWS loss in Xinjiang. This article provides valuable information for the water management over the arid NWC.


I. INTRODUCTION
W ATER scarcity is the main constraint to socio-economic development in arid and semiarid regions. Global warming has further heightened hydrological fluctuations and water resource uncertainties and ultimately controls the amount of water available for human and ecological consumption [1]. Revealing the regional change features and mechanisms of terrestrial water storage (TWS) is critical, because TWS is a significant and dynamic constituent of the hydrologic cycle and includes all forms of water stored above and below the land surface, such as surface water, soil moisture, groundwater, snow, ice, and biomass water. The in situ observation networks are not enough to directly monitor the watershed-scale TWS change features. Therefore, modeling approaches based on physical processes, including land surface models (LSMs) [2], [3] and global hydrologic models [4], [5], are applied to identify macroscale TWS changes. Though the models are based on rigorous physical principles, many factors will affect the model outputs and consequently misestimate the climate and human-induced changes in the TWS [6]. Fortunately, the Gravity Recovery and Climate Experiment (GRACE) and its follow-on (GRACE-FO) missions measure the monthly variations of gravity fields to track TWS changes at a global or regional scale [7]. Currently, GRACE-TWS itself or its combination with model data has been widely applied to retrieve the variation of TWS and its components (e.g., groundwater, soil moisture, or snow-mass) [8], [9], [10], [11], [12], evaluate the water budget [13], [14], [15], define drought [16], [17], and flood risks [18].
Arid Northwestern China (NWC), located in the hinterland of the Eurasian continent and the rain shadow of the Tibetan Plateau, has a continental and dry climate with an annual mean precipitation (P) of less than 200 mm [19] and an annual mean pan evaporation of generally higher than 2000 mm [20]. However, the NWC have been one of important agricultural production regions for China. By virtue of GRACE and LSM datasets, previous researches revealed TWS decline in some watersheds of NWC and, thus, imply the stress to ensure local water security, such as Tarim Basin [21], Hexi Corridor [12], Junggar Basin [22], [23], and Qaidam Basin [24]. However, comparison among watersheds is rare under the context of spatiotemporal heterogeneities in climatic factors and human activities in NWC.
The NWC region has been getting warmer and wetter in recent decades [19], [25] in the context of global warming and present spatiotemporal heterogeneity due to land-atmosphere interaction and orographic effects [26]. In the twenty-first century, the warming and wetting trends in NWC have become more pronounced [27]. It is noteworthy that the decrease in TWS occurs in the context of increasing P input. The P distribution over NWC present temporal and spatial differences [28], [29]. Increased external moisture and internal intensifying recycling are the primary and secondary causes of the P growth [30], [31]. In the mountain-oasis-desert direction, the P and pan evaporation presents a decreasing and increasing trend, respectively [32]. The P and evapotranspiration (ET) show their highest values in the summer over the region. The P growth mainly occurs in winter and summer [26], while the greatest increase in temperature occurs in winter [33]. Taking Xinjiang as an example, winter P increased significantly in Northern Xinjiang, while summer precipitation increased significantly in Southern Xinjiang [34].
The NWC region is an important food-and cotton-producing area, where agriculture mainly relies on irrigation. The rapid expansion of crop acreage along with industrial and population growth has raised concerns about the shortage of the water supply. According to China Water Resources Bulletin 2020, agricultural water consumption accounts for 87%, 85%, and 77% of the water supplies in Xinjiang Uygur Autonomous Region (Xinjiang), Qinghai province, and Gansu province, respectively, which include the NWC region. Additionally, the agricultural structures and farmland ratios are spatially different in the NWC region (e.g., cotton, wheat, and corn in Xinjiang, wheat and corn in the Hexi corridor of Gansu province, and an oil plant in Qinghai province with the least amount of cropland) [35].
Growing water demand, coupled with water scarcity, poses a major threat to arid and semiarid NWC. Groundwater, as one of the main freshwater resources, has been experiencing overabstraction in the region [12]. The ecosystem is also under threat of degradation because of its vulnerabilities to changing hydroclimatic conditions. The regional differences in climatic factors (i.e., P and ET) and human activities can cause the spatiotemporal heterogeneity of the TWS change, necessitating a thorough and systematic comparative analysis in watershed scale. This article seeks to examine the spatiotemporal characterization of hydrological (i.e., TWS and its components) and meteorological variables (P and ET) and human water utilization (HU) over NWC by combining GRACE/GRACE-FO, global land data assimilation system (GLDAS), IMERG, global land evaporation Amsterdam model (GLEAM), and statistical data of water resource bulletins. The seasonal trend decomposition using loess (STL) method, correlation analysis and water balance analysis are used to evaluate the variation features and relationships of all variables during 2002 and 2020. The study's findings will provide new perspectives and a better understanding of TWS changes across NWC, as well as guidelines for water source sustainability.

A. Study Region
The NWC region was subdivided into 10 enclosed watersheds (i.e., endorheic basins without outflow of rivers) in consideration of the spatiotemporal features of hydrological and meteorological variables and human activities (see Fig. 1) [36]. Jungar Basin (JG), the Manas River Basin (MN), and the Ebinur Lake Basin (EB) are located in North Xinjiang. The Tarim River basin (TR) and Turpan-Hami Basin (TH) are located in South Xinjiang. Watersheds in North and South Xinjiang are bounded by the Tianshan Mountains. The Hexi Corridor to the east of the Qilian Mountains contains the Shule River Basin (SL), Shiyang River Basin (SY), and Heihe River Basin (HH). The Hexi Corridor was mainly in the jurisdiction of Gansu province. The Qaidam River Basin (QB) and Qinghai Lake Basin (QH) are in the Eastern Tibetan Plateau (east of the Eastern Kunlun Mountains) and belong to Qinghai province. Vapors transported into the NWC region by the mid-latitude westerlies, high-latitude northerlies, and southwestern and southeastern monsoons mainly form P over the mountainous areas (i.e., the Tianshan Mountains, Western and Eastern Kunlun Mountains, and the Qilian Mountains; Fig. 1) and their neighboring piedmont areas [28], [29]. The water vapor of the southwestern monsoon passing over the Western Kunlun Mountains is rare. The water vapor from the westerlies is transported into the NWC region through the Tianshan Mountains [26]. The southeastern monsoon provides the water vapor for the P around the Qilian Mountains [37]. Surface water bodies (e.g., rivers and lakes) are generally expanding to accommodate the growth of P and/or glacier meltwater [26], [38], [39], [40], [41].

B. Data
Datasets of P, ET, HU, TWS, and TWS components are listed in Table I. 1) GRACE/GRACE-FO: The three monthly GRACE/ GRACE-FO datasets used in this article are processed using mass concentration (mascon) solutions at the Center for Space Research (CSR, University of Texas, Austin) [42], the Jet Propulsion Laboratory (JPL, NASA, and the California Institute of Technology, California) [43], [44], and the Goddard Space Flight Center (GSFC, NASA) [45]. The GRACE/GRACE-FO to conventional spherical-harmonic solutions, TWSAs from mascon solutions have better signal-to-noise ratios and are easily applied with little empirical post-processing requirements (e.g., destriping or spatial smoothing) [46]. The three datasets have been corrected with a glacial isostatic adjustment model and are released at a spatial resolution of a 0.25°× 0.25°, 0.5°× 0.5°a nd 1°× 1°grid-scale, respectively. To facilitate the easy calculation and plotting of the water storage changes, they are resampled to a 0.25°× 0.25°grid scale by nearest-neighbor interpolation.
2) Global Land Data Assimilation System: Soil moisture storage (SMS), snow water equivalent (SWE), surface water storage (SW), and canopy water storage (CW), as the components of TWS, were obtained from the LSMs of the GLDAS. The SW and CW were ignored for their lower values than other components of TWS in the arid NWC region. GLDAS LSMs are built on a detailed energy budget and have the same forcing datasets for each version. GLDAS LSMs have been commonly used as a substitution for in situ measurements due to the scarcity of the latter in the NWC region [12], [47]. GLDAS version 2.1, including NOAH, VIC, and CLSM, were chosen in this article. The NOAH is a 0.25-degree resolution. VIC and CLSM of 1-degree resolution were also resampled to a 0.25-degree resolution via nearest-neighbor interpolation. To be consistent with TWSA, the SMS and SWE values of 3 LSMs were converted into anomaly values (SMSA and SWEA) by subtracting the respective baseline averages during 2004-2009. More detailed information about the GLDAS-2.1 datasets can be found at https://ldas.gsfc.nasa.gov/gldas/.
3) Meteorological Data: Integrated multisatellite retrievals for the global precipitation measurement mission (IMERG) version 6 is the latest release that provides P estimates on a 0.1°× 0.1°grid beginning in 2000 [48]. Compared with other P products, IMERG is considered suitable for this study because it can capture micro and solid P more accurately [49], [50].
The P records from 39 meteorological stations, maintained by the China Meteorological Administration, were also used in this article. The GLEAM is a set of algorithms that separately estimate ET and its components [51]. GLEAM evapotranspiration products have a 0.25°resolution. Since their development, GLEAM evapotranspiration products have been widely used in hydrometeorology fields [52].

4) Water Resources Utilization Statistics:
Statistical data of water utilization from Xinjiang, Qinghai, and Gansu provinces were collected from statistical yearbooks and water resource bulletins during 2002-2020. The watershed-scale data of water utilization are the sum of statistical data from sites within it. The water utilization data of JG, MN, and EB in North Xinjiang cannot be obtained due to interbasin administrative divisions. Hence, the water utilization data in North Xinjiang was summed instead of JG, MN, and EB. The QB and QH in Qinghai province were missing data for 2002. North Xinjiang was missing data for 2002, 2003, and 2012.

C. Methods
Related methods and data used in this article are shown in Fig. 2.
1) Interpolation of GRACE Data: Due to the impossibility of obtaining large-scale field measurements of TWS to calibrate the GRACE/GRACE-FO datasets, the mean of the solutions sourced from different processing centers is commonly used to reduce the uncertainty of TWSA derived from GRACE/GRACE-FO datasets [53]. In this article, we use the uncertainty-weighted average of the CSR, JPL, and GSFC-sourced TWSAs, which is calculated as [54] [55]. First, TWSA after the removal of linear trends is used to calculate the climatology signals (mean annual cycle); second, TWSA after the removal of climatology signals is used to obtain climatology-free signals (interannual variation) by cubic spline interpolation; and finally, climatology-free signals and climatology signals are summed up to fill the missing values in each mission. One reconstructed dataset of Terrestrial Water Storage (RTWSA) is applied to fill the 11 months of missing values between the two missions [56]. RTWSA spanned from 7/1979 through 9/2019 and was constructed using TWSA from JPL and GPCP 2.3 terrestrial precipitation using the cyclostationary empirical orthogonal functions method. Considering that RTWSA only includes interannual and seasonal components, linear trend, and residual components acquired from the TWSA time-series by linear and cubic spline interpolation are added back to RTWSA to fill the 11-month data gap.
2) Signal Decomposition: Signal decomposition has been a common practice for addressing multiple issues [6], [17]. Taking the GRACE signal as an example, it can be decomposed in terms of water storage compartments by (3), water balance items from hydrological and meteorological variables and HU by (4), and temporal components by (5). Runoff is not considered in (4) since the 10 watersheds in NWC are endorheic basins without outflow rivers where GWSA is groundwater storage anomaly relative to the baseline average during 2004-2009, and ΔTWSA is the TWSA change during a time period. Like previous research, linear trends were used to perform the water balance analysis among hydrological and meteorological variables [15], [24]. The Mann-Kendall (MK) method was used to obtain linear trends. As for the monthly time series data of hydrological and meteorological variables, linear trends were obtained via the MK method after the STL [57]. The STL procedure is a robust decomposition method and has been demonstrated for hydrological and meteorological time series [53], [58]. As shown in (5), the original signal (X total ) is decomposed into long-term (X long ), seasonal (X seas ), and residual components (X sub ). The STL algorithm performs smoothing on the time series in two loops, an innerloop nested inside an outer loop (a minimization of the effect of outliers). X seas is calculated first and removed to calculate the trend component during the inner loop. The long-term component is further decomposed into a linear trend (X lin ) and an interannual variation (X inter ) using the MK method. For more detailed procedures for STL, refer to Cleveland et al. (1990) and Humphrey et al. (2016).

3) Correlation Analysis:
The Pearson correlation coefficient (R) is an extensively used indicator to measure correlation between time series. A temporal lag typically exists among hydrological and meteorological variables [53], [59]. The maximum R among variables is obtained by moving a time-series with a step of one month in the range of 11 months. As shown in (6), the R is the product of the covariance of the two variables divided by their standard deviations where n is the sample size of the variables; xi and yi are the values of two variables at a given moment (month or year);x andȳ are the means of two variables of the entire time series. The significance test of the correlation coefficient is passed when the P-value is less than 0.05 [60].

A. Spatial Features of Hydrological and Meteorological Variables
The spatial patterns of TWSA and GWSA matched very well (see Fig. 3(a)-(d)], i.e., the growth and decline areas of the GWSA showed spatial coincidence with the TWSA. Both TWSA and GWSA presented positive monthly mean and growth trends in the southern NWC (the QB, its neighboring southern part of TR, and western parts of SL and QH), and vice versa in northern NWC. Contrarily, the SMSA mainly presented positive monthly mean and growth trends in northern NWC, and vice versa in southern NWC [see Fig. 3(e) and (f)]. Of note, the monthly mean and linear trends of SMSA tend to be positive in basin areas rather than in the mountain areas (e.g., the Tianshan area refer to Fig. 1) [see Fig. 3(e) and (f)].  P mainly occurred in the northwestern and southeastern NWC, i.e., the Tianshan and Qilian Mountains and the adjacent areas [see Fig. 4(a)]. Generally, linear trends of meteorological stations were consistent with those of the corresponding raster cells in IMERG [see Fig. 4(b)]. Increasing P tended to be in basin areas rather than mountainous areas, especially in Xinjiang and QB [see Fig. 4(b)]. Furthermore, watersheds in Xinjiang had a larger proportion of area showing P decline than other watersheds in Qinghai province and the Hexi Corridor. The monthly mean of ET spatially matched very well with that of P, but the former had a wider extent with a higher value [see

B. Watershed-Scale Variations of Hydrological and Meteorological Variables 1) Long-Term Trends in Hydrological and Meteorological
Variables: TWSAs and GWSAs presented increasing trends in the watersheds of Qinghai province, decreasing trends in the watersheds of Xinjiang, and increasing trends in the watersheds of the Hexi Corridor, except for SL (see Fig. 5 and Table II). The watersheds in North Xinjiang, especially EB (−1.49 cm/year), experienced more intensive TWS/GWS loss. The TWSA and GWSA trends of all watersheds but SL passed the significance test. The GWSA and TWSA trends of all watersheds except QH were roughly equal, whereas the trend values of SMSAs in all watersheds but QH were at least one order of magnitude lower than that of TWSAs and GWSAs (see Table II). Only 3 watersheds showed decreasing SMSAs, that is, QB in Qinghai province, TR in South Xinjiang, and MN in North Xinjiang. The SMSAs in watersheds failing to pass significance tests exhibit a slighter change. The trend values of SWEAs were far less than other hydrological variables (generally an order of magnitude lower than SMSAs). Hence, GWSA was the main contributor to the TWSA variation.
The P and ET presented increasing trends in all watersheds except EB in North Xinjiang, which experienced a decreasing P. The trend values of P were generally lower than the ET in all watersheds but QH in Qinghai province and SY and HH in the Hexi Corridor. The P and ET difference was more obvious in the watersheds of Xinjiang, where the trend values of P were at least an order of magnitude lower than those of ET (the P in EB was even decreasing). Furthermore, the P (ET) trend values in Xinjiang were far lower (higher) than in Qinghai province and the Hexi Corridor. Hence, the P and ET difference supports more serious TWS/GWS loss in Xinjiang rather than in the Hexi Corridor and Qinghai province.
The TWSA variation is the result of the accumulative net P (P minus ET) for endorheic watersheds. Hence, curve patterns of TWSA and accumulative net P anomaly (CNP) generally  Table S1 in Supplementary Tables.   TABLE II  LINEAR Table S1). The linear trends of CNP are also calculated in the same time intervals (see Fig. 5; Table S1) [53], [61]. With the exclusion of linear trends that failed to pass the significance test and had opposite trends, the trend values of CNP and TWSA presented a strong positive correlation with the fitting slope of 1.69 (R 2 = 0.7). Therefore, the TWSA variations respond well to the accumulative effects of P input and ET output and verifies the validity of the IMERG and GLEAM (see Fig. 6).

2) Correlation Analysis Among Meteorological and Hydrological Variables:
The strong correlation between P and ET (R values are greater than 0.6 except the 0.4 in QH) agrees with the dominant summer precipitation in NWC (see Fig. 7). The strong correlations between TWSA and GWSA with all R values greater than 0.6 also imply that GWSA is the major contribution to TWSA variation (see Fig. 7). The 7-11 lag months between TWSA/GWSA and P/ET were greater than that between SMSA and P/ET (less than 5 months, except EB and QH), which can be ascribed to stronger time-lag of GWS responding to P/ET (see Fig. 7).
In contrast to the watersheds in Xinjiang, the correlations between P and TWSA/GWSA are obviously better than those between ET and TWSA/GWSA in the Hexi Corridor and Qinghai province with the exception of QH (see Fig. 7). This feature refers to the ET (P) playing more important roles in the TWSA/GWSA variation in the watersheds of Xinjiang (the Hexi Corridor and Qinghai province). The correlation of GWSA with ET presents higher R values than that of SMSA in Xinjiang (except TR), and vice versa in the Hexi Corridor and Qinghai province (see Fig. 7). This feature implies GWS was used more intensively in Xinjiang, which was finally lost due to evaporation.

C. Multiyear Change of Water Supply Statistics in NWC
On the whole, the groundwater extraction trends of watersheds in Qinghai province and the Hexi Corridor were much lower than that of watersheds in Xinjiang (see Table III). Only HH experienced an obvious increase of groundwater extraction (0.242 cm/year) in the Hexi Corridor. The SL had the maximum reduction of groundwater extraction (−0.737 cm/year). Other watersheds in Qinghai province and Hexi Corridor had slight changes of groundwater extraction (<0.02 cm/year). North Xinjiang and TR in South Xinjiang experienced higher increases of groundwater extraction (0.391 and 0.283 cm/year, respectively). It is noteworthy that watersheds with slight changes in groundwater extraction had decreasing agriculture water use with the exception of QB (see Table III). The increase in agriculture water use in QB (0.024 cm/year) was supplied by surface water (0.032 cm/year). However, watersheds (e.g., TR, North Xinjiang and SY) with intensive changes in groundwater extraction were closely related to intensive changes in agriculture water use with the exception of HH. The decreasing agriculture water use (−0.768 cm/year) can interpret all reductions in groundwater extraction (−0.737 cm/year) in SY. The increasing agriculture water use (0.289 cm/year) can interpret all groundwater extraction growth (0.283 cm/year) in TR. The increasing agriculture water use (0.156 cm/year) can interpret 40% of groundwater extraction growth (0.391 cm/year) in the North Xinjiang. As for HH, the agriculture water use decreased, and groundwater extraction growth (0.214 cm/year) was mainly related to a decline in surface water supply (−0.172 cm/year). In addition, the utilization of water resources in watersheds in Qinghai province was significantly lower than watersheds in the Hexi Corridor and Xinjiang based on the yearly average of groundwater extraction and surface water supply during 2002-2020 (see Table III). Based on the trend difference between groundwater extraction and surface water supply, it showed that groundwater was the main source to compensate for the water use increase (see Table III).

D. Water Balance Analysis Based on Mutisource Datasets
The water balance analysis was based on the linear trends and yearly values of hydrological and meteorological variables  Fig. 8. Correlation between TWSA trends (TWSA_Trend) and trend differences of P with ET (P-ET_Trend) and P with ET and HU (P-ET-HU_Trend). The correlation between GWSA trends (GWSA_Trend) and P-ET-HU_Trend (only considering the groundwater extraction in HU term). The green dots of QB and SY were treated as outliers and only black dots were used for the linear fitting operation. and the statistical data for the water supply during 2002-2020 (see Tables II and III). The variable names with suffixes of "_Trend" and "_YV" represent their linear trends and yearly values, respectively. The approach was widely used in previous research [15], [24]. HU was total water supply. The statistical data in North Xinjiang cannot be divided by watersheds in Water Resources Bulletins. For the convenience of comparison, the yearly trends and values of TWSA, GWSA, P, and ET were calculated in North Xinjiang. The trends were −0.955, −0.980, −0.042, and 0.483 cm/year, respectively. Yearly value series were listed in Table S2.
For endorheic watersheds, TWS changes are closely related to P and ET, which should be the only recharge and loss items in the water balance formula. That finding was also supported by the significantly positive correlation between TWSA_Trend and P-ET_Trend [see R 2 = 0.87, Fig. 8(a)]. However, the slope (1.82) of the best-fit line implies a large value gap between TWSA_Trend and P-ET_Trend (i.e., the P trend minus ET trend). The R 2 (0.95) and fitting slope (1.28) of TWSA_Trend and P-ET-HU_Trend (i.e., the P trend minus ET and HU trends) were higher and lower than that of TWSA_Trend and P-ET_Trend and suggested that the former had a better fitting result and a lower value gap [see Fig. 8(a) and (b)]. Hence, the human-induced TWS loss cannot be ignored in the water balance analysis. It is also inferred that GLEAM dataset may missed a large amount of human-induced ET loss since the water was eventually lost in terms of ET for endorheic watersheds. The SY seriously deviated from the fitting line of the water balance analysis may be ascribed to the long-term overextraction of groundwater with a largest multiyear mean of 22.581 cm, which cannot be easily recovered (TWSA has an accumulative effect). The TWSA started to increase in response to the CNP growth until 2017 in SY [see Fig. 5(d)]. As for QB, the deviation from the fitting line of the water balance analysis may be ascribed to an underestimation of P (specially more frequent extreme P in recent years) in areas with a high elevation [62], [63].
The high R 2 (0.96) and low value gap (slop of best-fit line was 1.1) between GWSA_Trend and P-ET-HU_Trend also testified to the important influence of groundwater extraction on GWSA changes. It should be noted that the QB and SY were treated as outliers in fitting operations among TWSA_Trend, P-ET_Trend, and P-ET-HU_Trend for their excessive deviation from the basic rule; that is, TWS changes should be positively correlated with P-ET or P-ET-HU in endorheic watersheds.
The negative correlation in QB and SY and uncorrelation in QH were also ascribed to the abovementioned reasons [see Fig. 9(a), (b), and (e)]. The P-ET_YV values were far less than HU, then negative correlation mainly reflected the synchronous decrease of TWSA_YV and HU_YV in SY [see Fig. 9(e)]. The fitting line between the yearly values of TWSA_YV and P-ET-HU_YV also indicate that a good water balance analysis can be obtained after considering HU in Xinjiang and HH in the Hexi Corridor [see Fig. 9(g)-(i)]. As for SL in the Hexi Corridor, the TWSA was similar to CNP and P-ET-HU in the curve patterns [see Figs. 5(f) and 9(c)]. Hence, the weak correlation of the yearly TWSA_YV and P-ET-HU_YV values may be related to an underestimation of P in an upstream area such as QB (see Fig. 1).
Water balance analysis was mainly based on four different types of datasets (i.e., TWS of GRACE, P of IMERG, ET of GLEAM and HU from water resource bulletins). Even through the good fitting between TWSA_trend and P-ET-HU_Trend, there are still residue exited among TWSA, P, ET, and HU for the >1 slope (1.28) and the intercept (0.07) in Fig. 8(b). All watersheds are plotted away from y = x line and indicate that every watershed shows residue in water balance analysis. The IMERG may underestimate the P in high altitude areas. The better R 2 value and lower slope of fitting line support that ET loss induced by human is not completely observed by the GLEAM [see Fig. 8(a) and (b)]. Quantitative analysis of the contribution of the four datasets to the residue in water balance analysis will be covered in the future work.

A. Regional Comparison of TWSA Variation
The NWC region has been getting warmer and wetter in recent decades [19], which is widely manifested as P growth [33], [64], increasing ET [31], lake expansion [40], [41], [65], glacier shrinkages [25], [38], runoff increase [25], and so on. Despite the obvious growth of surface water bodies (runoff and lake volume) and P, the TWS has been generally decreasing in Xinjiang and the Hexi Corridor during 2002-2020 (see Table II and Fig. 3). The growth of surface water bodies was simultaneously ascribed to P growth and meltwater growth from glacier/snow. Surface water bodies' growth from glacier/snow meltwater only reflected the spatial water redistribution. In addition, surface water bodies are easily impacted by human activities and only reflect one part of TWS in NWC and its neighboring arid regions [66]. Although lake expands in recent years, the water resource shortage has not changed in NWC. Once disorderly development occurs, the lakes will shrink again, as the Ebinur and Bosten lakes experienced in previous years [65]. Thus, the information from the TWS variation was more important for water resource management.
GWSA was the main contributor to the TWSA variation and matched well with TWSA in the spatial features or watersheds (see Fig. 3 and Table II). The spatial or watershed-scale difference was closely related to groundwater extraction. In the Hexi Corridor and Xinjiang, watersheds experienced stronger groundwater extraction (i.e., high multiyear trend or mean) and presented an obvious TWS decrease during 2002-2020 (except for SL; Figs. 3 and 5; Tables II and III). As for SL in the Hexi Corridor, the watershed-scale TWS growth was ascribed to mountainous areas with rare human activities [see Fig. 3(a) and (b)]. In Qinghai province, watersheds experienced slighter groundwater extraction (i.e., low multiyear trend and mean) and presented a TWS increase during 2002-2020 (see Fig. 3; Tables III and III).
The water balance analysis further testifies to the importance of considering groundwater extraction (see Fig. 7) and confirms the validity of the multisource data in the article. Furthermore, water balance items (i.e., P, ET, and HU) also presented a clear difference in watersheds. Watersheds in Xinjiang generally presented a lower P growth trend (even decline in EB), higher ET growth trend, and higher HU growth trend (except for TH). As a result, a high TWS loss trend was shown in the watersheds of Xinjiang. Taking the TH as an example, the P trend was 0.048 cm/year. And even without an obvious increase in groundwater extraction (0.02 cm/year), the ET trend was already 0.379 cm/year in TH, not to mention the rest of the watersheds in Xinjiang, where groundwater extraction increased significantly. However, the watersheds in Qinghai province and the Hexi Corridor generally presented a higher P growth trend, lower ET growth trend, and lower HU trend (except for HH). As a result, TWS growth was shown in Qinghai province, and a low TWS loss trend was shown in the Hexi Corridor (except for SY). The abovementioned regional difference was also supported by the correlation analysis of interannual signals. P (ET) played a more important role in the TWSA variations of watersheds of the Hexi Corridor and Qinghai province (Xinjiang). Compared with Qinghai province and the Hexi Corridor, ET was more related to GWSA than was SMAS and consistent with a slight P growth (even P decline) and strong groundwater extraction in Xinjiang (see Fig. 6; Table III).

B. Influence of Agricultural Water Use on TWSA Variation
Agricultural water accounted for >70% of the water supply in NWC during 2002-2020. The agricultural water use of QB and QH in Qinghai province was mainly from surface water (see Table III). In addition, flood irrigation using surface water was the main irrigation method and also recharged the groundwater in the watersheds of Qinghai province [67]. Both the Hexi Corridor and Xinjiang are important agricultural production sites. Groundwater had been an important component of the water supply in the Hexi Corridor and Xinjiang to compensate for the water demand of agricultural activities (see Table III). Although water-saving technologies were widely applied and irrigation efficiency in NWC have been greatly improved [35], the core conflict still exists between water demand from cropland expansion and agricultural structure and the insufficient water supply, especially for Xinjiang. In Xinjiang, the cropland expanded from 2.7 × 10 4 in 2002 to 4.9 × 10 4 km 2 in 2020 (from 0.9 × 10 4 to 2.5 × 10 4 km 2 for cottons; Statistical yearbook 2002-2020) and caused an agricultural water increase at a rate of 0.22 cm/year (see Table III). The cottons accounted for >70% of cropland expansion and >67% of total agricultural water [68], [69]. As a result, groundwater extraction (0.24 cm/year) was increased to compensate for agricultural water augmentation (see Table III) and caused a continuous decrease in TWSA and GWSA (see Table II). Besides the good fit of multiyear trends in water balance items (see Fig. 7), the yearly values of the water balance items also presented a good fit for Xinjiang [see Fig. 8(a) and (b)], indicating an effective method for managers to decide groundwater extraction amount based on GLEAM, IMERG, and yearly statistical data.
In the Hexi Corridor, the cropland expanded from 0.9 × 10 4 km 2 in 2002 to 1.2 × 10 4 km 2 in 2020. More than 60% of the crop planting area was occupied by cereal crops [70]. The total planted area of cotton had been decreasing from 2007 and was 166 km 2 in 2020 in the Hexi Corridor. In total, 31% of the irrigation farmland in the Hexi Corridor had implemented water-saving irrigation to deal with water scarcity [71], and the agricultural water use efficiency had been improving [72]. The decrease in agricultural water use was ascribed to the significant efforts in water sustainable management, some of which were already implemented in the Hexi Corridor before 2015 (see Table III) [71], [73]. Taking SY as example (experienced most serious TWS depletion or groundwater extraction), the TWS started to increase since 2017 after the most significant reduction of agricultural water use [see

V. CONCLUSION
This article examined the TWSA/GWSA variation in response to the warming and wetting trend and carried out a water balance analysis to provide information to support local water resource management in the NWC. The main conclusions are as follows.
GWSA was the main contributor to the TWSA variation and matched well with TWSA in spatial features or watershed-scale variations. Despite the differences in values, the watersheds in NWC can be grouped into three subregions (i.e., watersheds in Xinjiang, watersheds in the Hexi Corridor, and watersheds in Qinghai province) by hydrological and meteorological variables and the statistical data of water resource bulletins. Watersheds in Xinjiang generally presented a lower P growth trend (<0.1 cm/year; even decline in EB), higher ET growth trend (generally greater than 0.3 cm/year), and higher HU growth trend (generally greater than 0.2 cm/year except for TH). As a result, a high TWS loss trend was shown in the watersheds of Xinjiang (generally lower than −0.4 cm/year except for TR). However, watersheds in Qinghai province and the Hexi Corridor province generally presented a higher P growth trend (generally greater than 0.1 cm/year except for QB), lower ET growth trend (lower than 0.3 cm/year), and lower HU trend (generally present negative values and lower than 0.1 cm/year). As a result, TWS growth was shown in Qinghai province (greater than 0.1 cm/year) and a low TWS loss trend was shown in the Hexi Corridor (greater than −0.4 cm/year). The abovementioned regional difference was also supported by the correlation analysis of the interannual signals of the hydrological and meteorological variables. According to R values, P (ET) played a more important role in the watersheds of the Hexi Corridor and Qinghai province (Xinjiang). Compared with Qinghai province and the Hexi Corridor, ET was more related to GWSA than SMAS, which was consistent with a slight P growth (even P decline) and strong groundwater extraction in Xinjiang.
Water sustainable management and water-saving technologies were effective in the watersheds of the Hexi Corridor. Even for the SY in the Hexi Corridor (experienced most serious TWS depletion and groundwater extraction), the TWS started to increase in 2017 under the coupled function of P growth and water resource management. However, the watersheds in Xinjiang had been facing serious conflict between water demand from cropland expansion (especially for cotton planting area) and an insufficient water supply. In addition, P growth trends were at least an order of magnitude lower than ET growth trends in Xinjiang. To address water scarcity in Xinjiang, more stringent water management policies are needed to reduce agricultural water use, including improvement of water use efficiency, control of the farmland area, and the reduction of highly water-intensive crops.