Introduction

The Tibetan Plateau (TP) is known as the Asian water tower1,2. Major rivers, including the Indus, Ganges, Brahmaputra, Mekong, Salween, Yellow, and Yangtze, originate from the Asian water tower and provide freshwater resources for nearly 2 billion people3. The freshwater storage (e.g., groundwater storage (GWS), snow) on the TP and its surrounding mountain ranges is instrumental for ecosystem restoration, agriculture, livestock breeding, and downstream environmental sustainability4. As a critical component of the hydrologic cycle5,6, groundwater can provide baseflow to rivers and lakes7 and supply water to root-zone soils during dry periods, thus maintaining ecosystems8. Understanding changes in GWS is critical for water resources management, engineering design (e.g., dams and tunnels), and geothermal energy development9.

Groundwater movement on the TP is typically driven by the topographic gradient, recharged by precipitation and glacier meltwater, and discharged in valleys4,10. Accelerated warming and changes in precipitation regimes on the TP have led to glacier retreat, snowmelt, and permafrost degradation, affecting storage and discharge processes and thereby altering GWS dynamics2,11,12. In general, with little human disturbance, GWS changes in TP are primarily related to changes in climatic factors (changes in precipitation, temperature, and shortwave radiation resulting in altered groundwater recharge and discharge). Owing to the challenging climatic conditions of the TP, groundwater wells, meteorological, and hydrological observation stations are sparsely distributed. The lack of direct access to GWS from satellite observations and the simplification of groundwater modeling in large-scale hydrological models5 have posed substantial obstacles to the accurate assessment of GWS change. Most previous studies have predominantly focused on permafrost degradation, glacier shrinkage, lake expansion, and changes in terrestrial water storage (TWS) on the TP12,13,14,15. Little is known regarding GWS on the TP and its response to climate change. Previous studies estimated the GWS changes on the TP using the Gravity Recovery and Climate Experiment (GRACE) satellite data and removing the known hydrological components (such as soil moisture and snow water) based on the water mass budget16,17. Recent studies4,10 have estimated GWS changes on the TP and removed the effects of glacier ice melt and lake water storage changes, but precise separation of soil water storage (Fig. S1) from TWS remained elusive. In addition, GRACE products derived from different solutions vary widely on the TP12,18 (Fig. S2), and there is no consensus on the estimation range of GWS from these products for both the entire TP and the individual basin.

Furthermore, future changes in GWS on the TP in a warming climate are still uncertain. Most past studies have focused on the future evolution of the TWS, glacier mass and lake water storage12,19,20,21. The changing trend of TWS on the TP is expected to slow down and possibly reach a “new equilibrium” (balance between solid and liquid water resources)12, the potential net loss of TWS on the TP could be as high as 230 billion m3. Projections indicate a deceleration in the increase of lake water storage on the Inner TP, but continued lake expansion may lead to drainage reorganization and an elevated risk of hydrogeological hazards21. High mountainous regions in Asia are anticipated to lose over 50% of glacier mass by the end of this century compared to 2000, with one-fifth of small glaciers likely to disappear by mid-century19. These findings offer valuable insights for decision-making in the sustainable development of water resources on the TP. However, potential future changes in GWS on the TP remain unclear. Climate system model-driven hydrologic modeling, although providing an opportunity to predict hydrologic processes, faces challenges due to inherent uncertainties22. Machine learning, capable of constraining predictions from observations, emerges as an effective approach for predicting hydrologic variables12,20,23. Therefore, a comprehensive understanding and reliable benchmarks of the magnitude of GWS change and its response to future climate change based on machine learning are essential.

This study endeavored to bridge these knowledge gaps by providing estimates of historical (2002–2018) GWS changes and projections into the future (2019–2100) on the TP (Fig. 1). Based on the water mass budget and a deep learning approach, we employed high-resolution land surface modeling, satellite observations, and global climate model data to provide the historical and future GWS changes under different Shared Socio-economic Pathways (SSPs) and Representative Concentration Pathway (RCP) scenarios (including SSP1-RCP2.6, SSP2-RCP4.5, and SSP5-RCP8.5, abbreviated as SSP126, SSP245, and SSP585, respectively). Our investigation encompassed an exploration of the contributions of different climatic factors to GWS changes, identifying the profound impact of climate change on GWS in the TP.

Fig. 1: Trends in groundwater storage anomalies (GWSA) on the TP from 2002 to 2018.
figure 1

