Anthropogenic and Climate-Induced Water Storage Dynamics over the Past Two Decades in the China–Mongolia Arid Region Adjacent to Altai Mountain

1. Introduction

Terrestrial water storage (TWS), which includes water stored in rivers, lakes, reservoirs, soil, groundwater systems, snow, glaciers and vegetation, plays a crucial role in the global hydrological cycle and land–atmosphere interactions [1]. TWS is essential for determining hydrologic transport and water availability [2], serving as a key indicator of water resource vulnerability [3]. However, accurately monitoring the dynamic changes in TWS, particularly in data-scarce regions, remains challenging due to the limited availability of field observation gauges, the high costs associated with installing and maintaining instrument networks, lack of data digitization and sharing, and significant uncertainties in hydrological and land surface models [4,5]. The China–Mongolia arid region adjacent to Altai Mountain (CMA) is a typical Asian inland arid zone, representing a key area of the Silk Road. The Junggar Basin of Northwestern China (JBNWC) and the Great Lakes Depression region of western Mongolia (GLDWM), which are located in the southern and northern Altai Mountains, are highly sensitive to global climate change, and their ecological environment is fragile [6,7]. The JBNWC has experienced population growth and urban expansion, which have increased the need for irrigation[8], while the GLDWM is primarily used for livestock farming, with mineral extraction also expanding in the past two decades [9,10]. Therefore, TWS plays a critical role in determining water availability, ensuring water and ecological security, and maintaining ecosystem functions in the CMA. Under the combined influence of climate change and human activities, monitoring the TWS changes in the CMA is of great significance for addressing regional water resource management and sustainable development.

The Gravity Recovery and Climate Experiment (GRACE) satellites provide an effective means of detecting TWS changes over large-scale regions, particularly in areas with limited in situ data [11]. To date, GRACE has been widely used for various applications, including identifying TWS variations at regional and global scales [11,12], assessing groundwater depletion [13,14], monitoring droughts [15,16], and detecting land subsidence [17]. However, GRACE data have some limitations that hinder their application for long-term hydrological changes at local scales, such as data gaps [18] and coarse spatial resolution [19]. Notably, there is an ~1-year gap between the original GRACE mission (launched in 2002) and its follow-on (GRACE-FO) mission, launched in 2018, which affects the continuity of monitoring water storage changes [20]. Furthermore, caution is necessary when applying GRACE observations to regions with heterogeneous geomorphological features [21]. Specifically, the GRACE observation footprint is approximately 200,000 km² due to the orbit height (~450×450 km = ~200,000 km²), which is too coarse to meet the detailed needs of water resource managers [22]. As a result, improving the spatial resolution of GRACE data holds significant potential for broadening their application.

There are two commonly used methods for GRACE downscaling: model- (dynamic) and data-based (statistical) approaches [23]. Model-based downscaling methods integrate GRACE observations into land surface, hydrological, or groundwater models based on physical processes, enabling the simulation of TWS at finer spatial resolutions. However, a key drawback is the high computational cost and strong dependence on model parameters and process formulations [24]. In contrast, statistical downscaling methods offer a simpler framework that requires fewer data. The basic assumption is that the relationships between TWS and other hydrological variables remain consistent across different spatial resolutions. By using the established relationships between TWS and other hydrological variables at the original spatial resolution, high-resolution data can be used to generate the downscaled TWS data [14]. Several statistical methods are commonly employed, including multiple linear regression (MLR), artificial neural networks (ANNs), and random forest (RF) models. Compared to other methods, the latter can simultaneously handle multiple input variables and exhibits greater robustness against overfitting [25,26].

Global Hydrological Models (GHMs), which simulate regional and global water fluxes and storage, typically include global land surface models (LSMs) and hydrological and water resource models (GHWRMs). These models provide valuable insights into global terrestrial water storage [11,27]. However, their outputs are often accompanied by large uncertainties, primarily due to limitations in model structure, data uncertainty, restricted water storage capacity, and the exclusion of human interventions (e.g., groundwater withdrawals) [20]. Those limitations are likely to contribute to the discrepancy between the TWS derived from GRACE and those from global models. Alternatively, LSMs are considered reliable tools for comparing and validating the accuracy of GRACE products [28]. Since GRACE cannot directly decompose TWS into its individual components, e.g., surface water storage (SWS), soil moisture storage (SMS), groundwater storage (GWS), snow water storage (SWE), canopy water storage (CWS), and glacier water storage (GlacierWS) [29], combining LSMs with GRACE has become an effective method for partitioning GRACE-derived TWS into specific components [30]. This approach is increasingly used to capture changes in GWS by separating SMS and SWE components in TWS, using auxiliary data from LSMs or GHMs [31]. Recently, there has been growing interest in using GRACE to enhance the simulation performance of GHMs and LSMs through data assimilation.

Climate change and human activities both influence hydrological cycle processes, which, in turn, affect variations in TWS [42,43]. Therefore, along with evaluating the TWS and GWS changes from GRACE, identifying and quantifying the effects of dominant controlling factors is critical for enhancing our understanding of and predictions for hydrological processes [42]. Recent research on TWS changes has increasingly focused on understanding the relationship between TWS and its individual drivers. Climate variables, such as precipitation and evapotranspiration, directly affect the amount of water available for TWS [44]. With climate warming, glacier retreat and reduced snow cover have led to increased meltwater, exacerbating TWS loss across entire basins [45]. Moreover, climate change—characterized by more frequent and intense extreme events—has increased the variability and uncertainty in precipitation, soil moisture and surface water, further impacting TWS [46]. In addition to climate change, human activities have become an increasingly important factor influencing the water cycle. Anthropogenic interventions, such as agricultural irrigation and coal mining, have been shown to accelerate TWS depletion, leading to regional-scale water storage loss [47,48,49,50]. Despite these findings, few studies have quantified the relative contributions of climate change and human interventions to TWS variations. Identifying the key factors driving TWS changes is crucial for improving regional water resource management and promoting sustainable development.

