Introduction

Water resources are essential for ecosystem health and socioeconomic sustainable development, especially in arid and semi-arid sandy regions with limited supply1,2,3. In particular, the MuUs, Otindag, Horqi, and Hulun Buir Sandy Lands, as a vital ecological protective screen of China (Fig. 1a), are of considerable importance for the ecological security of northern China and even the whole country. Previous studies mainly focused on the land degradation and desertification measures in these regions4,5,6, with very limited understanding on water resources variability.

Fig. 1: TWS, groundwater storage, and soil moisture storage variations during 2002–2023.
figure 1

a geographical information of the ___domain of interest. The global map data was provided by Google Earth. b water volume change over 2002–2023 (water volume change =  trend slope × number of years × area of sandy land). Error bars represent 95% confidence intervals for each water storage trend. c–f trends of water storage over the sandy lands. The shaded blue area represents 1σ error of TWS. The dashed line represents the best fit linear trend. Slopes with and without asterisks indicate the slopes are significant and insignificant at 0.01 level, respectively.

Full size image

While existing studies have shown changes in lake area7,8, water volumes9,10,11,12, and groundwater levels13 in sandy land regions, they often focus on individual components of water resources, lacking a comprehensive understanding of total water resources changes. In recent years, some studies have attempted to investigate the patterns and possible causes of total water resources changes in sandy areas with the development of satellite technology14,15. However, two key issues remain. First, there is an incomplete consideration of all potential causes. For example, a study by Zhang et al. 14 found a decrease in Total Water Storage (TWS) in the Otindag Sandy Land and surrounding areas, linking it to human water use. However, this study did not take into account the water consumed by vegetation. Another recent study15 observed a significant TWS decline over MuUs Sandy Land, which was attributed to the ecological restoration programs. However, the study did not fully explore other potential factors, such as agricultural, industrial, and domestic water use, all of which have increased significantly in these regions during 2002–2023 (Supplementary Fig. 1), leading to unneglected impacts on TWS changes. Second, there is a lack of comprehensive analysis on future TWS changes in sandy regions, which serve as China’s ecological protective screen. There is no study that assessed the future projections of TWS changes, as well as the potential impacts on the ecological security and sustainable development. Although Zhao et al.15 used scenario analysis to explore whether the current TWS decline in MuUs sandy land would continue under different intensities of ecological restoration policies, they did not quantify the exact magnitude of TWS changes. Additionally, they did not thoroughly examine the potential impacts of TWS changes on ecological security and socio-economic development. Furthermore, research on future TWS changes in China’s other sandy regions remains limited.

To fill the above research gap, this paper aims to (1) reveal the changes in TWS and its key components across the China’s sandy lands from 2002 to 2023 by integrating multi sources dataset; (2) clarify the main drivers of TWS changes from both aspects of natural and anthropogenic disturbances; and (3) project the future TWS changes and elucidate its implications to sustainable development and ecological security, providing a scientific basis for timely water resources management and ecological restoration efforts towards sustainable development.

Results

Groundwater decrease reduced TWS

In last two decades (2002–2023), all the sandy lands, located in arid and semi-arid regions (MuUs, Otindag, and Horqin), experienced 55.97 billion m3 net water loss in total (approximately 2 times of the water storage of the Poyang Lake, the largest freshwater lake in China; Fig. 1b). The fastest decline is found in MuUs at a rate of −9.03 ± 0.22 mm yr−1 (P < 0.01, Fig. 1c), resulting in a water loss of 30.12 ± 0.74 billion m3 (Fig. 1b), followed by Otindag ( − 4.43 ± 0.15 mm yr−1, P < 0.01; Fig. 1d; 14.74 ± 0.50 billion m3, Fig. 1b) and Horqin (−4.67 ± 0.44 mm yr−1, P < 0.01, Fig. 1e; 15.54 ± 1.45 billion m3, Fig. 1b). Hulun Buir Sandy Land, located in the semi-humid area, was the only exception, where the TWS exhibited a slow increase (2.24 ± 0. 36 mm yr−1, P < 0.01; Fig. 1f), resulting in a total water gain of 4.43 ± 0.58 billion m3 (Fig. 1b).