a Spatial distribution of trends in GWSA, as estimated from multi-source observations and high-resolution land surface model simulation. Basin boundaries (red lines), lakes, and glaciers distributions are indicated. Triangles represent trends in GWSA for each basin (red and blue triangles represent decreasing and increasing GWSA trends, respectively; only basins with trends > 0.3 Gt yr−1 are shown). The size of the triangles represents the magnitude of GWSA trends (Gt yr−1). The boundary of TP is derived from SRTM DEM at an altitude of 2500 m. b Annual GWSA on the TP for each basin during 2002–2018. The black line represents the GWSA estimated based on monthly lake water storage data (50 lakes larger than 150 km2 and two lakes larger than 100 km2, see “Methods”), and blue dots represent data from 5-year intervals for 1132 lakes larger than 1 km2 (see “Methods”).

Full size image

Results

GWS evolution over the past two decades

Based on the water mass budget, compiled multi-source satellite data, and land surface modeling (Table S1, and Figs. S3–S6), we derived the GWS anomaly (GWSA) for the entire TP and each sub-basin from 2002 to 2018 (Fig. 1). Comparisons with observations illustrate the reasonableness of our estimated GWSA (Fig. S7). The endorheic region experienced a decreasing trend in GWS (−0.12 Gt yr−1, GWS changes in the endorheic and exorheic regions were apportioned by multiplying the percentage of its contribution to the TP GWS changes), especially in the Inner TP and Tarim river basin. In contrast, GWS increased in the exorheic region (3.63 Gt yr−1), with the most increase in the Yangtze river basin (Fig. 1). The spatial heterogeneity of GWS is closely related to glacier and lake patterns (Figs. S8–S10), with GWS increasing in areas of glacier degradation, decreasing in the Tarim river basin where glacier mass is increasing, and decreasing in the Inner TP. The main reason for this pattern correlation may be that widespread active normal faults (Fig. S11) are conducive to groundwater discharge to downstream lakes in the Inner TP and reduces GWS24. In addition, the rift valley complex (20 km wide and 300 km long) in the southern Qiangtang basin is likely to be an essential transport channel for groundwater seepage, contributing to the GWS reduction in the Inner TP25.

All GWS estimations from different combinations of GRACE and lake water products exhibited consistent trends across the entire TP, except for slight differences at the basin scale (Table S2). Only minor discrepancies in the magnitude of the GWS trends could be found for different lake water level observations, particularly in the Tarim and Indus river basins (Table S2 and Fig. 1). Our results showed a significant increase in GWS (3.51 ± 2.40 Gt yr−1) on the TP during 2002–2018 (Fig. 1 and Fig. S8). One of the main reasons was that there has been a significant increase in annual precipitation since the 1980s (Fig. S12). In particular, annual precipitation increased quickly between 2015 and 2018, which led to a marked increase (14.5 Gt yr−1) in GWS. Another reason was that glacier mass loss and permafrost degradation on the TP due to global warming increased surface recharge to GWS2.

As shown in Fig. 1, GWS changes displayed divergent variations in different regions of the TP (Fig. 1). The Yangtze, Qaidam, Amu, Indus, and Brahmaputra river basins showed a substantial GWS increase (0.33–1.96 Gt yr−1). In contrast, a significant reduction of GWS was found in the Tarim and Inner TP (−0.99 ± 0.14 and −1.16 ± 0.42Gt yr−1). The most significant increase in GWS was in the Yangtze river basin (1.96 ± 0.48 Gt yr−1), possibly attributed to abundant groundwater recharge from precipitation (Fig. S12). GWS estimations for the Yellow river basin derived from different GRACE products exhibited slight differences (Table S2). Center for Space Research mascon (CSR-M) showed decreasing GWS in the Yellow river basin, whereas other products had increasing trends. Owing to the dominant role of increasing precipitation in controlling lake water, most lakes on the Inner TP expanded rapidly during 2000–201626 (Figs. S10, S12), but GWS was found to significantly decrease based on the water mass budget (Fig. 1). This may be related to groundwater discharge to the downstream lakes under the influence of active faults in the Inner TP24 and groundwater seepage due to the rift valley in the southern Qiangtang basin25. In 2012–2014, earthquakes occurred near the north-south rift valley in the southern Inner TP (Fig. S13), which may have been the main reason for the significant decline around 2012 (Fig. 1).