In conclusion, water resources in eco-fragile inland arid regions are highly sensitive to both climate change and human activities. It is critical to examine TWS and GWS dynamics across distinct sub-regions and catchments of the CMA and determine the relative contributions of climatic and anthropogenic drivers. Notably, in previous studies on GWS in Xinjiang, the presence of glaciers was often overlooked when estimating groundwater reserves in glaciated areas, resulting in overlap [32,33]. At the same time, most previous studies on the Mongolian Plateau have focused on TWS changes [40,41], while little is known about GWS changes and response to climate change. Although several studies have explored TWS or GWS changes in Xinjiang and the Mongolian Plateau [41,51,52,53,54], systematic comparative analyses of TWS and GWS dynamics in the CMA are still lacking. In the JBNWC, human activities like population growth, urban expansion, and rising irrigation demand are intensifying, while the GLDWM is more gradually impacted by human interventions, such as mineral extraction and animal husbandry. At present, spatiotemporal changes in TWS and GWS, along with their driving mechanisms in the CMA, are still poorly understood. This knowledge gap limits our understanding of the underlying processes driving TWS variations and hampers effective water resource management in the region.

This paper advances a systematic comparison of TWS and GWS dynamics across distinct sub-regions and catchments of the CMA at fine spatial scales, separating the relative contributions of climate and anthropogenic activities to TWS changes. The primary contributions of the paper are four-fold: (1) we enhanced the resolution of GRACE-based TWS to a 0.1°×0.1° grid scale by combining residual correction with the RF model; (2) considered the impact of glacial meltwater on GWS variability and included GlacierWS changes in TWS decomposition; (3) utilized the downscaled GRACE data to assess spatiotemporal trends in TWS and GWS at the pixel level across various basins in the CMA; and (4) quantified anthropogenic and climate-induced water storage trends in the CMA. High-resolution maps of TWS and GWS decline hotspots can inform local decision makers and stakeholders in addressing pressing water management issues by improving irrigation practices, managing water supply and demand, and mitigating the impact of droughts and human interventions. Understanding the climatic and anthropogenic drivers of water storage changes provides valuable insights into the impacts of climate change and human activities on water resources. This knowledge can inform future research and guide policies aimed at ensuring long-term water security.

2. Materials and Methods

2.1. Study Area

The CMA (79°52′ E–99°1′E, 43°4′N–50°57′N) includes the JBNWC and GLDWM; bordering Kazakhstan to the west and Russia to the north, its total area reaches approximately 69.37×km2. This region is in the hinterland of the Eurasian continent and has complex terrain and geomorphic features including the Altai, Tianshan, Sayan and Khangai ranges [55,56]. The CMA is one of the most climate-sensitive areas in the world, and it has become warmer and drier over the past few decades [38,55]. The JBNWC has a continental arid and semi-arid climate, influenced by the westerlies of the North Atlantic, with annual mean temperatures ranging from 5.13°C to 7.80°C. Winter temperatures average around -13°C, while summer temperatures reach approximately 22.2°C [56]. The annual average precipitation is 200 to 300 mm, with precipitation mainly concentrated in the summer months (June to August). Runoff in the JBNWC is mainly influenced by rainfall and meltwater from snow cover and glaciers, with summer runoff accounting for about 70% of the annual total. The GLDWM experiences an extreme continental climate, with annual mean temperatures ranging from -8°C to 2°C. Winter temperatures typically range from -20°C to -25°C, while summer temperatures vary from 15°C to 22°C [57]. The depression receives 100-150 mm of annual precipitation, while the Altai mountains receive over 300 mm (Tsegmid, 1969). Runoff in the GLDWM is mainly replenished by rainfall (25-40%), snowmelt (35-55%), groundwater (15-25%), and glaciers (<5%)[58].

Industrial development in the southern and northern Altai Mountains differs across the study area. Between 2002 and 2023, the JBNWC saw population growth, the expansion of irrigated areas, and rapid urbanization [33]. The largest land use type is bare land, followed by grassland and cultivated land (Figure 1). Though the region is predominantly industrial, it is also a key grain-producing and agricultural base [59]. In the GLDWM, bare land also occupies the largest area, followed by grassland and forest (Figure1). The region focuses on livestock development, and mineral resource extraction has expanded rapidly in recent years [9,60]. The fragile ecosystem and increasing human activities have heightened the demand for water resources in the region. Figure 1 shows the distribution of the 8 major river basins in CMA, including the Irtysh (IRB), Ulungur (URB), Emin (ERB), the north slope of Tianshan (NTRB), Khovd (KRB), Bulgan (BRB), and Zaykhan (ZRB) River Basins, as well as the Uvs Lake Basin (ULB).

2.2. Datasets

2.2.1. GRACE RL06 Mascon Solutions

GRACE data consist of two products [61]: (1) spherical harmonic (SH) coefficient products, derived from traditional basis function processing, and (2) mascon solution products, which are calculated based on Earth's gravitational changes, which are primarily influenced by TWS [31]. Compared with the SH solutions, mascon solutions have lower leakage errors and no post-process [62]. The GRACE(-FO) mascon data used in this study are based on three RL06 mascon solutions: those from the Jet Propulsion Laboratory (JPL-M) [63], the University of Texas at Austin (CSR-M)[18], and NASA's Goddard Space Flight Center (GSFC-M) [64]. The RL06 mascon solutions reduce signal noise and errors compared to RL05, with improved parameters, algorithms, data editing, and gravity models [65]. Note that the GRACE(-FO) mascon provides TWS anomalies (TWSAs) relative to the long-term mean from 2004 to 2009 [61].