Analysis of key TWS components showed that the observed net water loss was mainly dominated by the groundwater decline. Specifically, groundwater storage decreased (Fig. 1c–f) at rates of −9.64 ± 0.21 mm yr−1, −6.61 ± 0.16 mm yr−1, −8.11 ± 0.27 mm yr−1, and −2.62 ± 0.18 mm yr−1, respectively, in the MuUs, Otindag, Horqin and Hulun Buir Sandy Lands, causing 32.14 ± 0.69, 21.99 ± 0.52, 26.97 ± 0.90 and 5.18 ± 0.30 billion m3 groundwater storage loss (Fig. 1 b). The substantial groundwater decrease was partially counteracted by the increase of soil moisture at rates of 0.68 ± 0.16 mm yr−1, 2.31 ± 0.19 mm yr−1, 5.17 ± 0.40 mm yr−1, and 3.82 ± 0.30 mm yr−1 (Fig. 1c–f), resulting in about 2.27 ± 0.52, 7.69 ± 0.63, 17.18 ± 1.32 and 7.56 ± 0.48 billion m3 soil moisture gains correspondingly (Fig. 1b).

Intensified human activities cause water storage loss

The attribution analysis (see Methods) indicates that, overall, the primary causes of water loss are human water consumption and vegetation evapotranspiration (ET) resulting from ecological restoration. Over the past 20 years, climate change has caused TWS in the study area to increase by 31.67 billion m3 (Fig. 2a). However, increased human water consumption (Supplementary Fig. 1)—including agricultural (49.56 billion m3), industrial (7.86 billion m3) and domestic water consumption (3.94 billion m3)—has resulted in a net water loss of 61.36 billion m3 (Fig. 2a), offsetting the climate-driven gain. Additionally, ET from ecological restoration projects has increased at a rate of 1.68 ± 0.29 billion m³ per year (Fig. 2b), contributing to a further water loss of 37.04 billion m³ (Fig. 2a). In conclusion, both direct human water consumption and ecological restoration have significantly contributed to TWS loss. However, the main causes of water loss differ by different sandy region, as detailed below.

Fig. 2: Changes of main driving factors and their contributions to TWS changes among sandy lands.
figure 2

a Contributions of each driving factor to water storage changes. b Changes of ET caused by ecological restoration over from 2002 to 2023.

Full size image

In the MuUs and Horqin Sandy Lands, agricultural water consumption is the primary driver of TWS decline, accounting for 0.53 and 1.66 times the TWS loss in each region, respectively (Fig. 2a). Over the past two decades, expanding irrigated areas (Supplementary Fig. 2) have significantly increased water consumption (Supplementary Fig. 1). In these arid and semi-arid regions with poor surface water and precipitation, groundwater pumping is a major agricultural water source, depleting groundwater and exacerbating TWS loss16. Moreover, vegetation water consumption from ecological restoration is the second-largest contributor to TWS loss, accounting for 0.29 and 0.83 times the net TWS loss in each region, respectively (Fig. 2a). Land use data show that vegetated areas expanded by 6107.25 km² in the MuUs and 8391.54 km² in the Horqin Sandy Lands, at the expense of cropland and bareland (Supplementary Fig. 3). This expansion is primarily attributed to the increase in vegetation (Supplementary Fig. 4) driven by ecological restoration projects4,17, which also contributed to higher ET(Fig. 2b). These findings contrast with previous studies15, which did not consider anthropogenic water use in the MuUs Sandy Land.

In the Otindag Sandy Land, we found that water loss was primarily driven by ecological restoration. Specifically, the increased ET from restoration accounted for 0.58 times the TWS loss, followed by agricultural water consumption, which accounted for 0.33 times the TWS loss (Fig. 2a). Significant forest expansion of 501.76 km² (20.76% of the revegetated area) occurred after the implementation of ecological restoration projects (Supplementary Fig. 3). The expanded forests, which consume 20–40% more water than grasslands, may increase deep-water use, accelerating groundwater depletion18,19. Our findings differ from Zhang et al. 14 due to two key reasons: (1) they did not consider vegetation water consumption and (2) they did not conduct a quantitative analysis, which may have led to biased conclusions. In contrast, this study fully accounts for vegetation water use and employs quantitative methods to assess water loss from each factor. And our results show that ecological restoration in the Otindag Sandy Land leads to greater vegetation evapotranspiration than anthropogenic water use.

Future TWS changes threatens sustainable development

Future TWS changes and potential water scarcity (2024–2060) depend on climate conditions and human water consumption under different scenario (Fig. 3).

Fig. 3: Projected changes in water storage and supply–demand balance under different emissions scenarios.
figure 3