Our GWS estimates closely align with those from previous studies4,10,27,28,29,30,31 in TP and its subregions (Inner TP, Qaidam, Yellow, and Yangtze river basins, Fig. 2), where the range of GWS largely overlaps. For the early 21st century (2003–2009), the GWS in the Inner TP showed an increasing trend (~0.6 Gt yr−1), consistent with findings in previous studies31 (5.01 ± 1.59 Gt yr−1). Discrepancies in magnitude may be related to differences in GRACE products. This previous study31, utilizing the multi-basin inversion method, reported a TWSA trend of 11.06 ± 1.2 Gt yr−1, whereas the GRACE averaged value used in our study indicated only 5.0 Gt yr−1. In addition, discrepancies in GWS estimates for the Inner TP resulted from different lake datasets and data processing method. For example, GWS data derived from 1132 lake annual water volume change data were processed as anomalies relative to the time-averaged baseline for the 2005–2010 period owing to the 5-year interval of available data. Different GRACE products have been shown to exhibit disparities in the magnitude and trend of TWS in the TP. Thus, our approach of estimating GWS using the average value of multiple GRACE products contributes to reducing uncertainty. Our results are consistent with a recent study10, primarily due to the utilization of multiple GRACE products for estimation. For the Yangtze and Yellow river basins, the differences between the results are mainly related to the extent of the study area. Overall, the results of this study are based on multiple GRACE products, state-of-the-art atmospheric data, and high-resolution land surface simulations, and comparisons with observations and previous studies show the accuracy and reasonableness of the estimates in this study.

Fig. 2: Estimations of GWS change rate from previous studies4,10,27,28,29,30,31 (colored) and this study (black).
figure 2

Superscripts 1 to 7 represent the GWS estimations from Zou et al.10, Xiang et al.4, Zhang et al.31, Qiao et al.30, Bi et al.27, Jiao et al.29, and Chao et al.28. For Inner TP, results of Xiang et al4. are for the central part of Qiangtang Nature Reserve. For Qaidam, results from Xiang et al4. are for the eastern part of the Qaidam basin, which covers about 68% of the Qaidam basin. For Yellow river basin, results from Xiang et al.4 are for the Yellow river source basin. For Yangtze river basin, 2* and 2** denote results for Yangtze river source basin and Jinsha river basin from Xiang et al.4, and results from Chao et al.28 are for Jinsha river basin.

Full size image

Response of GWS to climate change

To investigate how climate change affected GWS changes on the TP, three climatic variables (air temperature, precipitation, and solar radiation) were chosen to investigate their effects on GWS changes over each basin on the TP (Figs. S14–S15). Although the northern part of the Yangtze river basin showed an increasing trend in annual precipitation (3.1 mm yr−1), increasing temperature (0.02 °C yr−1), insignificant runoff increase (0.4 mm yr−1) and significant evapotranspiration increase (2.0 mm yr−1) coupled with a groundwater lateral drainage affected by topography led to a slight decrease in GWS across the northern part of the Yangtze river basin (−0.49 mm yr−1, Fig. S8). In the glacier-distributed areas (west of the Yangtze river basin), the significantly increasing shortwave radiation and temperature led to accelerated melting of glaciers, favoring groundwater recharge (Fig. S14). Consequently, GWS in the Yangtze river basin exhibited a significant increase (1.96 ± 0.48 Gt yr−1) related to increased precipitation, melting of soil ice (−1.28 Gt yr−1) due to higher temperatures, and glacier ablation (-0.79 Gt yr−1) caused by increased radiation. The results of the regression analysis (See “Methods” for details) showed that the relative contributions of precipitation, temperature, and shortwave radiation to the changes in GWS across the Yangtze river basin were comparable. Therefore, the GWS in the Yangtze river basin increased under the combined effects of increased precipitation and warming-induced glacial recharge.

Meltwater is the primary recharge source of surface runoff for endorheic basins in arid and semi-arid regions2. Shortwave radiation decreased at a rate of 9 W m−2 yr−1 in the Tarim river basin (Fig. S15), and there was a decrease in winter temperatures (Fig. S14) and an increase in glacier mass (Fig. S9) in the Karakorum. This may explain the decrease in GWS in the Tarim river basin. In the Qaidam basin, although the glacier area is relatively small, accounting for approximately 0.5% of the area in the basin compared to 8% in the Tarim river basin (Table S2), the region experienced a decrease in shortwave radiation (−14 W m−2 yr−1) and an increase in precipitation (1.55 mm yr−1), resulting in an increase in GWS. The regression analysis results also indicated that precipitation and shortwave radiation caused changes in GWS in the Qaidam basin (Fig. 3 and Figs. S16, S17). The increased winter temperatures and shortwave radiation (Figs. S14, S15) led to glacier loss (Fig. S9) and increased GWS in the Amu Darya river basin. The results of the regression analysis also indicated that temperature was the main factor causing changes in GWS in the Amu Darya river basin (Fig. S17). In the Inner TP, the spatial distribution of GWS changes (Fig. S8) is roughly consistent with that of precipitation (Fig. S14). However, active normal faults are distributed in most of the lake basins on the TP (Fig. S11), which are more favorable for groundwater discharge to the downstream lakes and cross-basin groundwater transport24; accompanying frequent earthquakes occur along the faults (Fig. S11) may cause significant leakage (Fig. S11), which may explain the decrease in GWS in the Inner TP. In particular, in the western endorheic region, the large amount of groundwater recharge substantially changed the seasonal and interannual lake water level24. For the glacier-dominated exorheic region, increased precipitation in the Indus river basin, and increased winter temperatures and shortwave radiation, resulted in a significant decrease in glacier mass. In particular, the Indus and Amu Darya river basins had similar glacier coverage (8% and 7.5%, respectively). Although annual precipitation was more abundant in the Indus river basin, the magnitude of the GWS increase was comparable in the two basins (0.76 Gt yr−1 and 0.73 Gt yr−1, respectively), which may be related to groundwater lateral drainage in the Indus river basin.