However, different gravity field solutions increase data uncertainty. Some researchers have proposed that the simple arithmetic means of the three products can effectively reduce this uncertainty [32]. Considering that the grid resolution of the mascon solutions of JPL and GSFC was 0.5° × 0.5°, and while the CSR mascon products offer a grid resolution of 0.25° × 0.25°, the solution was resampled to a consistent resolution of 0.5° × 0.5° using TWSA. Then, we directly calculated the simple arithmetic mean of the three products. Due to its efficient, universal, and widespread use in handling missing data, cubic spline interpolation was used to fill in the missing data for 31 months from 2003 to 2023. After preprocessing steps, including format conversion, resampling, cropping, interpolation, and averaging, a TWSA dataset spanning 253 months was created.

2.2.2. Water Storage Components

The Global Land Data Assimilation System (GLDAS) can take advantage of satellite-and ground-based measurements, advanced land surface modeling, and data assimilation techniques in order to generate optimal land surface state and flux estimates. The GLDAS provides four land surface models (i.e., Noah, Mosaic, Variable Infiltration Capacity (VIC), and Catchment Land Surface Model (CLSM)), and is widely used in many fields. However, the Mosaic data are not available in the latest version of GLDAS. Therefore, in this study, we employed a suite of LSMs (Table 1) that are primarily derived from GLDAS versions 2.1 (Noah-3.6, VIC-4.1.2 and CLSM-OL) and 2.2 (CLSM-DA) to capture GWS changes by separating SMS and SWE components in TWS. Considering the inconsistent resolution of the four LSMs, to facilitate GRACE-derived TWS decomposition, the spatial resolution of the gridded data was resampled to 0.1° using a bilinear interpolation method. Furthermore, we used the LSM outputs to analyze the temporal variations in components. Note that CWS makes a small contribution to TWS compared to SMS in the SWE and SWS in the Asian inland arid region (Figure S1). In addition, the GlacierWS was estimated using the complete and consistent global dataset on glacier mass change [66]. These above components were transformed into anomaly data relative to the January 2004 to December 2009 time–mean baseline, which were then utilized to disaggregate and downscale GRACE-derived TWSA estimates. Meanwhile, TWS changes simulated by GLDAS during the 2001-2016 period were used for comparison with GRACE-derived TWS. More information about GLDAS can be found at https://disc.gsfc.nasa.gov/datasets.

2.2.3. Climate Variables and Anthropogenic Factors