a TWS anomalies (relative to the 2004–2009 time-mean value) during the historical (2002–2023) and future periods (2024–2060). b–d Projected changes in storage supply capacity and corresponding acceptable maximum growth of total water demand guaranteeing water security. The corresponding upper limits of human water use variability are specified. Notably, because GRACE estimated TWS anomalies rather than absolute TWS values, we used the relative changes of projected TWS anomaly relative to the historical anomaly to characterize future water storage changes. Hence, water gain (or water loss) means that the future TWS anomaly is higher (or lower) than the historical anomaly.

Full size image

Under the low emission scenario (SSP1-2.6), a slight water loss (0.48 mm, Fig. 3a) is observed only in the MuUs (0.5% reduction in water supply capacity, Fig. 3b). To avoid a worsening water supply–demand, human water consumption should decrease by 30% compared to historical levels (Fig. 3b). In contrast, the other three regions show modest water gains (4.45–6.76 mm, Fig. 3a), increasing water supply capacity by 5.70% to 8.71% (Fig. 3b), allowing slight increment in human water consumption (30–40%, Fig. 3b).

Under the medium emission scenario (SSP2-4.5), water loss will intensify. MuUs is expected to lose 6.34 mm TWS (6.54% reduction in water supply capacity, Fig. 3a, c). Accordingly, a 60% reduction in human water consumption is essential to prevent an imbalance in the water supply–demand relationship. Although a slight water gain in the other three Sandy Lands (2.76–5.68 mm, Fig. 3a) induced a limit increase in water supply capacity (4.30%–7.32%, Fig. 3c), rises in human water consumption may not exceed 15%–25% (Fig. 3c).

In the high emission scenario (SSP5-8.5), water loss will become more severe. All four regions of MuUs, Otindag, Horqin and Hulun Buir Sandy Land will experience water loss (18.49, 4.25, 1.70, and 0.29 mm, respectively, Fig. 3a), leading to reductions in water supply capacity of 19.07%, 5.44%, 2.19%, and 0.51% (Fig. 3d). To avoid a water crisis, anthropogenic water use must be strictly controlled, requiring reductions of 135%, 125%, 20%, and 20% (Fig. 3d).

For MuUs, the above reduction by 30%–135% in human water consumption would be highly challenging, especially in the background of potential population growth. Excessive water constraints would threaten food security in arid and semi-arid zones20 and reduce national economic outputs21. Therefore, in addition to tight control of local human water consumption, a feasible solution to this region is to strengthen water supply capacity through various ways, for example, water reuse and recycling22, rainwater harvesting and storage23, and improvement of hydrological connectivity between regions24. In the Otindag, Horqin, and Hulun Buir Sandy Lands, although small water gains are projected under SSP1-2.6 and SSP2-4.5, water scarcity will remain a risk unless human water use is controlled. It is also noted that these projections and budgets are based on the assumption that the vegetation coverage in the future remains unchanged at the current level. If the ecological restoration is further performed in the future in those regions, the human water use increment has to be squeezed out by the significant ecological water requirement.

Discussion

The ecological restoration in sandy regions is critical to local environmental protection and national ecological security. Remarkable success has been achieved in past decades, with a significant increase in vegetation cover5,25. However, as vegetation growth is sensitive to changes in water storage (Fig. 4a), the significantly reduced water resources may jeopardize vegetation that has been planted during ecological restoration26. Potential vegetation degradation may lead to desertification, worsening dust storm in northern China, threatening national ecological safety. Therefore, a comprehensive and persistent monitoring of water resources capacity is imperative during desertification control and ecological restoration. Additionally, strengthening local water supply capacity, such as water reuse and recycling22, rainwater harvesting and storage23, and improving hydrological connectivity between regions24,27,28,29 may be envisaged to guarantee the water supply and balance the groundwater deficit.

Fig. 4: Potential impacts of water storage on vegetation and agriculture.
figure 4

a Sensitivity of NDVI in the Four Major Sandy Lands to TWS. b Proportion of Agricultural and Related Industry Revenue in GDP from 2002 to 2023.

Full size image

We found that the large-scale artificial forests during ecological restoration processes in arid and semi-arid sandy areas have increased ET (Fig. 2c) and aggravated water shortages (Fig. 2a)30. In view of these, ecological restoration measures are suggested to be further optimized. For example, plant water-efficient tree species for revegetation to balance ecological restoration with water resource consumption, towards a healthy, stable, and sustainable ecosystem.