Fig. 3: Contributions of climate factors (precipitation, temperature, and shortwave radiation) to GWS changes on the TP from 2002 to 2018.
figure 3

Climate data sourced from a high accuracy atmospheric dataset48. Contributions were calculated for each basin based on a regression method (see “Methods”).

Full size image

Overall, precipitation, temperature, and shortwave radiation collectively affected the change of GWS over the whole TP, and an increase in GWS in the exorheic region is associated with warming-induced glacial melting. The regression analysis results showed that in approximately 70% of the basins with GWS changes of more than 0.3 Gt yr−1 (Fig. S17), these three climatic factors (precipitation, temperature, and shortwave radiation) had nearly equal contributions to the GWS changes.

Future projection of GWS

We projected the future GWS using a deep learning approach, climate model outputs, and bias correction methods (see “Methods”, Figs. S18–S22). The GWS in the TP will persistently increase for different future climate scenarios, but the growth rate is slowing down (Fig. 4, 0.14 Gt yr−1 for 2079–2100 under SSP585). This slowing increase trend in GWS under the high emissions scenario at the end of the 21st century (2079–2100) may be related to the pronounced warming (0.09 °C yr−1, almost twice as in 2002–2018) and unsustainable glacier meltwater32, which reduces groundwater recharge, even the precipitation increase rate is comparable to historical periods. A slight increase in GWS in the TP under the low emissions scenario between 2019 and 2040 is associated with slightly increased precipitation (1.34 mm yr−1, Fig. S19). Although the GWS increase rate for the entire TP is higher in the far future than in the near future in the SSP245 scenario, the increase rate slows down relative to the historical period.

Fig. 4: Projected GWS trends under different climate scenarios on the TP.
figure 4

Bars represent projected GWS trends in the near (2019–2040) and far future (2079–2100) under three SSP scenarios. Projections were derived from a deep learning approach using bias-corrected climate projections from CMIP6 model outputs.

Full size image

Under all three forcing scenarios (SSP126, SSP245, and SSP585), GWS projections in endorheic basins of the TP showed substantially increasing trends both in the near (2019–2040) and far future (2079–2100) (Fig. 4). The Amu Darya river basin will likely experience a sustained loss of GWS, with the fastest rate of loss in 2019–2040 and a weakening rate of loss in 2079–2100. The projected GWS trend in the Indus river basin over the next decade (2019–2040) under the SSP585 scenario (0.92 Gt yr−1) is comparable to the 2002–2018 trend (0.73 Gt yr−1). However, the subsequent decades show little change, particularly under the SSP245 scenario, with a decreasing trend of less than 0.1 Gt yr−1. The projected GWS trends for the Inner TP indicate a slight increase in GWS over the next decade, with GWS increasing more in the Inner TP under the SSP126 scenario compared to the other two scenarios owing to insignificant warming (Fig. 4 and Figs. S23–S25). Subsequently, the rate of increase gradually slowed under the SSP245 and SSP585 scenarios, and there will be a small decrease under the SSP126 scenario, particularly towards the end of the century. The projected GWS change was slower (80% of basins had a trend of less than 0.4 Gt yr−1) relative to the historical period (2002–2018) on the TP, which is consistent with the conclusion of the separate projections of TWS12, indicating that the TP is approaching a new equilibrium. For basins with a high proportion of glaciers (such as Amu Darya river basin, Indus river basin, and Tarim river basin), the GWS projections basically show a more minor change rate in the far future than the near future, which is also consistent with previous predictions of unsustainable glacier meltwater32. The GWS in the Tarim Basin, one of the basins with a high proportion of glaciers, is projected to increase in the near future, which is probably mainly due to the buffering effect of glacial increases in the historical period.