To more accurately reflect the climate characteristics of the study area, we compared multi-source climate datasets and selected the best-performing products (Table 1, Figure S2). The temperature (Temp) data, with a spatial resolution of 0.1◦ and monthly temporal scale, are from the Climate Research Unit (CRU TS v4.05) at the University of East Anglia (https://crudata.uea.ac.uk/cru/data/hrg/). The precipitation (Pre) data used are mainly from the Global Precipitation Climatology Centre (GPCC) with a spatial resolution of 0.25◦ and monthly temporal scale (https://opendata.dwd.de/climate_environment/GPCC/html/fulldata monthly_v2022_doi_download.html). The evapotranspiration (ET) was derived from GLEAM satellites (https://www.gleam.eu/). In addition, other auxiliary data are also used in the downscaling of TWS, including land surface temperature (LST); potential evapotranspiration (PET); vegetation indicators such as Leaf Area Index (LAI) and Normalized Difference Vegetation Index (NDVI); topographic conditions such as Digital Elevation Model (DEM) and slope; soil texture variables such as sand, silt, and clay; and hydrological variables such as soil moisture (SM), snow water equivalent (SWE), Runoff (R), surface runoff (SurR), and subsurface runoff (SubR). The NDVI was established using MOD13C2 Version 6.1 with a monthly scale and 0.05◦ resolution (https://lpdaac.usgs.gov/products/mod13c2v061/). DEM and slope were measured using Shuttle Radar Topography Mission (SRTM) 90 m (https://earthexplorer.usgs.gov/). Soil texture variables were provided by the Harmonized World Soil Database version 2.0 (HWSD v2.0) (https://gaez.fao.org/pages/hwsd). The other auxiliary data used are mainly from ERA5-land reanalysis data (https://cds.climate.copernicus.eu/datasets/reanalysis-era5-land-monthly-means).

In addition, this study analyzes the impact of human factors on TWS by using indicators such as human water intensity, population, sown area, number of livestock and GDP. The human water intensity data, with a spatial resolution of 1km and an annual temporal scale, are from https://doi.org/10.6084/m9.figshare.10012559.v1. The data on population, sown area, number of livestock and GDP were gathered from the statistical yearbook and National Statistics Office of Mongolia at an annual temporal resolution (https://www.1212.mn/en/statistic/statcate/573052/table-view/DT_NSO_0500_007V1).

2.3. Methods

2.3.1. Partial Least Squares Regression

While the quality of predictors plays a crucial role in the performance of downscaling models, there is no standardized criterion for selecting input variables. To address this, we first compiled a comprehensive list of predictors used in previous GRACE downscaling studies [69,70,71]. These input variables can be broadly categorized into several types, including climatic indices (e.g., Pre and Temp), water storage components (e.g., SM and SWE), and topographic factors (e.g., lithology and hydraulic conductivity), among others.

Partial least squares regression (PLSR) was used to assess the contribution of various predictors to GRACE-based TWSA in the CMA. PLSR is a multivariate data analysis technique based on principal component analysis (PCA). Unlike PCA, which focuses solely on explanatory variables, PLSR incorporates information from both the explanatory and response variables, allowing for the effective extraction of the most influential predictors. For instance, Yuan et al. used PLSR to identify key factors driving TWS changes in different global regions [72]. To measure the importance of each predictor, the variable importance in projection (VIP) score was calculated. A higher score indicates a greater influence on TWSA, with predictors with values greater than 0.8 typically considered highly influential in explaining the dependent variable [73]. The VIP score for the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ variable was calculated as follows [19]:

$${\mathrm{V}\mathrm{I}\mathrm{P}}_j={\left[{\displaystyle \frac{p·{\sum }_{m=1}^{h}R\left(Y,t_m\right)·w_{mj}^2}{{\sum }_{m=1}^{h}R\left(Y,t_m\right)}}\right]}^{{\displaystyle \frac{1}{2}}}$$

(1)

where h is the number of components, p is the number of predictors, ${\mathrm{t}}_{\mathrm{m}}$​ is the ${\mathrm{m}}^{\mathrm{t}\mathrm{h}}$ component, R(Y,${\mathrm{t}}_{\mathrm{m}}$) is the explanatory power of ${\mathrm{t}}^{\mathrm{m}}$​ for Y, and ${\mathrm{w}}_{\mathrm{m}\mathrm{j}}$ is the marginal contribution of the predictor.

2.3.2. Random Forest Model

Random forest (RF) is an ensemble method based on decision trees, enhancing prediction accuracy and stability [74]. In regression tasks, it creates multiple subsets via bootstrapping, trains a regression tree for each, and combines their predictions into the final model. The algorithm used in this study is summarized as follows:

Step 1: Model parameter optimization and low-resolution modeling. For each month, 80% of the low-resolution (0.5°) TWSA and auxiliary data are used for training, and 20% for validation. Two key random forest parameters are optimized: the number of trees (ntree, ranging from 50 to 1000 in steps of 50) and the number of input variables per node (mtry, ranging from 1 to 10). Other parameters were set to default values. By iterating over ntree and mtry 360 times, the parameter combination with the smallest prediction error was selected for optimization. Based on this optimal combination, a random forest regression model for the low-resolution scenario of each month was established.

Step 2: High-resolution scenario prediction. Assuming that the relationships between variables do not change with the spatial resolution, the statistical regression model established under low-resolution conditions is also applicable to high-resolution scenarios. The high-resolution (0.1°) independent variables are input into the random forest regression model to initially obtain the TWSA at 0.1° resolution.

Step 3: Residual correction. The independent variables at 0.5° are inputted into the model to obtain the predicted values of TWSA. The difference between the actual and predicted values of TWSA is the residual, representing the portion that the model cannot explain (mainly due to human activity). Using bilinear interpolation, the residual at 0.5° is resampled to obtain that at 0.1°. The 0.1° TWSA obtained in step (2) is then added to the 0.1° residual, yielding the residual-corrected 0.1° TWSA.

2.3.3. Decomposition of TWSA

We derived GWS anomalies (GWSAs) from GRACE TWSA by subtracting GLDAS land water storage and GlacierWS anomalies (GlacierWSA); the latter were calculated using glacier mass change data. The GRACE-derived GWSA was validated through comparison with groundwater well observations. Compared with in situ groundwater level measurements, the overall correlation coefficient (CC) and root mean squared error (RMSE) of 0.1° GWSA with residual correction are 0.437 and 0.087, respectively. More descriptions of GWSA are included in the supporting information (Text S2).

$$GWSA=TWSA-SMSA-SWEA-SWSA-GlacierWSA$$

(2)

where SMSA, SWSA, SWEA, GWSA and GlacierWSA denote anomalies in soil moisture, surface water, snow water equivalent, groundwater and glacier water storage, respectively.

2.4.4. Quantifying the Relative Contributions of Climate and Human Activities to Changes in TWSA

To quantify climatic and anthropogenic impacts on TWS changes, this study applied the seasonal-trend decomposition procedure based on Loess (STL), one of the effective methods for refining and highlighting variations in time data series at different levels. The STL divides a time data series into (1) a trend component, (2) a seasonal component, and (3) a residual term [75]. This method has been used to clarify the properties of the TWS and examine the relative importance of the resulting features of GRACE-derived TWS [76]. Additionally, to identify the small variations in monthly TWSA, we adapted the STL as recommended [77]. Lastly, the time series are represented as a sum of a seasonal cycle (${\mathrm{S}}_{\mathrm{i}}$), residual term (${\mathrm{R}}_{\mathrm{i}}$), linear trend (${\mathrm{T}}_{\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{i}}$) and inter-annual variability (${\mathrm{T}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}}$), whereas the linear trend and inter-annual items were derived from the long-term trend component, as summarized in Eq. (3). Both the trend component and linear trend were used to perceive the magnitude of change in TWSA and its components. More details can be found in references [75] and [77]. According to the authors of [78], we defined the inter-annual and seasonal components of TWSA as the natural TWSA, while defining the linear trend and sub-seasonal residuals of TWSA as the anthropogenic TWSA.

$$Y_i=T_i+S_i+R_i=T_{{linear}_i}+T_{{inter}_i}+S_i+R_i\left(i=\mathrm{1,2},\cdots ,N\right)$$

(3)

3. Results

3.1. Evaluation of Downscaled Terrestrial Water Storage Anomalies

3.1.1. Selection of Representative Predictors for Downscaling Models

Following the initial screening, variables identified in previous GRACE downscaling studies [19,79,80] were chosen as the primary drivers influencing the downscaling process, as illustrated in Figure 2. It is important to note that human activities have a significant impact on TWSA within the CMA. However, due to the lack of comprehensive, high-resolution datasets on agricultural management over extended periods, the auxiliary data utilized in this study consist solely of natural variables. The effects of human activities are indirectly accounted for through adjustments made during the residual correction in the downscaling procedure.

According to [81], this study selected 17 natural variables, including meteorological variables such as Temp, LST, Pre, ET, and PET; vegetation indicators such as LAI and NDVI; topographic conditions such as DEM and slope; soil texture variables such as sand, silt, and clay; and hydrological variables such as SM, SWE, R, SurR, and SubR. The seventeen selected variables, along with GRACE-derived TWSA, serve as the input and target variables for the PLSR model. The significance of each variable in relation to TWSA is presented in Figure 2. From this, it is evident that the VIP values for SM, SWE, Temp, LST, NDVI, DEM, Pre and SubR exceed 0.80, indicating that these eight variables are the most influential for the CMA. As a result, these eight variables are chosen as the input features for the RF-based downscaling models.