In addition to ecological significance, sandy regions are home for more than 13 million people who mainly live on high-water consuming agriculture31. In these regions, irrigated agriculture and related industries account for over 25% of total income (Fig. 4). As this sector heavily relies on water, reduced water availability may threaten both agricultural development and socio-economic sustainability16,32. Therefore, promoting efficient agricultural water use is crucial. For example, precision irrigation techniques (e.g., drip and sprinkler irrigation) can reduce water wastage and improve water efficiency33. Additionally, drought-resistant crops and native plants can reduce water consumption while maintaining crop yields. Besides, with the growing population in the future, it is possible for further reclamation of arable land. In view of low precipitation, scarce water resources, and large agricultural water requirements, water resource carrying capacity should be taken as a veto index for the land reserve investigation, avoiding serious water shortages and unsustainable agricultural production.

This study highlights that sustainable water management is crucial not only for China’s sandy regions but also for arid and semi-arid areas globally, including the western United States34, Mexico35, India36, and parts of Africa. In these regions, agricultural development and ecological restoration must be integrated with water conservation to avoid depletion, ecosystem degradation37, and food insecurity38. As water scarcity intensifies, these findings offer key lessons for mitigating future crises in water-stressed regions.

In summary, we found a decline in TWS anomalies over China’s sandy regions due to the decrease in groundwater storage attributed to human activities (including increased agricultural water consumption and ecological restoration) over the past two decades. In the next 40 years, the socioeconomic development of the study area could face severe water threats, requiring either a significant reduction or a limited increase in human water use to ensure water security. Given the projected water stress, it urges immediate scientific regulations on water resources, such as enhancing local water supply capacity through water reuse and rainwater harvesting, in addition to optimizing anthropogenic activities including agricultural irrigation and ecological restoration. This study provides a scientific basis for the sustainable management of water resources in sandy regions, offering valuable insights for achieving water-related Sustainable Development Goals.

Methods

Calculation of terrestrial water storage anomalies

The GRACE/GRACE-FO mass concentration (mascon) solutions, based on the latest versions of Release Number 06 (RL06) produced by the Jet Propulsion Lab, was used to estimate the monthly TWS anomaly during the historical period (2002–2023, Supplementary Table 1). All GRACE data reported here are anomalies relative to the 2004–2009 time-mean baseline39. This product is represented on a 0.5° × 0.5°grid, but the native resolution is 3 × 3 degree equal-area mascon block. Our study area spans more than two mascons (Supplementary Fig. 5). To capture more accurate signals of TWS, we define a 3° × 3°, 5° × 3°, 5° × 3°and 4° × 3.75°region for MuUs, Otindag, Horqin and Hulun Buir, respectively (Fig. 1a), which best covers our study area while maintaining the native resolution of a single mascon. We calculate the TWS time series of our defined region by averaging the covered mascons. We calculate the error of GRACE TWS by propagating that of JPL mascon solutions. Missing monthly TWS anomaly data during the GRACE period (April 2002–June 2017) caused by instrument failure and the gap period between the GRACE and GRACE-FO missions (July 2017–May 2018) were reconstructed employing the Random Forest method which is described in Section of Reconstruction and projection of TWS anomaly.

Calculation of lake water storage

We adapted the methodologies proposed by Khazaei40 and Mohamed41 to estimate the lake water storage in the study area from 2002 to 2023. First, we identified large lakes (over 100 km²) within the study area based on the HydroLAKES global lake dataset. It was found that there are only three large lakes, namely Dalinor (in Otindag Sandy Land), Hulun Lake (in Hulun Buir Sandy Land), and Buir Lake (Hulun Buir Sandy Land). Next, we obtained lake depth data using the Global Lakes Bathymetry Datasets (GLOBathy) and established the relationship between lake area, depth, and volume using the method of Khazaei40. Then, we extracted monthly lake area from 2002 to 2023 using the MNDWI index and Otsu’s method based on surface reflectance products derived from Landsat-5, Landsat-8, and Landsat-9 provided by the Google Earth Engine (GEE) platform. Finally, we estimated the lake water storage based on the established area-depth-volume relationship. To ensure consistency with the processing of GRACE-TWSA, we converted water storages to anomalies by removing the time series mean for January 2004-December 2009. It should be noted that the lake water storage anomalies, expressed in cubic meters, must be spatially distributed in order to estimate the water equivalent depth anomalies. For each sandy area, individual lake storage anomalies were summed and scaled by the sandy land area41.

The estimated results (see Supplementary Fig. 6) show that although the lakes water storage has significantly decreased, their contribution to the TWSA change is very small (see Supplementary Fig. 7). Therefore, lake water storage changes are not analyzed separately in our study.

Calculation of groundwater storage anomaly