Although the projected GWS over the whole TP shows a slow increase trend, the GWS loss in individual basins was still not negligible. The projected groundwater loss in the Amu Darya river basin was approximately 3.4–4.6 Gt for 2019–2040, and the total groundwater loss in the Yellow and Yangtze river basins was approximately 15.5 Gt for 2079–2100 for the SSP 245 scenario. Although the projected GWS trend on the TP was an increase of 2.4–3.3 Gt, the glacier-dominated Amu Darya river basin and exorheic basins (Yangtze and Yellow river basins) may face threats to groundwater resources.

Changes in climatic factors, including precipitation, temperature, and shortwave radiation, are consistent with the projected changes in GWS. In the Ganges and Brahmaputra river basins, where precipitation increases significantly (Figs. S23–S25), projected GWS increases in the near future for all three scenarios. The Amu Darya river basin had a slight increase in winter precipitation, but the region is likely more sensitive to increases in temperature and shortwave radiation. The rate of decrease in GWS in the Amu Darya river basin in 2019–2100 under the SSP585 scenario was approximately six times that under the SSP126 scenario and twice that under the SSP245 scenario. A relatively higher increase in precipitation in the Indus river basin may account for the slight increase in GWS. The increasing temperature trend in the Inner TP is smaller in the SSP126 scenario than in the SSP245 and SSP585 scenarios. Consequently, the projected increase in GWS was greater in the low-emission scenario than in the high-emission scenario. Overall, GWS losses are less under the low emission scenario, suggesting that climate change mitigation efforts may improve groundwater availability on the TP in the future.

Discussion

This study investigated historical changes and future projections of GWS on the TP and their response to climate change, providing insight into alpine groundwater resources that were important for downstream water supply. Increased GWS on the TP (3.51 ± 2.40 Gt yr−1) during the historical period, mainly in the exorheic region, was related to warming-induced glacial meltwater. Divergent implications of climate change result in the spatial heterogeneity of GWS change on the TP. The reduction of melting glaciers caused by the decrease in shortwave radiation in the Tarim river basin led to a reduction in GWS. Warmer temperatures have led to glacier reduction, dominating the increase in GWS in the Amu Darya river basin. Widespread active normal faults are conducive to groundwater discharge to downstream lakes and rift valley complex-induced groundwater seepage in the last two decades have resulted in a significant decrease in GWS in the lake-dense Inner TP (−1.16 ± 0.42 Gt yr−1). The most pronounced increase in GWS in the Yangtze river basin (1.96 ± 0.48 Gt yr−1) was related to the combined effects of increased precipitation and warming-induced soil ice and glacial recharge. Based on future GWS projections under the three different scenarios (SSP126, SSP245, and SSP585), the GWS in the TP will be persistently increasing for different future climate scenarios but the growth rate is slowing down. The Amu Darya, Yangtze, and Yellow river basins are likely to experience GWS losses.

Future reductions in GWS can have significant implications for water supply in densely populated areas downstream, complicating water resources management policies. Predictions suggest that the Amu Darya river basin may experience a decrease in GWS over the 21st century (Fig. 4). Glaciers and snowmelt in the upper Amu Darya river basin are important water sources for the rivers and alluvial aquifers33. Thirty percent of the Amu Darya river basin’s aquifers are transboundary, with the majority being extensively exploited33. As early as 1999, Uzbekistan, downstream of the Amu Darya river basin, utilized about 99% of the authorized GWS for drinking and irrigation33. The over-exploitation of groundwater not only negatively affects localized eco-hydrological processes but can also exacerbate tensions and potential conflicts. For example, the shared aquifer between Uzbekistan and Tajikistan has been a major topic of dispute due to the lack of proper groundwater accounting33. The stable water source in the Amu Darya river basin is disappearing as warming intensifies glacial retreat19, and with the capacity of the reservoirs already exceeding precipitation34, adapting options to future changes in water storage are very limited. Therefore, effective water conservation measures and water diversion projects in the downstream countries of the Amu Darya river basin are necessary to sustain groundwater resources and maintain national water security.