3.1.2. Comparison of Downscaled TWSA and GWSA Using Residual Correction

Figure 3 suggests the spatial variability of TWSA and GWSA trends before and after downscaling during the 2002–2023 period. The overall TWSA trend of the downscaled data and that after residual correction are consistent with the original GRACE observations, with CC metrics of 0.82 and 0.91, respectively (Figure 3a–3c). The downscaling scheme can capture downtrends in the JBNWC, which are consistent with original GRACE data. However, without residual correction, the TWSA decline rate in the JBNWC will be significantly underestimated (Figure 3b). The downscaling result after residual correction is highly consistent with the original data in both numerical size and spatial distribution pattern (Figure 3c), indicating that the downscaling model after residual correction in this paper is relatively reliable. It is worth noting that the downscaling scheme has a better ability to describe the local details of TWSA changes. The downscaling results identified a large rate of decline in the GLDWM (Figure 3b–3c), which was not evident in the original GRACE observations. There are glaciers in the GLDWM (Figure 1), so the downscaling results more finely identify glacier retreat-induced declines in TWSA.

The spatial distribution of GWSA trends in the JBNWC closely mirrors that of TWSA trends, with a CC metric of up to 0.95, although the trend amplitudes differ. This suggests that the inter-annual variations in TWSA in the JBNWC are primarily driven by GWSA, which aligns with findings from previous studies [32]. In contrast, the GWSA trends in the GLDWM exhibit less similarity to those of the TWSA, with a CC of 0.78. This indicates that human activities have less impact on TWSA and GWSA in the GLDWM.

3.2. Temporal Variations in TWSA and GWSA

3.2.1. Multiscale Characteristics in TWSA and GWSA in the CMA

Using downscaled data derived from the RF model with residual correction, the temporal variations in TWSA and GWSA from 2002 to 2023 were obtained, as shown in Figure 3. The monthly TWSA in the CMA experienced a significant decline (-3.32 mm/a, p < 0.01) over the study period, with an accelerated decrease observed after 2017 (-6.36 mm/a, p < 0.01) (Figure 4a). In terms of intra-annual variability, spring snowmelt, which peaks in April, rapidly replenishes surface water reserves, resulting in the highest TWSA during this period. In autumn, due to less precipitation, higher evapotranspiration and the continuous increase in human water withdrawal, TWSA declines to its lowest point in October. This pattern corresponds to the seasonal variations in TWSA observed in the CMA (Figure S3), which is a result of the joint consequences between natural recharge mechanisms and human-induced water consumption. Meanwhile, Figure 4b illustrates the monthly GWSA variations in the CMA. The inter-annual trend of GWSA (-2.71 mm/a, p < 0.01) is slightly lower than that of TWSA over the entire study period. However, since 2017, the rate of GWSA decline (-9.38 mm/a, p < 0.01) has exceeded that of TWSA, suggesting that natural recharge processes can no longer meet the growing demand for water, leading to an increasingly severe overexploitation of groundwater resources. In addition, due to the delayed response of groundwater to seasonal variations and the inherent stability of the system, the intra-annual variability of GWSA remains relatively steady.

3.2.2. Temporal Variations in TWSA, GWSA and Its Components in the JBNWC and GLDWM

We further analyzed the trends in TWSA, GWSA, and their components across the JBNWC and GLDWM regions, as shown in Figure 5. Between 2002 and 2023, TWSA in the GLDWM declined at a rate of -1.11 ± 0.12 mm/yr (equivalent to 6.67 Gt). The slight negative correlation between precipitation and TWSA suggests that precipitation is not the primary driver of TWSA changes. Interestingly, peak precipitation coincided with low TWSA values, likely due to high evapotranspiration during the summer months (Figure 5c). The TWSA in JBNWC decreased by 77.2 mm (15.66 Gt) from 2002 to 2023, with the highest TWSA value in April 2002. TWSA decline from April 2002 to May 2023 occurred in three phases. From April 2002 to October 2009, TWSA declined significantly at -5.78 ± 0.45 mm/year (p < 0.05), with a total decrease of 58.4 mm. Between May 2010 and September 2015, the TWSA decline continued, though at a slightly lower rate (-4.38 ± 0.70 mm/yr, p < 0.05), leading to a reduction in equivalent water height by 52 mm. Finally, from March 2017 to May 2023, TWSA showed a slight increasing trend, with an equivalent water height decrease of 53.25 mm at a rate of -9.31 ± 0.79 mm/yr (p < 0.05).

Since snow water, glacier meltwater, soil moisture, surface water and groundwater affect TWSA (with canopy water being negligible in the CMA (Figure S2)), we estimated each component's contribution to TWSA in the JBNWC and GLDWM regions separately. In both the JBNWC and GLDWM, GWSA is the primary contributor to TWSA variance, accounting for 78.16% and 90.28%, respectively, followed by SWEA at 18.61% and 6.04%. This indicates that TWSA variation is largely driven by fluctuations in groundwater reserves. A comparison of the GWSA, SWEA and TWSA time series reveals the three-month lag effect of GWSA and SWEA on TWSA in the JBNWC; this is due to delayed effects from snowmelt and groundwater recharge processes. In the GLDWM, although SWEA occasionally aligned with TWSA, a time lag of approximately 1–2 months was identified. In contrast, no lag was observed between GWSA and TWSA, as limited surface water resources in the GLDWM indicate that TWSA changes primarily reflect groundwater fluctuations. Our results showed that glaciers in the JBNWC and GLDWM experienced accelerated shrinkage during the study period, with glacier mass balance rates of -0.001 m w.e./yr and 0.01 m w.e./yr (Figure 8b-Figure 9b), respectively, which also explains the negative TWSA. These estimates were derived from a globally complete and consistent analysis of glacier mass change based on satellite observations, aligning with findings from previous research [57,82].