The key components of TWS anomaly, namely soil moisture (summed over the 0–10, 10–40, 40–100, and 100–200 cm layers), canopy water equivalent, snow water equivalent, and runoff were obtained from GLDAS-NOAH during the period from 2002 to 2023. GLDAS-NOAH data were provided by NASA, with a spatial resolution of 0.25° × 0.25° and a temporal resolution of 1 month (Supplementary Table 1). Time series of each variable were obtained by averaging all grid over each sandy land. For each variable, anomalies were calculated by subtracting the 2004–2009 time-mean baseline to be consistent with GRACE/GRACE-FO data.

Groundwater storage anomaly was calculated based on the water mass balance equation (Eq. 1):

$${GWSA}={TWSA}-{CWEA}-{SWEA}-{SMSA}-{RA}-{LWSA}$$
(1)

where GWSA is the groundwater storage anomaly, TWSA is the terrestrial water storage anomaly (mm) obtained from GRACE/GRACE-FO, and CWEA, SWEA, SMSA, and RA represent canopy water equivalent anomaly, snow water equivalent anomaly, soil moisture storage anomaly, and runoff anomaly, respectively, as provided by GLDAS 2.1-NOAH. LWSA is the lake water storage anomaly (mm) estimated by the above method (See Section of Calculation of lake water storage). Given the considerably smaller magnitude of SWEA, CWEA, RA, and LWSA in comparison to TWSA (Supplementary Fig. 8), these were not analyzed in our study.

Trend analysis of time series datasets

The Theil-Sen Median trend analysis method (TS) and Mann-Kendall test (MK), which are widely used trend analysis methods for hydrological and meteorological variables, were used to analyze the long-term trends of water storage change. The TS method was used to estimate trends since it is robust and resistant to outliers42. As a common non-parametric test, the MK test was used to detect whether significant water storage changes have occurred over the study period.

Simulation of climate-driven TWSA

To quantify the impact of climate change on TWSA, we applied the method proposed by Humphrey et al. 43 to simulate the climate-driven TWSA from 2002 to 2023. This method has been widely used in studies assessing the impact of climate change on TWSA44,45,46. The basic formula is as follows (Eq. 2):

$${{TWS}}_{t}={{TWS}}_{t-1}\times {e}^{-\frac{1}{\tau \left(t\right)}}+P\left(t\right)$$
(2)

where t is a daily time vector, \({TWS}\) is daily TWS, \(P\) is the daily precipitation and \(\tau\) is the residence time of the water storage calculated based on daily temperature data. Further details on how this residence factor is calculated using temperature data and how the model is calibrated can be found in the original paper.

The meteorological data used to simulate climate-driven TWSA were derived from the fifth-generation ECMWF Reanalysis (ERA5) dataset, which includes daily temperature and precipitation records spanning from 1985 to 2023. This dataset is provided by the European Space Agency (ESA). Comprehensive details are available Supplementary Table 1.

Calculation of evapotranspiration caused by ecological restoration area

To quantify the impact of ecological restoration on TWS changes, we first identified areas of vegetation improvement resulting from ecological restoration and then assessed evapotranspiration (ET) in these regions.

We utilized land use/land cover (LULC) data for 2000 and 2018, provided by the Chinese Academy of Sciences’ Resource and Environmental Science Data Center, to delineate areas where ecological restoration led to vegetation improvement. The LULC dataset, with a spatial resolution of 30 m and a classification accuracy exceeding 94%47, includes 25 secondary land cover categories. Using the provided documentation47, the dataset was grouped into seven primary land cover categories: cropland, forest, dense grassland, sparse grassland, water, settlement, and bareland. A land-use transition matrix was then computed for the period 2000–2018. Based on ecological restoration policy information (Supplementary Table 2)4, regions where cropland and bareland were converted to forest and grassland (including dense-, and sparse grassland) or where sparse grasslands grassland was converted to dense grassland were classified as areas of vegetation improvement due to ecological restoration.

Additionally, ET during the growing seasons from 2002 to 2023 were obtained for the restoration areas using two data sources, namely MOD16A3 and the Global Land Evaporation Amsterdam Model (GLEAM). MOD16A3, a global ET dataset derived from MODIS satellite data, has a spatial resolution of 500 m and a temporal resolution of 8 days. GLEAM has a spatial resolution of 0.25° and a temporal resolution of 1 month. To ensure data quality, only “Good quality” pixels from the MOD16A3 and GLEAM datasets were retained according to quality control procedures. Annual ET was calculated by summing the growing season ET values. The mean of the two ET datasets was computed for each sandy area to represent the final ET.