In addition, the projection results also indicate that the GWS in the Yangtze and Yellow river basins will decrease significantly by the end of this century (2079–2100). Groundwater is the primary source of domestic water and industrial production in many rural areas in the Yangtze river basin35. Excessive groundwater depletion during the historical period (2003–2009) triggered groundwater droughts36. Although reservoir leakage can recharge groundwater, the storage capacity of reservoirs hardly offsets the negative effects of groundwater extraction37. Groundwater extraction in the Yellow river basin has remained relatively high from 2003 to 2019, roughly three times the annual extraction in the Yangtze river basin. Agricultural irrigation is the main source of groundwater consumption in the Yellow river basin. Despite implemented water conservation measures, the declining trend in GWS has not been reversed38. In addition, due to the rapid economic development of the Yellow river basin, industrial and domestic water consumption has significantly increased in recent years39. The continuous population increase and economic development were projected to raise water demand40, posing a challenge for groundwater resource management in the Yangtze and Yellow river basins, where achieving a reasonable trade-off between water resources, food production, and industrial output value becomes increasingly challenging.

Methods

GWS estimations for historical period

Based on the water mass budget31, TWS is the sum of surface (including lake water storage, snow water, and glacier-specific mass changes) and subsurface water storage components (including soil water and groundwater). Therefore, the GWS estimation can be written as

$${GWSA}={TWSA}-P{WSA}-{SWA}-{SWSA}-{CWSA}-{LWSA}-{GLSA}$$
(1)

where GWSA indicates the GWS anomaly, TWSA denotes the TWS anomaly; PWSA denotes the water storage anomaly in permafrost. SWA, SWSA, and CWSA represent the snow, soil, and canopy water storage anomalies, respectively, and LWSA and GLSA indicate the lake water storage and glacier specific mass change anomalies.

Previous studies have shown that there were differences in the magnitude and trends of TWS in the TP between different GRACE products12. We collected the Jet Propulsion Laboratory mascon (JPL-M) product41, the CSR-M product RL0642, and the German Research Centre for Geosciences (GFZ)43, for GWS estimation. The average of the three products was used to estimate the temporal evolution in GWS. Moreover, we collected two reconstructed TWSA data sets44 based on precipitation and combined the CSR RL06 and JPL RL06 mascon products for comparison in China. Lake water level observations are hardly characterized by both large numbers and high temporal resolution31,45, so we compared the differences between GWSA obtained from two sets of lake water level observations. One was the monthly lake water level data for 52 large lakes (50 lakes larger than 150 km2 and two lakes larger than 100 km2) derived from Landsat archives and altimetry data45. Another was the 5-year interval lake water level data for 1132 lakes larger than 1 km2 derived from digital elevation, optical satellite imagery, and radar altimetry data13. A comparison of GWS obtained based on different GRACE solutions and lake data is presented in Table S2.

The uncertainties in GWS were calculated following Xiang et al.4, combined with uncertainties in TWS and individual hydrologic components. The uncertainty in TWSA was represented by the standard deviation of the three products (JPL-M, CSR-M, and GFZ). The uncertainty in LWSA was assumed to be 15%46. For GLSA, the uncertainty was estimated following Hugonnet et al.47. The overall uncertainties of the calculated GWS trend are the summation of the uncertainties in the linear regression and uncertainties calculated by Monte Carlo for 1000 times.

GLSA was derived from satellite and airborne elevation data47. PWSA, SWA, SWSA, and CWSA were calculated from a high-resolution land surface model simulation driving by a high-accuracy precipitation dataset48. Before the GWS estimations, we collected observations to verify the accuracy of the hydrological components using the land surface model (Text S1, Table S1, and Figs. S1, and S4–S6). Because groundwater level observations on the TP were sparse, only 19 stations were available. The observations, belonging to the same 0.5° grid, have been averaged for comparison with GWS estimates. In general, the GWS estimation in this study is consistent with groundwater well observations (Fig. S7).

Lake and glacier water storages

To calculate the glacier-specific mass changes for each basin on the TP, we multiplied the glacier elevation change47 by the glacier area and then by a density conversion factor:

$$\Delta {GM}=\mathop{\sum }\limits_{i=1}^{n}\Delta {h}_{i}\times {A}_{i}\times \rho$$
(2)

where \(\triangle {GM}\) is the glacier mass changes for a specific basin, \(\triangle {h}_{i}\) is the glacier elevation change, \({A}_{i}\) is the glacier area, the subscript i represents the ith glacier located in this basin, and \(\rho\) is the density conversion factor for converting volume change to mass change (850 kg m-3 and an uncertainty of 60 kg m-3 in this study).

The per-glacier elevation change data were derived from a global estimate of glacier mass change during 2000–2019 at 100 m resolution, which was charted based on satellite and airborne elevation datasets. The glacier polygons were derived from the Randolph Glacier Inventory (RGI 6.0)49.

Two satellite-derived lake water storage datasets were used in this study. One was the water volume, area, and water level changes of 1132 lakes larger than 1 km2 between 1976 and 2019, with a 5-year sampling interval obtained based on digital elevation data, optical satellite imagery, and radar altimetry data13. Based on the lake water level and area, the change of lake volume was calculated:

$$\Delta V=\frac{1}{3}({H}_{2}-{H}_{1})({A}_{2}{+A}_{1}+\sqrt{{A}_{2}{A}_{1}})$$
(3)

where \(\Delta V\) is the change of lake volume, \({H}_{2}\) and \({H}_{1}\) are the lake water level at the beginning and end dates of a 5-year interval, and \({A}_{2}\) and \({A}_{1}\) are lake area at the beginning and end dates of a 5-year interval. Then we converted the lake volume change to mass change using a water density of 1000 kg m-3. Changes in lake water storage for 2000–2019 were assumed to be the same as those for 2002–2018.

We also collected the water level changes of 52 lakes derived from multi-source altimetry and remote sensing data for 2000–201845. All data have been re-processed as the anomalies relative to the time-mean baseline from January 2004 to December 2009 and presented as the equivalent water-thickness. The ratio between lake area and grid cell area was used to transform the variation in lake water storage data for individual lakes into gridded data representing a grid-averaged water thickness.

Land surface model simulations

We conducted a high-resolution land surface simulation to estimate changes in soil, snow, and canopy water storage. The model used in this study was the land surface model of the Chinese Academy of Sciences (CAS-LSM)50,51, which is the land component of the Flexible Global Ocean-Atmosphere-Land System model grid-point version 3 (CAS-FGOASL-g3)52. It considers groundwater lateral flow (GLF) and changes in the depth of frost and thaw fronts (FTF) as optional configurations based on CLM4.553. More details concerning these new parameterizations can be found in Xie et al.50,51. This model considers the FTF module, which improves the simulations of soil temperature and water phase change processes51,54. We used the change of soil ice to represent the effect of frozen soil instead of the empirical formulae4,10 (Text S1.3 and S2.1). The GWS modeled by CAS-LSM was not directly used in this study due to the lack of consideration of complex surface-groundwater interactions and active faults influences.

The simulation was established from 2000 to 2018 and run at a resolution of 0.025° (approximately 3 km at the equator) for both latitude and longitude over the TP region, and the time steps were set to 1800 s. The meteorological forcing dataset was obtained from the third pole meteorological forcing dataset (denoted as TPMFD)48,55. It merges the atmospheric simulation-based ERA5_CNN with in situ precipitation observations from more than 9000 stations. The surface air pressure, 10-m wind speed, 2-m air temperature, and 2-m specific humidity are generated similarly to the precipitation data, with the downward shortwave generated by merging ERA5 with a satellite-based solar radiation product. The initial value of the groundwater table depth was obtained from Fan et al.8. Land cover data were extracted from the Multi-source Integrated Chinese Land Cover map available at a 1 km resolution56. The soil texture, including the percentage of soil sand and clay in the model was replaced with the soil characteristics database, which was developed for use in the land surface modeling with a 1-km resolution57. The leaf area index with a 1-km resolution was obtained from Yuan et al.58. To reach a model equilibrium state (annual changes in model outputs for terrestrial hydrologic variables were less than 0.1%), we conducted a spin-up run before the formal simulation for 100 years by recycling the meteorological forcing.

In addition, we collected in situ observations of soil moisture, snow depth, and active layer thickness to validate the model performance. The results showed that the hydrological components simulated in this study agreed with observations well (see Text S1.3 for detailed discussion).

Regression method for estimating contributions

The regression method, used to calculate the contributions of precipitation, temperature, and shortwave radiation to GWS, can be written as follows:

$$Y={a}_{1}{X}_{1}{+a}_{2}{X}_{2}+{a}_{3}{X}_{3}+b$$
(4)

where \(Y\) is the standardized annual GWS, and \({X}_{{\rm{i}}}\) is the standardized annual precipitation, temperature, and shortwave radiation at each grid cell. \({a}_{{\rm{i}}}\) and \({b}\) denote the regression coefficients and intercept. The contributions were calculated with regression coefficients ranging between 0 and 1 as follows:

$${{Con}}_{i}=\frac{\left|{a}_{i}\right|}{\mathop{\sum }\nolimits_{i=1}^{3}\left|{a}_{i}\right|}$$
(5)

The variance inflation factor (VIF) for each explanatory variable (Table S3) was calculated to quantify whether multicollinearity exists:

$${V{IF}}_{i}=\frac{1}{1-{R}_{i}^{2}}$$
(6)