3.3. Spatial Dynamics of TWSA and GWSA

The spatial distribution of TWSA in the CMA, derived from three GRACE mascon solutions, reveals a consistent pattern, with the lowest values located in the southwestern region of the study area (Figure S4). As shown in Figure 6a, the spatially averaged annual mean TWSA derived from GRACE(-FO) indicates that the majority of the CMA (89.11%) experiences negative TWSA over the study period, reflecting TWSA deficits. The NTRB of the JBNWC exhibits the lowest TWSA at -147.34 mm. A significant decreasing trend in TWSA was detected across 98.61% of the CMA during this period (Figure 6b-6c, Figure S5), with the NTRB of the JBNWC showing the most substantial decline at a rate of -28.39 mm/yr. Compared to the JBNWC, the TWSA in the GLDWM exhibits a relatively uniform spatial distribution, with minimal spatial heterogeneity. Notably, the glacier-covered area in the KRB shows the most significant decreasing trend, at 13.37 mm/year. These findings suggest that glacier retreat in Mongolia's Altai Mountains has accelerated in recent years, likely driven by climate warming [83]. In addition, positive TWSA values are observed in the IRB, URB and ERB of the JBNWC, as well as in the KRB and BRB of the GLDWM, with the IRB of the JBNWC showing the highest TWSA at 30.38 mm. Notably, only the IRB and BRB exhibit significant increasing TWSA trends, at rates of 3.67 mm/yr and 2.84 mm/yr, respectively (Figure 6b-6c). Unlike TWSA, the maximum values and increasing trends of GWSA are observed in the KRB and BRB, whereas the GWSA in the IRB and URB shows negative values with a declining trend (Figure 6d–6f). This indicates that the rapid socioeconomic development in the JBNWC has led to a higher water demand, resulting in a significant reliance on groundwater resources. In contrast, in the GLDWM, glacier meltwater significantly recharges groundwater, resulting in an increase in groundwater storage.

As presented in Figure 7, the trends in TWSA and GWSA were analyzed across four seasons—spring, summer, autumn and winter. NTRB in the JBNWC experienced a significant decrease in TWSA and GWSA during spring, with the lowest values of -28.14 and -33.35 mm/yr, respectively, followed by summer and autumn, and decreased the least in winter. This highlights that agricultural irrigation in the JBNWC requires substantial water resources. However, the limited glacier-snow meltwater and precipitation are insufficient in meeting irrigation demands, making groundwater the primary source. Compared to the JBNWC, the GLDWM exhibits smaller seasonal variations in TWSA. However, the increases in GWSA within the KRB and BRB are most pronounced during winter and spring, followed by summer and autumn. This is attributed to the combined effects of glacier-snow meltwater and lower evaporation during winter and spring, which enhance groundwater recharge. In contrast, higher evaporation rates and increased human water demand during summer and autumn lead to relatively lower groundwater recharge.

3.4. Climatic and Anthropogenic Contributions to TWSA Dynamics

3.4.1. Climatic Factors

To assess the influential factors of climate-driven TWSA, multiple climate variables in the JBNWC and GLDWM have been analyzed. Figure 8 and Figure 9 depict the temporal variations in monthly TWSA, in comparison to precipitation, glacier and ET. The TWSA peak in the JBNWC occurs in April, while in the GLDWM, it remains relatively high throughout the winter and spring. In both regions, TWSA peaks occur earlier than the precipitation peak in July. Additionally, the summer evapotranspiration peak aligns with the precipitation peak, offsetting much of the water recharge from precipitation. This suggests that summer precipitation contributes only limited direct recharge to TWSA in the CMA region. In contrast, TWSA variations show a strong correlation with changes in glacier mass balance, indicating that meltwater from glaciers and snowmelt serves as the primary source of TWSA recharge in the CMA region[37,84]. Notably, the TWSA in JBNWC remains negative throughout the year, indicating a persistent water deficit in the region. This suggests that winter precipitation, glacier and snowmelt recharge are insufficient in meeting the region's water demand, thereby exacerbating groundwater depletion. Moreover, climate warming accelerates glacier retreat and increases evapotranspiration, further exacerbating the potential threats to regional water resources. Compared to the JBNWC, the TWSA in the GLDWM exceeds zero during the winter and spring, indicating that the region temporarily achieves water resource balance during these seasons. However, this balance is highly dependent on glacier meltwater recharge, and the long-term decline in glacier mass balance threatens to disrupt this seasonal equilibrium, posing significant risks and challenges to regional water resources in the future.

Figure 8. Spatial distribution of TWSA trends (left panel) and GWSA trends (right panel) over four seasons: (a) and (e) MAM, (b) and (f) JJA, (c) and (g) SON and (d) and (h) DJF.

Figure 8. Spatial distribution of TWSA trends (left panel) and GWSA trends (right panel) over four seasons: (a) and (e) MAM, (b) and (f) JJA, (c) and (g) SON and (d) and (h) DJF.

Figure 9. Interannual and intra-annual variations of TWSA, precipitation, Glacier mass balance, and ET in JBNWC. GMB refers to glacier mass balance.

Figure 9. Interannual and intra-annual variations of TWSA, precipitation, Glacier mass balance, and ET in JBNWC. GMB refers to glacier mass balance.

3.4.2. Anthropological Factor