Acquisition of human water consumption

Annual human water consumption data from 2002 to 2023, including agricultural, industrial, and domestic water usage, was used to represent the direct human impact on water resources. These data were obtained from the Water Resources Bulletins published by the Ministry of Water Resources the People’s Republic of China, which is widely used in examining supply–demand relationships27,48,49. As a water-scarce country, China has implemented a rigorous guideline for measuring human water withdrawals to ensure data reliability (http://www.mwr.gov.cn/zwgk/gknr/202004/t20200426_1441177.html). Specifically, (1) surface water withdrawals from rivers, lakes, reservoirs, and ponds are measured at the intake points; (2) groundwater abstraction from wells is measured at the wellheads; (3) water sourced from rainwater, seawater desalination, and recycled water is measured at the respective intake points; and (4) withdrawals from the water supply network are calculated based on the amount of water distributed to consumers. We converted these data to water equivalent thickness, dividing the amount of water by the area of each sandy land.

Quantitative attributions of terrestrial water storage loss

Climate change and human activities contributed to the water resources changes of China’s Four major Sandy Regions. Analysis based on ERA5 dataset found that the climate over sandy land regions has undergone drying and warming trends (blue line and red bar in Supplementary Fig. 4). Moreover, the increased human water consumption (Supplementary Fig. 1) and extensive vegetation greenness (indicated by NDVI derived from multi satellite data, green line in Supplementary Fig. 4) caused by ecological restoration may also affect water resources through evapotranspiration50. Consequently, the GRACE-TWSA should reflect the combined effects of climate-driven TWSA variation, human water consumption, and evapotranspiration due to ecological restoration, as represented in the following Eq. (3):

$${{TWSC}}_{{GRACE}}={{TWSC}}_{{climate}}+{HWC}+{ETC}$$
(3)

where \({{TWSC}}_{{GRACE}}\) represents the TWS changes observed by GRACE from 2002 to 2023 (Fig. 1b). \({{TWSC}}_{{climate}}\) refers to the climate-driven TWS changes calculated using the method of Humphrey et al. (see Section on Simulation of climate-driven TWSA). ETC represents changes in evapotranspiration over ecological restoration areas. HWC represents human water consumption over the same period.

Then, we quantified the contributions of climate change, human water consumption, and evapotranspiration due to ecological restoration to the changes in TWSA, as shown in the following formulas:

$${\theta }_{{climate}}=\frac{{{TWSC}}_{{climate}}}{{{TWSC}}_{{GRACE}}}\times 100 \%$$
(4)
$${\theta }_{{human}}=\frac{{HWC}}{{{TWSC}}_{{GRACE}}}\times 100 \%$$
(5)
$${\theta }_{{restoration}}=\frac{{ETC}}{{{TWSC}}_{{GRACE}}}\times 100 \%$$
(6)

where \(\theta\) represents the contributions of driving factors.

Calculation of vegetation greenness

Two NDVI datasets for the growing season (May–September) from 1985 to 2023 were used to indicate vegetation greenness, namely the latest Global Inventory Modeling and Mapping Studies 3rd generation (GIMMS-3g) and the Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 MOD13A3. GIMMS-3g NDVI was retrieved from the National Oceanic and Atmospheric Administration’s Advanced Very High-Resolution Radiometer (AVHRR) data available bi-monthly from 1985 to 2015 with a spatial resolution of 0.083°. This NDVI product has been preprocessed for radiation correction, geometric correction, and image enhancement. MOD13A3 NDVI is available on a monthly basis from February 2000 to the present with a spatial resolution of 1 km. To ensure the quality and consistency between the datasets, we have implemented the following data processing steps: (1) Only high-quality pixels with pixel reliability of ‘good data’ or a value of ‘0010’ according to the user manual were used. (2) The maximum value composite method was applied to create monthly GIMMS-3g NDVI time series. (3) Linear interpolation was performed for missing pixel values by using nearby time points in both MODIS and GIMMS-3g datasets. (4) Both NDVI datasets were masked using the Land Use and Land Cover (LULC) data to eliminate the impact of crops on vegetation greenness. (5) The GIMMS-3g NDVI was calibrated to MODIS NDVI using a linear regression between their monthly values from 2002 to 2015 to obtain a consistent NDVI record over the entire study period (Supplementary Fig. 9)51.

Reconstruction and projection of TWS anomaly

We employed Random Forest (RF), a machine-learning-based approach, to reconstruct missing TWS anomaly for each sandy land and extend the predictions until 2060. RF is a non-parametric ensemble learning algorithm that builds regression trees from bootstrapped samples of the training data. Then each node is further split based on randomly selected predictor variables from the full array until the node’s size reaches the minimum value. Therefore. three of its main hyperparameters, including the number of trees (ntree), the number of predictor variables (mtry), and the minimum node size (mnode) are required in the RF model. Forcing data for RF models included bias-corrected precipitation, air temperature, and surface short-wave radiation derived from five Global Climate Models (GCMs) in Coupled Model Intercomparison Project Phase 6 (CMIP6) and simulated TWS from the latest CMIP6 CESM2-WACCM model under the r1i1p1f1 format52,53. Details of the GCMs dataset are shown in Supplementary Table 1. In order to avoid overfitting caused by relatively small number of target data (GRACE and GRACE-FO observations of 228 months in total; April 2002–December 2023), we strictly set small number of trees (from 1 to 150) and nodes (from 1 to 5), which could effectively reduce the risk of overfitting from the model structure52. For each sandy land, the forcing and target data from April 2002–December 2022 were randomly divided into training and testing sets, accounting for 70% and 30% of the samples (216 months total), respectively. The remaining data (January 2023-December 2023) were used to further validation of TWSA projection. A detailed workflow of model training and testing is presented in Supplementary Fig. 10. Model performance was evaluated using root mean square Error (RMSE), correlation coefficient (CC), and Nash-Sutcliffe efficiency coefficient (NSE).

Supplementary Fig. 11 shows the RF model performance based on the testing sets from different GCMs. The results indicate that the RF model successfully reproduced TWSA with RMSE ranging from 0.63 to 2.91 mm, CC values from 0.73 to 0.96, and NSE between 0.65 and 0.96. Furthermore, the independent validation based on the latest GRACE-FO observations during January 2023‒December 2023 also confirms the good performance of the model (Supplementary Table 3). Time series analysis of the TWSA changes derived from both GRACE observations and the ensemble mean of model simulations (Supplementary Fig. 12) shows consistent trends across all sandy land regions, further demonstrating the model’s reliability in simulating TWSA. Given the strong agreement between the ensemble mean of multi-model simulations and GRACE observations, we used the ensemble mean to fill in missing GRACE data and to project future TWS anomalies through 2060.

Correction of forcing data derived from CMIP 6 models

The climate forcing data of the RF models, including temperature, precipitation, and surface short-wave radiation, were derived from five GCMs with a nominal resolution of 100 km at daily timescales. These data were obtained from CMIP6 Scenario Model Intercomparison Project covering both historical (2002–2014) and future (2015–2060) periods under three different radiative forcing scenarios, namely low-forcing scenarios (SSP1-2.6), moderate-forcing scenarios (SSP2-4.5), and high-forcing scenarios (SSP5-8.5). Details are described in Supplementary Table 1.

Due to inherent biases in the GCMs, a bias correction was applied to the datasets using the quantile mapping (QM) method, with ERA5 data from the European Centre for Medium-Range Weather Forecasts (ECMWF) serving as the reference. The correction was performed on daily-averaged data, which were then aggregated (or averaged) to monthly timescales. The historical and prediction periods were set as 2002–2023 and 2024–2060, respectively.

The QM method involves three steps54: (1) defining the cumulative distribution functions (CDFs) for both the reference and model data during the historical period, (2) calculating the model data quantiles based on the CDFs, and (3) using these quantiles to generate the bias-corrected data.

Quantification of water demand and supply capacity

Changes in water storage impact water security directly. Therefore, our study assessed future water security based on TWS anomaly and total water use. Here, storage supply capacity corresponds to TWS anomaly, while total water demand is the total water use, which is the sum of human water use and ecological water requirement. Because GRACE estimated TWS anomalies rather than absolute TWS values, and total water use is an absolute value, a direct comparison between these two variables is challenging. Hence, we assumed a tight balance in the water supply and demand from 2002 to 2023 across the study area, which describes a state where the water supply and demand are nearly balanced, but water resources are scarce. Given this assumption, we focused on the change in future (2024–2060) water supply capacity and water demand relative to historical (2002–2023) water demand52. The percentage change in storage supply capacity and TWD was calculated as:

$${{SSC}}_{C}=\quad \frac{{\overline {{SSC}_{2022 \sim 2060}}}-{\overline {{SSC}_{2002 \sim 2023}}}}{{\overline {{TWD}_{2002 \sim 2023}}}}$$
(7)
$${{TWD}}_{C}=\quad \frac{{\overline {{TWD}_{2022 \sim 2060}}}-{\overline {{TWD}_{2002 \sim 2023}}}}{{\overline {{TWD}_{2002 \sim 2023}}}}$$
(8)

where SSCC and TWDC represent the percentage changes of future storage supply capacity and total water demand relative to historical (2002–2023) total water demand, respectively. When TWDC is higher than SSCC, water security is threatened, and vice versa. TWD is calculated using Eq. (9):

$${TWD}={EWR}+{HWC}$$
(9)

where EWR and HWC represent ecological water requirement and human water consumption, respectively. EWR is estimated using Eq. (10) and eq (11)55:

$${ERW}= {k}_{c}\times {PET}$$
(10)
$${k}_{c}=\quad -0.0012\times {t}^{3}+0.0265\times {t}^{2}-0.0825\times t+0.2166$$
(11)

where PET is potential evapotranspiration calculated using the method proposed by Thonthwaite56; t represents the number of ten-day periods of vegetation growth, which is set to 1555. The historical HWC was provided by the water resources bulletin. Due to the absence of accurate future HWC values, we assumed that future HWC will increase (0–100%) or decrease (−100% to 0%) by a percentage of historical HWC.

Calculation of vegetation sensibility to TWSA

Vegetation sensitivity to TWSA was quantified using an autoregressive model that integrated the normalized difference vegetation index (NDVI), TWS, and meteorological data for the growing season (May–September) from 2002 to 2023. Following the methodology of De Keersmaecker et al. 57. TWSA was introduced into the autoregressive model as follows:

$${{NDVI}}_{t}=\alpha \times {{{NDVI}}}_{t-1}+\beta \times {{TWSA}}_{t}+\gamma \times {{prec}}_{t}+\delta \times {{temp}}_{t}+\varepsilon$$
(12)

where \({{NDVI}}_{t}\) and \({{NDVI}}_{t-1}\) represent the \({NDVI}\) at times t and t − 1, respectively. \({{TWSA}}_{t}\), \({{prec}}_{t}\) and \({{temp}}_{t}\) correspond to the total water storage anomaly, precipitation, and temperature at time t. The parameter \(\beta\) reflects the sensitivity of vegetation to TWSA, with a larger absolute value indicating greater sensitivity.

NDVI data used to calculate vegetation sensitivity to TWSA were obtained from MOD13A3. Comprehensive description and the preprocessing steps were detailed in Section of Calculation of vegetation greeness. Meteorological data were from the ERA5 dataset, with full details provided in Supplementary Table 1. Daily meteorological variables were aggregated or averaged to produce monthly timescale data. Before calculating sensitivity, the monthly NDVI, TWSA, and meteorological data were standardized using Eq. (13) to remove the climatological signals from the original time series.

$${z}_{i,j}=\frac{{x}_{i,j}-{{\bar{x}}_{j}}}{{\sigma }_{j}}$$
(13)

where i represents the year (2002–2023), j denotes the month (May–September), \({\overline{x}_{j}}\) and \({\sigma }_{j}\) are the mean and standard deviation of time series x for month j58.

Limitations and uncertainties

Our study is subject to several limitations and uncertainties related to observational data, attribution analysis, and model predictions. Specifically, the low spatial resolution of the GRACE satellite data introduces potential uncertainties, particularly when monitoring water resources at finer spatial scales59. In the attribution analysis, we analyzed the causes of water loss using a simplified approach, as detailed in the “Quantitative Attributions of Terrestrial Water Storage Loss” section. However, this method does not adequately account for the interaction effects between different drivers. For instance, we adopted ET from ecological restoration areas to represent the water loss due to ecological restoration. However, this approach may also reflect changes in ET driven by climate factors, introducing uncertainty in isolating the specific impact of restoration. Moreover, increased vegetation cover resulting from ecological restoration can alter local and regional hydrological cycles. These interactions were not captured in our current methodology. Future research should focus on the interactions between various drivers of water loss and provide more precise attribution analyses. For future TWS projections, uncertainties arise from the meteorological input data, which may be influenced by model biases, the absence of explicit human activity simulations, and the limited availability of data60,61. As a result, discrepancies remain between the machine learning model outputs and satellite observations. Since water supply capacity changes were estimated based on projected TWS, uncertainties in meteorological data under different scenarios will further contribute to the overall uncertainty in water supply capacity projections.