where \({{{R}}}_{{{i}}}^{2}\) is the coefficient of determination obtained by regressing the \({X}_{{\rm{i}}}\) on the remaining \(X\) variables. 0 < VIF < 10 indicate that there is no multicollinearity; 10 ≤ VIF < 100 indicate that there is strong multicollinearity; VIF ≥ 100 indicate that there is severe multicollinearity. Results show no multicollinearity exists between explanatory variables in Eq. (4).

GWS projection using deep learning and climate model

In this study, a long short-term memory neural network (LSTM) method coupled with a hybrid empirical model decomposition (EMD) was used to predict the GWS changes on the TP. LSTM is a type of recurrent neural network used to predict the time-dependent variables59. The EMD was used to separate the GWSA time series into different components, and the LSTM was used to predict each component extracted by the EMD. Then, each predicting component was combined as the predicted GWS on the TP.

In this study, three climate factors (precipitation, temperature, and shortwave radiation) derived from 11 CMIP6 climate models were chosen to predict GWS based on existing studies12. Owing to the challenge of reliable evapotranspiration estimation on the TP12 and the scarcity of evapotranspiration flux towers, as well as the dependence of evapotranspiration on precipitation and radiation, solar radiation was chosen.

Monthly GWSA under different climate scenarios were evaluated with two or three steps (months) ahead of predictions23. The EMD-LSTM used in this study comprised two LSTM layers of 64 neurons, followed by a dropout layer to prevent overfitting. A fully connected output layer was set atop the LSTM layer to improve the learning ability of the model. The loss function was the mean squared error. We leveraged the Adam optimizer to find optimal values of other hyperparameters (epochs and batch size) by running a 5-fold cross-validation grid search. For each basin, we randomly selected 60% of the target data for training, 20% for validating, and the remaining 20% was used for testing. The performance of EMD-LSTM is shown in Fig. S22. The EMD-LSTM model can accurately reproduce the GWSA on the TP during the historical period (2002-2018), and the Nash efficiency coefficient (NSE) value between the EMD-LSTM results and estimated GWSA is greater than 0.7 for 80% basins (Fig. S22). The EMD-LSTM model was implemented using Python 3.10, the deep-learning framework TensorFlow, and its Keras API.

Bias correction of climate model simulations

Precipitation, temperature, and surface shortwave radiation derived from 11 climate models (Table S4) were used as the input to the EMD-LSTM model. Because of the insufficient knowledge of the atmospheric physical processes and the simplification of the parameterization scheme, the outputs of climate models included unavoidable biases, requiring bias correction. Daily precipitation, temperature, and surface shortwave radiation data were obtained from CMIP6 in both historical (2002–2018) and prediction periods (2019–2100) under future scenarios (combination of the Shared Socio-economic Pathways (SSPs) and the Representative Concentration Pathway (RCP)). Three scenarios (SSP1-RCP2.6, SSP2-RCP4.5, and SSP5-RCP8.5, abbreviated as SSP126, SSP245, and SSP585, respectively) were selected to represent the low-medium-high level development with low, intermediate, and high vulnerability, respectively.

The delta method assumed the local (for example, a grid cell) bias was constant over time and was used for temperature and shortwave radiation correction:

$${T}_{{bc}}={T}_{m}-\left(\overline{{T}_{m,h}}-\overline{{T}_{{ref},h}}\right)$$
(7)

where bc is the bias-corrected daily data (temperature or surface shortwave radiation); m and ref represent daily data from the climate model simulation and reference data from TPMFD, respectively; h represents the historical period; and an overbar represents the mean value of the historical period. The quantile mapping method was used to correct the distribution bias of daily precipitation data:

$${P}_{{bc}}={{F}_{{ref}}^{-1}({F}_{m}(P}_{m}))$$
(8)

where P is the daily precipitation; F is the cumulative distribution function, and F−1 is the inverse cumulative distribution function; bc and m represent bias-corrected precipitation data and climate model simulations, respectively; and ref is reference data from TPMFD. For each month, we used all daily data from the month to generate the cumulative distribution function. The correction was performed separately for the historical and predicted periods.

In comparison with the reference data TPMFD, the delta method reduced the root mean square error (RMSE) of simulated temperature and radiation (Fig. S20). A relative decrease in RMSE reached 24% and 58% for temperature and radiation over the whole TP. Because the quantile mapping method corrected the distribution of precipitation data, we also calculated the correlation coefficient (R) between the bias- and unbiased-corrected precipitation data and the reference precipitation data TPMFD (Fig. S21). The quantile mapping method significantly improved the R value between simulated and reference precipitation data for the Indus and Tarim river basins (increasing approximately 68 and 126%) and greatly reduces the RMSE of precipitation data on the southern TP (approximately 41–65% for Mekong, Salween, Brahmaputra, and Ganges river basins).