To identify the anthropogenic drivers of TWSA variations, we analyzed the differences between GRACE and GLDAS data for the JBNWC and GLDWM. Figure 10a and Figure 10b illustrate these differences in the two regions, with RMSE values of 32.5 and 16.4, respectively. Notably, the TWSA trends from GRACE and GLDAS show relatively good consistency in both the JBNWC and GLDWM prior to 2010 (pre-2010), compared to the period thereafter (post-2010) (Figure 10a and Figure 10b). Since GLDAS primarily simulates TWSA variations driven by natural factors [100], this suggests that natural factors were the dominant contributors to TWSA changes before 2010. However, a significant increase in the inconsistency between GRACE and GLDAS is observed post-2010, with the RMSE increasing by 88.5% in the JBNWC and 34.5% in the GLDWM (Figure 10c). This substantial rise in discrepancy implies a shift in the dominant drivers of TWSA variations, possibly indicating the growing influence of anthropogenic factors. In the pre-2010 period, the TWSA trends derived from GRACE and GLDAS were -3.87 mm/yr and -0.25 mm/yr in the JBNWC, and -2.31 mm/yr and -1.42 mm/yr in the GLDWM, respectively (Figure 10d and Figure 10e). The underestimation of the declining TWSA trend in GLDAS is likely due to its exclusion of anthropogenic influences and its limited representation of water storage capacity [85]. This suggests that the observed TWSA decline during this period was primarily driven by climatic factors. In contrast, after 2010, the TWSA trends derived from GRACE showed a continued decrease (-3.95 mm/yr and -0.42 mm/yr) in both the JBNWC and GLDWM, while GLDAS showed an increasing trend (0.91 mm/yr and 1.1 mm/yr). Given that GLDAS primarily reflects climate-driven changes in TWSA, these results imply that the TWSA decline after 2010 was likely more influenced by anthropogenic activities than by natural factors.

Additionally, we analyzed the contributions of water withdrawals to shifts in the driving factors of TWSA changes. Figure 10f and Figure 10g show that water withdrawals increased by 10.19% in the JBNWC and 14.79% in the GLDWM from pre- to post-2010. This significant rise in water withdrawals further underscores the growing dominance of anthropogenic activities in the decline of TWSA, suggesting a shift from climate- to human-driven factors. It is important to note that, following the government’s implementation of efficient water-saving agricultural practices [86], water withdrawals in the JBNWC have been effectively controlled and gradually reduced since 2010. In contrast, although water withdrawal intensity in the GLDWM remains relatively low, the intensity of agricultural development, particularly in terms of water and land resource use, has been steadily increasing in recent years [87], leading to a continuous rise in water withdrawal intensity. Given the region's ecological vulnerability, the combined impacts of climate change and human activities are likely to exacerbate the conflict between water supply and demand in the future.

3.4.3. Relative Contributions of Drivers to the TWSA Trend

In this study, we employed STL decomposition to separate the inter-annual and seasonal components of TWSA as natural TWSA, while defining the linear trend and sub-seasonal residuals of TWSA as anthropogenic TWSA. Figure 11 illustrates the relative contributions of climate change and human activities to TWSA trends in the JBNWC and GLDWM over the 2002–2023 period. During this time, the natural TWSA variations in the JBNWC displayed fluctuating patterns, with climate change contributing 38.12% to the overall TWSA trend. In contrast, anthropogenic TWSA showed a marked declining trend (-0.26 mm/month), with human activities contributing 61.88% to the TWSA reduction. This significant decline highlights excessive groundwater depletion due to continued over-extraction, which was not adequately compensated by runoff from precipitation or glacier meltwater. In the GLDWM, natural TWSA variations were relatively smaller, with climate contributing 50.30% to the TWSA trend. The anthropogenic TWSA in this region also exhibited a downward trend (-0.09 mm/month), with human activities accounting for 49.70% of the TWSA reduction. Further analysis of the periods before and after 2010 revealed that, prior to 2010, TWSA decline in the JBNWC was driven by both climate change and human activities. However, in the GLDWM, climate change was the dominant driver (62.05%), with a lesser contribution from human activities (37.05%). Since 2010, the contribution of human activities to TWSA reduction has increased significantly, accounting for 70.36% in the JBNWC and 57.98% in the GLDWM, solidifying anthropogenic influences as the main factor driving TWSA decline (Figure 11c, f).

4. Discussion

4.1. Effect of Glacier Cover on GWSA

In glacier-covered regions, most previous studies estimating GWSA using the water storage components equation have directly ignored the presence of glaciers, resulting in groundwater storage change signals that are confounded with glacier change signals. To investigate the impact of glacier coverage on GWSA estimates, we examined the inter-annual trends of GWSA in different basins before and after glacier removal. There is no glacial development in the ERB, the Gurbantunggut Desert of the JBNWC, or the ULB of the GLDWM, so the trends in GWSA before and after glacier removal are identical (Table 2). For other basins, without removing glaciers, the GWSA estimates in glaciated areas actually represent the combined changes in groundwater and glaciers. For the NTRB of the JBNWC, since groundwater extraction is primarily concentrated in the plain and basin areas, the groundwater storage changes in the non-glaciated areas largely reflect those across the entire basin. In recent years, the Tianshan Mountains have experienced a significant loss in glacier mass, with melt causing a continuous decline in terrestrial water storage in the region [43,88]. As a result, the absolute rate of change in GWSA in the NTRB is slightly higher than that observed after removing glaciers, leading to an overestimation of 1.85%. In contrast, for the IRB of the JBNWC and the KRB of the GLDWM, located in the southern and northern slopes of the Altai Mountains, groundwater extraction is relatively low, and glacier recharge compensates for part of the groundwater depletion, resulting in underestimations of 24.98% and 17.78%, respectively. Additionally, although the glacier coverage in the URB of the JBNWC and the BRB and ZRB of the GLDWM are relatively small, there are still underestimations of 1.73%, 4.48%, and 4.42%, respectively. In terms of spatial variability, the spatial distributions of GWSA with and without glacier removal are nearly identical. Compared to the former, the GWSA decline rate in the glaciated areas of the NTRB is lower, while the increase rate in the KRB is higher, both highlighting the contribution of glaciers to groundwater recharge (Figure 12).

Figure 12. The components of TWSA decomposed by STL in (a) JBNWC, (d) GLDWM; the anthropogenic activities dominated–TWS changes and climatic features dominated–TWS changes in (b) JBNWC, (e) GLDWM; Climatic and anthropogenic contributions on TWSA for the period 2002–2023 and subperiods 2002–2010 and 2011–2023 in (c) JBNWC, (f) GLDWM. EWT refers to equivalent water thickness.

Figure 12. The components of TWSA decomposed by STL in (a) JBNWC, (d) GLDWM; the anthropogenic activities dominated–TWS changes and climatic features dominated–TWS changes in (b) JBNWC, (e) GLDWM; Climatic and anthropogenic contributions on TWSA for the period 2002–2023 and subperiods 2002–2010 and 2011–2023 in (c) JBNWC, (f) GLDWM. EWT refers to equivalent water thickness.

Figure 13. Spatial distribution of interannual trends in GWSA before and after removal of glaciers.

Figure 13. Spatial distribution of interannual trends in GWSA before and after removal of glaciers.

4.2. Drivers of Water Storage Changes

Anthropogenic and climatic influences on changes in water storage are complex and non-linear, we compared the relative contributions of these two factors to TWS before and after 2010. Although annual precipitation in JBNWC is increasing (1.19 mm/yr), a larger rise (2.09 mm/yr) in evapotranspiration offsets the trend toward wetness (Figure 8a, e), leading to a warmer, drier climate [55]. With the implementation of China's Western Development Strategy, the economy of the JBNWC began to grow rapidly. From 2000 to 2010, the region's GDP was approximately CNY 1469.031 billion, and the sown area was 911,000 hectares (Figure S6). Economic expansion was largely driven by the extraction of abundant natural resources, including oil, natural gas and minerals, while agriculture and animal husbandry also played a supportive role in economic development. During this period, increased evaporation-induced drought played a dominant role in TWS decrease, with a climate-induced contribution rate of 55.71% (Table S1), relative to that of human activities (31.84%). From 2010 to 2021, the JBNWC experienced substantial economic growth, with GDP soaring by approximately 308.64% and sown area nearly doubling. Water withdrawal intensity rose by 10.19% (Table S1). To prevent the over-extraction of groundwater, the government implemented a series of policies and regulations, including shutting down illegal or non-compliant groundwater wells. Despite these efforts mitigating the decline in groundwater levels, the surging economy may continue to escalate water demand. Consequently, human activities contribute 70.36% to TWS reduction. Meanwhile, from 2010 to 2021, intensifying droughts accelerated glacier melting, replenishing regional water resources and reducing the contribution of climate change to TWS reduction to 27.54% (Table S1). The annual precipitation increase in GLDWM (1.42 mm/yr) slightly exceeds that of ET (1.34 mm/yr), but since annual ET surpasses annual precipitation (Figure 9a,e), the GLDWM is also becoming warmer and drier [89]. Compared to the JBNWC, the GLDWM's economic development remains in its early stages, relying on livestock farming and resource extraction. From 2000 to 2010, climate change drove 65.21% of TWS reduction (Table S1). To ensure food security, Mongolia has expanded arable land significantly since 2010, with sown areas and fodder crop cultivation increasing alongside rising livestock numbers to address extreme climate conditions. Industrial growth has also accelerated due to intensified mineral extraction. Since 2010, human activities have dominated TWS reduction, contributing 58.22%, which is slightly above climate change's 42.68% (Table S1).

Under the Belt and Road Initiative, the strategic importance of the CMA is growing, but water scarcity continues to be a significant constraint on the region's economic development. A thorough understanding of the water shortages in the area is crucial, as is the need to propose effective water resource management strategies. For the GLDWM, potential solutions include regulating livestock populations, promoting industrial diversification, and expanding the use of water-saving irrigation technologies. In the JBNWC, building a water-conserving society and implementing inter-basin water transfer projects offer feasible strategies to mitigate the water crisis.

4.3. Limitations and Future Work

This study does not explicitly account for the effects of mineral extraction and reservoir operations on TWS and GWS variability, which may introduce mass loss interference and lead to an overestimation of TWS changes. The lack of long-term, high-resolution datasets on mineral resource extraction limits our ability to accurately assess the spatiotemporal variation in TWS on a monthly scale. While the impact of mineral extraction on TWS trends is relatively minor, removing its influence could improve the accuracy of seasonal estimates. Future studies could explore potential vulnerabilities by projecting TWS using advanced climate models (e.g., CMIP6) to gain a deeper understanding of how water storage changes might affect food and water security in the CMA.

5. Conclusions

We examined GRACE-based water storage trends and possible drivers of variation to gain insights into how climate and human factors are linked to water storage variability. We used the RF model combined with residual correction, which exhibited superior performance compared to the RF model alone in predicting TWS and GWS anomalies. While TWS and GWS have durably decreased since 2002 over the majority of the CMA, the declining trends were more pronounced in areas with higher agricultural and human settlement impacts. We quantified the relative contributions of human activities and climate change to TWS dynamics and examined the influences of their driving factors. Although increased precipitation positively impacts TWS, the rise in evapotranspiration may accelerate its decline, partially offsetting these benefits, while glacier melting driven by climate change further accelerates TWS depletion. In the JBNWC, human activities such as increased sown area, water consumption, and settlements played dominant roles in TWS decrease. Compared with climate change (38.12%), the human-induced contribution rate was 61.88%, which increased to 70.36% after 2010. In the GLDWM, the primary driver of TWS decline shifted from climatic (2002-2010) to anthropogenic factors (2011-2020), while the contribution of human activities—such as the expansion of mineral extraction and sown area—increased from 37.95% to 57.98%. In addition, this study assessed the impact of glacial meltwater on GWS estimates, showing that neglecting it can lead to over- or underestimations across basins. Introducing advanced technology and controlling livestock in the GLDWM, while developing a water-saving society and implementing inter-basin water transfers in JBNWC, are viable strategies in addressing the water resource crisis. The presented framework is applicable to other regions for assessing local water resource vulnerabilities and developing effective management strategies.