Past and future adverse response of terrestrial water storages to increased vegetation growth in drylands

Abstract

The response of terrestrial water storages to dryland vegetation growth remains poorly understood. Using multiple proxies from satellite observations and model outputs, we show an overall increase (decrease) in vegetation growth (terrestrial water storages) across drylands globally during 1982–2016. Terrestrial water storages in greening drylands correlate negatively with vegetation growth, particularly for cropland-dominated regions, and such response is pronounced when the growth rate of vegetation productivity is high. Reduction in terrestrial water storage is dominated by precipitation and evapotranspiration variability rather by than runoff. We predict reduction in terrestrial water storage of 41–84% by 2100, accompanying expansion of drylands by 4.1–10.6%. Our findings, which indicate sustained adverse response of terrestrial water storage to vegetation growth in drylands, highlight the need for concerted planning for balanced ecological restoration, agricultural management, and water resource utilization that will affect 5.17 billion people, 64% of whom live in developing countries.

Introduction

Drylands are generally characterized by water shortage, infertile soil, and scant coverage of vegetation. Such regions include water-scarce environments such as deserts, croplands, grasslands, and shrublands that account for approximately 40% of the global land area¹. Drylands provide food, energy, water, and ecosystem services; consequently, their long-term sustainable development is a critical issue. Furthermore, intensive human interventions and extreme climatic regimes accompanying global climate change are expected to render dryland ecosystems even more sensitive to the growing dangers of extreme events such as drought and heatwave^2,3.

Crop and livestock productions are the principal economic activities in dryland regions, sustaining over 38% of the global human population. Therefore, regional vegetation growth is critically important⁴, and high-quality long-term spatiotemporally consistent records of vegetation growth are of considerable relevance for assessing such environmental ecosystems. While ground measurements and modeling evidence indicate reduced vegetation productivity in arid and semi-arid regions over recent decades^5,6, earth observation studies have revealed a more nuanced situation with signs of both decrease and increase^7,8. In particular, despite massive swaths of deforestation globally, such as in the dry woodlands of Africa and Asia⁹, satellite observations have revealed earth’s greening in recent decades as a consequence of substantial human interventions¹⁰. For example, China, India, Australia, and regions of the Sahel in North Africa have exhibited a broad trend of greening^11,12. Collectively, the vegetation dynamics of drylands reported in the literature appear inconsistent, reflecting the facts that areas of greening and browning are highly geographically and temporally heterogeneous, and that worldwide vegetation dynamics are complicated by the nexus among environmental conditions, climatic elements, and human activities. Consequently, little consensus exists currently regarding the complete picture of how and to what extent dryland vegetation trends are altered globally in the context of continuous human interventions and climate change.

Sustainable growth of vegetation generally requires a continuous supply of water; therefore, the pressure on water resources demanded by anthropogenic and ecological activities continues to grow. Vegetation restoration projects on the Loess Plateau, in China, and in northern India have nearly reached the limit of the soil water carrying capacity for vegetation^13,14. The underlying mechanisms of the hydrological cycle in drylands are complex, involving multipath interactions, feedbacks, and sensitivities between vegetation, atmosphere, and hydrology^15,16,17. The question of whether vegetation variation affects soil drying or wetness is a subject of active debate in relation to global drylands. Several studies found that vegetation greening increases evapotranspiration (ET), resulting in excessive soil water demand and more ecological drought^18,19, while others reported that vegetation has ecological functions of water conservation and water use efficiency regulation, and that precipitation can compensate for soil moisture (SM) loss due to increased ET^20,21,22. This controversy relates to the variety of assessments based on different available datasets from different regions and on different spatio-temporal scales. Furthermore, a knowledge gap remains regarding the response and shift of terrestrial water storage (TWS) in relation to future vegetation dynamics, particularly under the background of global climate change^23,24,25, which restricts our ability to manage agricultural forestry and animal husbandry.

Herein, we examine the trends of vegetation growth and terrestrial water storages of earth’s drylands over the past four decades, and assess the response of terrestrial water storages to vegetation growth. Drylands were confined to the regions with the values of aridity index (AI) (ratio of annual precipitation to annual potential evapotranspiration) smaller than 0.65²⁶. Hyper-arid (AI < 0.03), arid (0.03 ≤ AI < 0.2), semi-arid (0.2 ≤ AI < 0.5), and sub-humid (0.5 ≤ AI ≤ 0.65) regions were further identified (Supplementary Fig. 1). Five proxies derived from satellite-retrieved vegetation products (i.e., the normalized difference vegetation index (NDVI), leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FPAR), gross primary productivity (GPP), and vegetation optical depth (VOD)), and an ensemble of satellite- and model-derived terrestrial water storages (i.e., terrestrial water storage anomaly and SM), are used to delineate ecohydrological regimes (“Methods”). We also quantify the contributions of precipitation, ET, and runoff variability to the terrestrial water storages that are manifested as changes in vegetation growth. Finally, we project the response of terrestrial water storages to vegetation dynamics by combining observational datasets and earth system model simulations based on an ensemble machine learning approach to identify how such responses are likely to evolve toward the end of the twenty-first century.

Results

Trends of vegetation growth and terrestrial water storages

We quantified the trends of vegetation growth proxies and terrestrial water storages during 1982–2016. From the NDVI, we can infer that vegetation growth increased during this period over a considerable fraction of drylands (Fig. 1a), i.e., 74.35% (39.78% at the 95% confidence intervals) of drylands experienced a trend of increase, while 5.46% experienced a significant decrease. Such a trend of increase is supported by a series of other satellite-derived variables, i.e., LAI, FPAR, GPP, and VOD (Supplementary Fig. 2). These proxies indicate that vegetation growth increased over 41.57% (16.54–59.75% range) of dryland areas and decreased over 6.76% (2.15–9.34% range) at the 95% confidence intervals. Upward trends with large magnitude were detected by VOD and FPAR, whereas the magnitude of the trend captured by GPP was low because vegetation productivity did not increase in proportion to greening, probably owing to climatic stress²⁷. Hotspots of positive trends in vegetation growth occurred in regions of South and East Asia (e.g., northern India, and northern China), Americas, Sahel, and high latitude areas bordering the subarctic regions where rising CO₂ fertilization and large-scale human interventions are evident^28,29.

Fig. 1: Variations and trends of vegetation growth proxies and terrestrial water storages.

a–d Spatial distribution of trends regarding the monthly normalized difference vegetation index (NDVI), terrestrial water storage anomaly (TWSA), ERA root zone soil moisture (RZSM), and ESA surface soil moisture (SSM), respectively, based on monthly Seasonal and Trend decomposition using Loess (STL)-decomposed components (see Supplementary Figs. 2 and 3 for complementary vegetation growth inferred from other indicators). The magenta contours indicate the dryland boundaries, and the cross mark denotes values significant at the 95% level. The white indicates no available data or bare lands. In each trend distribution map, the boxplots in the bottom left show the magnitude of trend over different types of drylands where HA, A, SA, and DSH are short for hyper-arid, arid, semi-arid, and dry sub-humid, respectively; the green or blue lines denote medium values, and boxes and whiskers represent 25−75% and 5−95% percentiles. The plot in the right shows the mean (solid line) and the standard deviation (shading) of zonal mean trend of dryland. e, f Decadal trends of vegetation growth and terrestrial water storages based on ensemble empirical mode decomposition (EEMD)-decomposed components of region-specific time series. The trend was measured with the Sen’s slope and Mann–Kendall test. Only values significant at the 95% level are shown.

Regarding terrestrial water storages, clear trends of decrease dominated dryland areas (Fig. 1b), and 56.45% of dryland areas witnessed a decreasing tendency of TWS anomaly (TWSA) (34.66% at the 95% confidence intervals). Negative trends were observed for SM (Fig. 1c, d), with 59.27% of dryland areas affected in terms of root zone SM (RZSM) (e.g., western and northern China, northern India, the Mediterranean, and western U.S. and southern America), despite 17.34% experiencing a significant trend of increase in monsoon regions such as South Africa and Australia, which was mostly attributed to atmospheric moisture cycling from oceans³⁰. Such a trend of decrease was manifested in the records of the Global Land Evaporation Amsterdam Model (GLEAM) RZSM and surface SM (SSM) (Supplementary Fig. 3). Meanwhile, we observed that approximately 31.91% (12.54–40.87% range) of dry regions exhibited a significant positive trend in vegetation growth and a negative trend in terrestrial water storages. This proportion increased to 33.79% (16.34–47.65% range) for semi-arid and dry humid regions, particularly in northern latitudes, consistent with the expectation that dry–humid transitioning regions could experience both water- and energy-limited conditions^18,31. It implies a potential response of terrestrial water storages to vegetation growth, which could be corroborated by the satellite-observed sensitivity of the TWSA and SSM to LAI³² (–3.5 mm and –0.01 m³ m⁻³ decade⁻¹ per m² m⁻², respectively) over global drylands (Supplementary Fig. 4).

An increasing trend in vegetation growth (i.e., NDVI, FPAR, LAI, VOD, and GPP) and a decreasing trend in terrestrial water storages (i.e., TWSA, RZSM, and SSM) were captured in region-specific interannual variability (Supplementary Fig. 5). Focused on the annual trend over individual types of drylands, the contrasting pattern in trend between vegetation growth and terrestrial water storages (Fig. 1e, f), aligns with the derived patterns of trend at the grid level. We observed that vegetation growth increased by 34.4% (12.3–43.8% range) relative to the baseline level during 1982–1986, while the TWSA and SM declined by 89.2% and 9.3–21.6%, respectively.

Response of terrestrial water storages to vegetation growth

The response of terrestrial water storages to vegetation growth over greening drylands (defined as significantly increasing trend of vegetation growth at the 95% confidence level) was quantified. We used a sliding window strategy to provide the maximum Pearson correlations accounting for time-lagged effects (“Methods”). The negative response of terrestrial water storages to increased vegetation growth can be essential over greening drylands, with correlation of –0.25, –0.22, and –0.19 for TWSA, RZSM, and SSM, respectively (Fig. 2a). A general additive model (GAM), used to verify such adverse response, showed consistency with the sliding window strategy in terms of the strength of the correlation. In addition, such correlations in hyper-arid regions were small or even insignificant, hinting at the limited role of vegetation growth in water consumption across these regions. A more complete understanding of the causation between vegetation growth and terrestrial water storages can be further provided by Granger causality analysis (Supplementary Fig. 6), which demonstrated that vegetation growth accounts for 54.71% of the overall causality of terrestrial water storage across global greening drylands.

Fig. 2: Response of terrestrial water storages to vegetation growth.

Regional-specific maximum Pearson correlations between the vegetation growth proxies (i.e., NDVI, leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FPAR), gross primary productivity (GPP), and vegetation optical depth (VOD)) and terrestrial water storages (terrestrial water storage anomaly (TWSA), ERA root zone soil moisture (RZSM), ESA surface soil moisture (SSM), GLEAM RZSM and GLEAM SSM) based on monthly Seasonal and Trend decomposition using Loess (STL)-decomposed components over the a greening drylands and b non-greening drylands. The circle indicates the adjusted determination of coefficients (R²) calculated with a generalized additive model (GAM). c Trends of water storages as a response to vegetation growth trends over greening drylands. The shading areas represent one standard deviation around the mean. The vertical dashed line identified by turning point analysis indicates that the negative response becomes stable primarily on its right. Vegetation growth trends in (a–c) are significant at the 95% level. d Spatial distribution of regions substantially impacted by water stress under vegetation greening based on the thresholds identified in (c) (vertical dashed line). The blue box indicates the focused aggregated regions.

To detect whether the response of terrestrial water storages to vegetation growth change is causal, we examined such responses over non-greening drylands and found an obvious reduction in the correlations between vegetation growth and terrestrial water storages when compared with those over greening drylands (Fig. 2b). This reduction means that the likelihood of impacts from other unrevealed confounders, such as rising CO₂ and extreme climatic events, on terrestrial water storage variations is inconspicuous over the global scale, thereby signifying the substantial effects of increased vegetation growth over greening drylands.

For elaborated understanding of how terrestrial water flux variability is subject to dryland vegetation dynamics, we inspected the interconnections of the trends. As illustrated in Fig. 2c, decreased terrestrial water storages occurred mostly at the high end of the increased trend in vegetation growth, where water availability (e.g., precipitation and irrigation) is relatively abundant for maintaining vegetation growth^33,34. These regions, characterized by noticeable adverse response of terrestrial water storages to vegetation growth (Fig. 2d), are mostly distributed in areas close to transitioning between dry and humid conditions such as northern China, northern India, north western United States, and southwestern South America. Conversely, the trends of terrestrial water storages exhibited pronounced fluctuations toward the low end of the vegetation growth trend (left of the dashed cyan line in Fig. 2c), where regions are severely constrained by water supply, hinting at the absence of any uniform effect on water supply from the observed greening. Such fluctuations could be partly related to the complex photosynthetic biochemical activity and water use efficiency regulation controlled jointly by water and energy availability alongside ecosystem feedbacks^33,35,36, as reflected in variables such as the vapor pressure deficit and solar radiation (Supplementary Fig. 7). This finding corroborates earlier studies and underlines the necessity of involving water–vegetation–climate feedback in explaining hydrological response³¹. Collectively, greening-induced decline in TWS is subject to a progressive constraint, but is substantial and distinct over different dryland regions.

Attribution of terrestrial water storages

Increased vegetation growth is expected to alter the global and regional water cycles because it generally requires water to maintain stomatal thermodynamics and photosynthesis³⁷. Therefore, how the hydrological cycle components relating to increased vegetation growth contribute to terrestrial water flux change desires a comprehensive investigation. Given the complex mechanism connecting to the regional atmospheric circulation and land–atmosphere interactions, the trends of hydrological cycle components were first investigated using multiple datasets³⁸. An overall trend of decrease in precipitation since 1982 has been observed over global drylands, particularly in central South America and South Africa (Fig. 3a and Supplementary Fig. 8). Meanwhile, the trends of increase in ET have been pronounced (Fig. 3b and Supplementary Fig. 9), accounting for 59.67% (39.54–72.16% range) of global drylands, particularly in western North America and northern Asia, which experience additional water consumption for plant growth in addition to the continued impact of climate change. Regarding runoff, an overall trend of decrease was found over drylands but its magnitude was low, partly owing to its small absolute change (Fig. 3c and Supplementary Fig. 10). Nuanced information derived from annual time series (Fig. 3d and Supplementary Fig. 11) revealed the trends of region-specific precipitation, ET, and runoff to be –1.095, 0.674, and –0.394 mm per decade, respectively, over global drylands.

Fig. 3: Trends and variations of hydrological cycle components and their contributions to terrestrial water storages.

a–c Spatial distribution of trends of CRU precipitation, GLEAM evapotranspiration (ET), and ERA5 runoff based on monthly Seasonal and Trend decomposition using Loess (STL)-decomposed components (see Supplementary Figs. 8–10 for other complementary datasets). The magenta contours indicate the dryland boundaries, and the cross mark denotes the significance at the 95% level. In each trend distribution map, the boxplots in the bottom left show the magnitude of trend over different types of drylands where HA, A, SA, and DSH are short for hyper-arid, arid, semi-arid, and dry sub-humid, respectively; the blue or green lines denote medium values, and boxes and whiskers represent 25−75% and 5−95% percentiles. The plot in the right shows the mean (solid line) and the standard deviation (shading) of zonal mean trend of dryland. d Region-specific time series of annual precipitation, ET, and runoff for dryland, with the quantitative trends provided in Supplementary Fig. 11. e Partial correlation between region-specific precipitation, ET, runoff, and the terrestrial water storages (i.e., terrestrial water storage anomaly (TWSA) ERA root zone soil moisture (RZSM), ESA surface soil moisture (SSM), GLEAM RZSM, and GLEAM SSM), respectively, for the greening drylands represented by different vegetation growth proxies (i.e., normalized difference vegetation index (NDVI), leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FPAR), gross primary productivity (GPP) and vegetation optical depth (VOD)). All the correlations in (e) are significant at the 95% level.

Spatially, we observed that 42.56–69.88% (13.23–22.86% at the 95% confidence intervals) of greening dryland areas exhibited both reduced precipitation and increased ET. It implies a link between greening and the hydrological cycle. Such potential connection over greening dryland regions could be preliminarily characterized by the Pearson correlation and GAM analysis (Supplementary Figs. 12 and 13), giving confidence that the attribution of terrestrial water storage changes could be demonstrated with hydrological cycle components. Nevertheless, hydrological cycle variability may be modified by vegetation growth through water vapor transpiration and atmospheric moisture circulation^39,40 (Supplementary Fig. 14). We therefore used a partial correlation approach to quantify the contribution of TWSA, RZSM, and SSM to precipitation and ET (Fig. 3e), which are apparent over greening drylands, with average standard coefficients of 0.49, 0.27, and 0.23 for precipitation, respectively, and of –0.37, –0.29, and –0.22 for ET, respectively. Lower standard coefficients are found in hyper-arid lands, consistent with the constrained response of terrestrial water dynamics to vegetation growth.

Notably, other possible mechanisms such as the underlying catchment type, hydraulic diversity, legacy effects, and soil properties could also shape terrestrial water dynamics over greening drylands. For example, afforestation in some dry sub-humid regions such as China and Europe could induce increased precipitation, thereby intensifying TWS and the water-retaining capacity of soils^20,21,41. Therefore, the derived standard coefficients between terrestrial water storages and runoff that are relatively small may plausibly hold for all runoff datasets. Increase in ET is likely to result in a reduction of local runoff, and some studies have reported substantial correlation between TWS and runoff in monsoon regions^42,43; however, our work suggests that the role of runoff in modulating water depletion trends might be overstated or at least it might be region-dependent. This could be connected to the fact that changes in local runoff in dryland areas are minor and less responsible for climate variations⁴⁴, which also could be offset by regional atmospheric circulation effects^37,45.

Responses regarding to specific vegetation types

Given that vegetation type can shape ecohydrological patterns through physiological functions, plant community composition, and human management practices^46,47, we conducted formal analysis stratified by vegetation types. Increasing trends in vegetation growth proxies were found for most vegetation types, especially in semi-arid and dry sub-humid lands (Fig. 4a), with the most pronounced trends in cropland-dominated regions. Such positive trends align with the fact that approximately 33.44% (28.75–37.63% range) of greening drylands were from cropland at the 95% confidence level, mostly in the Northern Hemisphere. This finding coincides with earlier studies, such as Gao, Liang⁴⁸, which highlighted that cropland greening contributes nearly one-third of Northern Hemisphere greening, and that such a pattern is more pronounced in less developed regions than in developed regions. Conversely, trends of decrease in TWSA were found for most vegetation types over dryland regions (Fig. 4b), except forest that has hydraulic memory of groundwater systems, buffering the climatic and anthropogenic influences on terrestrial water dynamics⁴⁶. As for SM, a trend of decrease was observed for nearly all vegetation types, with larger declines for cropland. These trends mean that cropland and herbaceous vegetation tend to be more related to hydrological regimes, whereas wood vegetation might be more energy-limited, specifically evidenced by the co-occurring of terrestrial water deficit and vegetation greening in croplands, occupying 24.1–37.6% of the greening areas and spatially distributed mostly in the western United State, northern China, Pakistan, and northern India.

Fig. 4: Trends and correlations of vegetation growth and terrestrial water storages under varying vegetation types.

a Region-specific trends of vegetation growth proxies (i.e., normalized difference vegetation index (NDVI), leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FPAR), gross primary productivity (GPP) and vegetation optical depth (VOD)). b Region-specific trends of terrestrial water storages (i.e., terrestrial water storage anomaly (TWSA), ERA root zone soil moisture (RZSM), ESA surface soil moisture (SSM), GLEAM RZSM and GLEAM SSM) under five primary vegetation types: cropland (CRO), forest (FOR), shrubland (SHR), grassland (GRA) and sparse vegetation (SPV). c Region-specific maximum correlation between vegetation growth and terrestrial water storages over greening dryland based on monthly Seasonal and Trend decomposition using Loess (STL)-decomposed components. All the trends and correlations are significant at the 95% level.

The opposing trends between vegetation growth and terrestrial water storages under specific vegetation types suggest an existing intrinsic connection. Therefore, further analysis was conducted to quantify the strength of the correlation between them across greening regions (Fig. 4c). We found negative correlations in drylands regarding most vegetation types, with correlation of –0.44, –0.32, and –0.29 for the TWSA, RZSM, and SSM, respectively; this correlation was more pronounced for semi-arid and dry sub-humid regions. A consistent pattern in magnitude was found for the adjusted determination of coefficients (R²) using the GAM (Supplementary Fig. 15). Lower correlations in forest areas could be partly explained by the regional water availability potential as a response to atmospheric cycling or lateral transport of evaporated water supplying additional water³⁸. By comparison, higher correlations were identified in cropland areas, implying that such regions are subject more to dry conditions. This discrepancy between forests and crops is intrinsically related to the lower water coping capability of cropland in comparison with that of woody vegetation, attributable to shallower roots and hence less access to water in deeper soil⁴⁹ (Supplementary Fig. 16). Such an explanation could be corroborated by the maximum correlation coefficient from the sliding window strategy, that is, the terrestrial water storage decrease was highly synchronous with vegetation growth increase in cropland regions with time lags ranging from 1 to 5 months, whereas a varying time lag of up to 13 months was asynchronous in forest areas.

In-depth analysis with regard to the simulation dataset from the Community Earth System Model Version 2 (CESM2) might provide a direct means for verifying the identified negative impact of increased crop productivity on terrestrial water storages. We observed a visible positive trend of crop productivity for the two scenarios in CESM2: one with human land management (i.e., irrigation) switched on and the other with it switched. Comparison of these two experiments revealed that irrigation could lead to greater decline in terrestrial water storages (by 3.4–78.1%) over croplands (Supplementary Fig. 17), confirming the adverse response of decreased terrestrial water storages to increased productivity.

Future response of terrestrial water storage

Global change might result in ecosystem water limitation and shifts over dry regions³¹. By leveraging four future emissions pathways from the Coupled Model Intercomparison Project phase 6 (CMIP6) (“Methods”), we observed that the total area percentage of global drylands is projected to increase by 7.2% (4.1–10.6% from low- to very high-emission scenario) during 2015–2100 relative to those during 1982–2014 (Fig. 5a). Spatially, approximately 1.83 million km² of land will be converted to drylands by the end of twenty-first century, affecting more than 5.17 billion people. In this sense, water constraint and depletion in drylands are expected to be major concerns in the future, having implications for issues such as land degradation, food scarcity, and biodiversity reduction. Furthermore, we found that the vegetation growth proxies—LAI and GPP—are projected to increase during 2015–2100 (Fig. 5b, c). Concurrently, decreased water yield (precipitation minus ET) will persist and dominate despite increasing supply from precipitation (Supplementary Fig. 18). It hints that the adverse response of terrestrial water storages to greening could be maintained in the future, especially given that human activities (e.g., irrigation water withdrawal) will become increasingly intensive in many areas such as North China, South Asia, the Middle East, and North America^50,51.

Fig. 5: Projections of vegetation growth and terrestrial water storages in dry regions.

a Area percentage of various climate zones from hyper-arid to humid during the historical period (1982–2014), and future scenarios (2015–2100) based on the multimodel ensemble mean of CMIP6 simulations. b, c Observations and multimodel ensemble mean of the leaf area index (LAI), gross primary productivity (GPP) anomaly. d–f Observations and machine learning-predicted root zone soil moisture (RZSM) anomaly, surface soil moisture (SSM) anomaly, and terrestrial water storage (TWS) anomaly based on multimodel ensemble mean. The boxplots show the decadal trends based on the ensemble empirical mode decomposition (EEMD) applied to the four emissions scenarios. The circle above the boxplot delineates the trends based on the least absolute shrinkage and selection operator (LASSO)-predicted outcomes. The shading represents standard deviation. All the trends are significant at the 95% level. g, h Schematic showing interactions between vegetation growth and terrestrial water storages. Numbers in parentheses in (g) correspond to the percentage change in the period 1987–2014 relative to the level of the baseline period (1982–1986), based on observations; numbers in parentheses in (h) correspond to the percentage change relative to the historical level in 1982–2014 based on CMIP6 simulations.

To quantify the future response of terrestrial water storages to vegetation greening, we predicted the variations and trends of RZSM, SSM, and TWS during 2015–2100 using a machine learning approach and multimodel ensemble mean (“Methods”). Despite the inherent uncertainty in earth system models, basic agreement in the hydrological response among the four emissions scenarios was found, that is, noticeable decline in terrestrial water storages, resembling that determined from observations (Fig. 5d–f). Specifically, the anomalies of RZSM, SSM, and TWS during 2015–2100 showed sustained decrease of –84%, –41%, and –67%, respectively, in comparison with those of the baseline historical period (1982–2014). This decrease is plausible because increased vegetation resulting from elevated CO₂ fertilization and favorable energy availability could contribute to enhanced canopy transpiration, soil evaporation, and surface runoff. However, such trends of decrease trend might slow owing to certain factors such as increasing precipitation, as evidenced by the negative trends of anomalies of RZSM, SSM, and TWS during 2015–2100 accounting for 24.8% (11.3–28.6% range), 79.3% (49.2–121.9% range), and 37.2% (22.3–48.4% range) of the historical trends (1982–2014), respectively.

The negative trends in terrestrial water storages derived from the ensemble machine learning approach were further corroborated by results from least absolute shrinkage and selection operator regression, which also showed an obvious decrease in the projected terrestrial water storages (Fig. 5d–f and Supplementary Fig. 19). Thus, the uncertainties of model projections might influence the magnitude but not alter the signs of the trends in terrestrial water storages. This study established the capacity for reflecting the projected variability in terrestrial water storages, and revealed that the effects of increasing vegetation growth are projected to amplify reduction in terrestrial water storages over global dryland regions (Fig. 5g, h).

Discussion

The spatial and temporal heterogeneities and inconsistencies of observational and modeling records pose major barriers to their use in assessing ecohydrological response; therefore, unlike previous related studies, we levered multiple satellite-derived proxies of vegetation growth and a robust ensemble of terrestrial water storages to derive credible conclusions. Through use of a comprehensive ecological and hydrological gridded dataset, we found an overall trend of increase in vegetation growth over the previous four decades in drylands globally, with a concurrent trend of decrease in terrestrial water storages. On the basis of this correspondence, we detected that a considerable proportion of drylands experiencing terrestrial water decline could be attributed to vegetation growth increase, such as in South and East Asia, the Mediterranean, and in the Americas. Although such a pattern has previously been reported on regional and global scales using a variety of data sources, the extent to which the dryland ecosystem might shift remains obscure. Here, we provided evidence suggesting that increased vegetation growth could explain 19–25% of the decline in terrestrial water storages.

Fluctuation in the response of terrestrial water storages to vegetation growth increase is observed in drylands when vegetation growth change is low; that is, less water availability might be supplied. Although we cannot confirm such an underlying mechanism, greening-induced intensification of decline in terrestrial water flux is likely to be subject to a progressive constraint. It aligns with the results of earlier studies associating the sensitivity in ecohydrological regimes with the dominant role of CO₂ fertilization and nitrogen uptake^52,53. Furthermore, a threshold might exist regarding control of the interplay between vegetation growth and terrestrial water storages, corresponding to the identified turning point in structural attributes along climatic and aridity gradients^54,55. It is thus possible that once an aridity level is reached, small increases in aridity may lead to pronounced interaction among vegetation and soil properties (e.g., nutrient cycling and microbial communities), potentially modifying hydrological patterns such as water use efficiency. Although there is lack of knowledge regarding the timeframe over which an established threshold might be crossed, these thresholds might be region-specific and their values might become smaller as the water content increases. Given the complex mechanisms dominating vegetation hydrology response, research efforts that focus on physical and process models with multiple datasets are encouraged to provide more multivariate perspectives.

On the basis of confirmed trends and relationships, we infer that the substantial decrease in terrestrial water storages of distinct vegetation types could be accompanied by increase in vegetation growth, which is more pronounced for cropland-dominated regions. Such negative effects in cropland are robustly manifested in earth system models by disentangling irrigation activities. Our results imply that conflicts regarding water availability between ecosystem water demands and human food necessity could become a matter of concern, which is particularly worrisome for identified hotspots such as northern China, northern India, and the western United States, where increased vegetation growth might override any positive effect of water availability, especially considering that areas of cropland and use of irrigation water have increased during the previous four decades (Supplementary Fig. 20 and Fig. 21). Our work goes beyond current knowledge by underscoring the intense adverse response of terrestrial water storages to increased vegetation growth over croplands in global dryland areas, and by highlighting the need for adaptation measures to minimize water depletion associated with agricultural management.

The potential connection between vegetation growth and terrestrial water storages not only enhances our current understanding, but also contributes to the developments yet to come. Our findings suggest that earth system model simulations integrated with machine learning approaches could be used to establish an empirical manifestation of the future response of terrestrial water storages to vegetation growth over global drylands. From a regional aggregated perspective, we projected that TWS and SM will decrease persistently to the end of this century. Therefore, the predicted effects of increased vegetation growth over drylands are highly likely to make terrestrial water storages vulnerable to climate change and human interventions, which are less considered in current literature. Notably, the assumed input–output correlation in the machine learning-based prediction was constant; nevertheless, it could be influenced by uncertainties associated with cropland water extraction, extreme drought, and land use management, which occur extensively. Consequently, terrestrial water storages are expected to decrease, especially considering that global drylands will expand by approximately 1.83 million km2 and that the region of croplands is highly likely to increase^56,57. Thus, the importance of the regional resource carrying capacity for vegetation in agricultural management and ecological restoration cannot be overstated, and specific measures such as water withdrawal constraints and water-efficient irrigation should be encouraged, particularly in relation to developing countries that occupy 64% of the population impacted by expansive drylands.

Notwithstanding the critical implications derived, several potential caveats should be considered when interpreting our work. First, the spatial scales between the vegetation growth proxies and the terrestrial water storages are biased owing to the systematic errors and sensor configurations of the collected products. Specifically, the biases in vegetation growth can be 10–21% between the Moderate Resolution Imaging Spectroradiometer (MODIS) observations (NDVI, FPAR, LAI, and GPP) and the used dataset when focusing on the period of 2002–2016 (Supplementary Fig. 5). Nevertheless, consistent trend patterns are identified, manifesting the rationality of our results for the period 1982–2016. However, the bias in the magnitude of the model simulations could be the primary source of dataset sensitivity^58,59, and uncertainty is inevitable in the modeling of water storages and fluxes, such as RZSM, ET, and runoff. Given the above, the greening-induced hydrological responses might not be fully delineated in our study and improved representation is needed. Second, despite being known critically, the implications of other climatic and anthropogenic factors such as increases in atmospheric CO₂, aerosols, and greenhouse gas emissions on terrestrial water are not disentangled in our study. Elevated CO₂ levels could potentially ameliorate the reduction in water availability by enhancing leaf water use efficiency^60,61. This is also indirectly confirmed by the decrease in terrestrial water storages subjected to a constraint by limited water supply in arid regions. It also means that the attribution of the decline in terrestrial water fluxes to vegetation growth might be overestimated in our study. Thus, critical analysis of terrestrial water dynamics is demanded using coupled hydrological–vegetation models to disentangle the implications from elevated climatic variations and anthropogenic interventions. In addition, vegetation can impact TWS by modifying the water cycle, while terrestrial water dynamics also affect vegetation growth. Here, we only provide a preliminary overview of the response of terrestrial water storages to vegetation productivity change; further in-depth studies of their interactions are demanded. Despite these appropriate caveats and limitations, our work justifies the use of multiple global ecological and hydrological proxies in providing a scientific basis for further investigation.

Methods

Vegetation growth proxies

The period of 1982–2016 was selected in our study to record long-term vegetation growth and terrestrial water storages. We collected multiple proxies of vegetation growth from satellite-derived observations. Daily Advanced Very High Resolution Radiometer (AVHRR) NDVI, LAI, and FPAR data with the resolution of ~0.05° were used to depict vegetation growth. The Ku-band VOD, collected from multiple satellite microwave observations during 1988–2016 and owned a spatial resolution of 0.25°, was also used to depict vegetation productivity⁶². Daily NDVI, LAI, and FPAR values, averagely composited into monthly datasets, and daily VOD values, maximally composited into monthly datasets, were calculated for a given grid cell only when data availability was ≥80%. Global monthly GPP modeled at 8 km spatial resolution is provided by an improved Monteith light use efficiency model⁶³. To confirm the robustness of the long-term satellite-observed trends, Moderate Resolution Imaging Spectroradiometer (MODIS) based monthly NDVI (average of MOD13A3 and MYD13A3), LAI and FPAR (MCD15A3H), and GPP (average of MOD17A2 and MYD17A2) dataset with a spatial resolution of ~1 km were selected from the period of 2002–2016 for depicting vegetation growth based on the Google Earth engine platform.

Terrestrial water storages

Monthly terrestrial water storage anomaly (TWSA) with a spatial resolution of 0.25° was obtained from the average values of two Gravity Recovery and Climate Experiment (GRACE) mascon-based products: the Center for Space Research (CSR) and the Jet Propulsion Laboratory (JPL). The missing pixels in the available GRACE were interpolated with one linear approach. GRACE dataset are unavailable before 2003, and therefore the TWSA dataset for 1982–2002 was reconstructed using one ensemble machine learning approach that involved a neural network, multiple linear regression, and the random forest method, following the strategies employed in Liu et al.⁶⁴. We also used the GRACE‐FO TWSA during 2017–2019 as reference for evaluating the accuracy of the reconstructed TWSA and the results demonstrated that the performance of the reconstructed model was acceptable and that the model could be applied for further analysis (Supplementary Fig. 22).

Daily root zone SM (RZSM) and surface SM (SSM) were used for representing terrestrial water storages, provided by the reanalysis products of European Centre for Medium-Range Weather Forecasts (ECMWF) ReAnalysis 5 version (ERA5) and the satellite product of the European Space Agency’s Climate Change Initiative (ESA CCI), respectively. For ERA5, with a spatial resolution of ~0.1°, the RZSM of the soil layers is weighted between 0 m and 1 m. Regarding satellite ESA CCI, with a spatial resolution of 0.25°, only the thin SSM is sensed. We also collected daily RZSM and SSM from the Global Land Evaporation Amsterdam Model (GLEAM)⁶⁵. The composite method based on average values was used to convert daily ERA, ESA, and GLEAM SM into monthly values.

Precipitation, evapotranspiration (ET), and runoff

A variety of reanalysis and model simulations during the study period were collected to depict the hydrological cycle components, i.e., precipitation, ET, and runoff. We collected monthly precipitation from four sets of products, i.e., Climatic Research Unit Timeseries (CRU TS), Global Precipitation Climatology Center (GPCC), ERA5, and TerraClimate⁶⁶. Three sets of monthly ET were collected from the products of GLEAM, FluxCom⁶⁷, and TerraClimate. The ET dataset from Numerical Terradynamic Simulation Group (NTSG)⁶⁸ covering 1982–2013 was also used. As for the runoff, we collected three sets of monthly products from the ERA5, Global Runoff Reconstruction (GRUN)⁶⁹, and TerraClimate dataset.

Earth system model simulations

We collected earth system model simulations from eight CMIP6 models (Table S1), which could provide monthly LAI, GPP, precipitation, ET, runoff, SM, and TWS, as well as air temperature, short solar radiation and wind speed for terrestrial water storages predicting and arid index calculating. These earth system models cover the historical (1850–2014) and future (2015–2100) periods, and here the period of 1982–2014 was selected as reference periods, during which corresponding observation datasets were available. Four future emission and land-use scenarios, i.e., the SSP1-2.6 (low forcing sustainability pathway), SSP2-4.5 (medium forcing sustainability pathway), SSP3-7.0 (medium- to high-end forcing sustainability pathway), and SSP5-8.5 (high-end forcing sustainability pathway) were used, representing radiative forcing stabilized at 2.6, 4.5, 7.0, and 8.5 Wm⁻² by 2100, respectively. All the selected CMIP6 outputs were interpolated to 0.5° × 0.5° spatial resolution to match the available observations.

We employed the multimodel mean of the CMIP6 simulations to overcome the uncertainty among the CMIP6 models. The ensemble mean values were calculated by weighting eight selected model members based on the root mean square error and independence scores of each model simulation relative to reference data, following earlier studies^3,70. The AVHRR LAI, light usage efficiency model-derived GPP, ERA RZSM, GLEAM SSM, and GRACE TWS were used as references. We used the ensemble mean of available precipitation, ET, and runoff products as references. The references of air temperature, short solar radiation, and wind speed were obtained from the ensemble mean of CRU, ERA5, and TerraCliamte datasets.

To check the potential effects of irrigation activities on terrestrial water flux changes over croplands, an additional dataset from the Community Earth System Model Version 2 (CESM2)⁷¹ simulations was used. The CESM2 model is built on the Community Land Model version 5 (CLM5) and incorporates the dynamics of crop growth and management practices to realistically simulate land surface fluxes in agricultural regions. Two CESM2 simulations driven using two different land management strategies were analyzed: one was conducted with a normal land management configuration that included irrigation activities, while the other was conducted with the same configuration but with cropland treated as unmanaged grassland. Comparing these two simulations with the same preprocessing procedure can disentangle the terrestrial water storage trends arising from irrigation.

Aridity index (AI)

The 30 arc seconds spatial resolution AI was provided by the Global Aridity Index and Potential Evapotranspiration Database (Global_AI_PET) that averages the global hydro-climatic data during 1970–2000. We also calculated the AI during 1982–2100 based on CMIP6 simulations to check future dryland expansion. The PET was calculated using FAO Penman-Monteith method:

$${\rm{PET}}=\frac{0.408\triangle \left({R}_{n}-G\right)+\gamma \frac{900}{T+273}{u}_{2}({e}_{s}-{e}_{a})}{\triangle +\gamma (1+0.34{u}_{2})}$$

(1)

where $\triangle$ represents the slope of the saturation vapor pressure curve versus temperature; $T$ is the mean air temperature; ${e}_{s}$ and ${e}_{a}$ are the saturation and actual atmospheric water vapor, respectively; ${R}_{n}$ and $G$ are the surface net radiation and the all-wave ground heat flux, respectively; $\gamma$ is the psychrometric constant; and ${u}_{2}$ is the surface wind speed.

We utilized CMIP6 dataset to estimate the population affected by global dryland expansions, through overlaying population density data onto dryland areas. To investigate the potential impacts of climate change on dryland populations, we analyzed the spatio-temporal changes in population density within these dryland areas over time. The global annual population is obtained from Inter-Sectoral Impact Model Comparison Project (http://clima-dods.ictp.it/Users/fcolon_g/ISI-MIP).

Auxiliary datasets

Global long-term annual datasets from the ESA Climate Change Initiative Land Cover (CCI-LC) are used to identify primary vegetation types. Considering this product only provides the dataset from 1992, we collect the land cover dataset from the Global LAnd Surface Satellite (GLASS)⁷² dataset during 1982–1991. Here, we focused on five primary vegetation types: cropland (CRO), forest (FOR), shrubland (SHR), grassland (GRA), and sparse vegetation (SPV) (Supplementary Fig. 23). Global monthly water withdrawals dataset⁷³ during 1982–2010 was collected. We collected the monthly vapor pressure deficit from TerraClimate dataset, as well as plant rooting depth and ground water table depth for further analysis.

Trend calculation

We used Sen’s slope estimator to calculate the trends of the vegetation growth proxies and terrestrial water storages. The Mann–Kendall test was used to determine the significance level. The trend (${\rm{Tr}}$) is the median slope that is calculated based on the differences among all possible pairs of month or year (i and j, i < j),

$${\rm{Tr}}={\rm{median}}\left(\frac{{V}_{j}-{V}_{i}}{j-i}\right)$$

(2)

We first calculated the pixel-level trend of the monthly time series, and only pixels occupying the time series with more than 80% completion were considered. To eliminate seasonal signals, we applied a Seasonal and Trend decomposition using Loess (STL) to each time series. Specifically, the original time series was separated into three components after decomposition: trend, seasonal signal, and the remainder. The pixel-level trend was calculated based on the STL-decomposed time series that excludes the seasonal signal and the remainder components.

We also calculated the region-specific trend of interannual variability of the vegetation growth proxies and terrestrial water storages, which was implemented on the time series generated by averaging the pixel data over each dryland region. An ensemble empirical mode decomposition (EEMD) method was used in conjunction with Sen’s slope, and the Mann–Kendall test was used to detect region-specific trends; that is, a region-specific trend calculation was implemented on the decomposed secular trend component. EEMD is an adaptive method for separating a time series into a finite number of intrinsic mode functions (${\rm{IMF}}$) and a secular trend, as the following:

$${Y}_{j}\left(t\right)={\sum \limits_{i=1}^{n}}{{\rm{IMF}}}_{{ij}}\left(t\right)+{R}_{j}\left(t\right)$$

(3)

where ${{\rm{IMF}}}_{{ij}}\left(t\right)$ and ${R}_{j}\left(t\right)$ denote the ith IMF component and residue from the noisy pseudo times series ${Y}_{j}\left(t\right)$ which is generated by adding realizations of white noise to original time series $x\left(t\right)$:

$${Y}_{j}\left(t\right)=x\left(t\right)+{g}_{j}\left(t\right),\,j=1,2,\ldots m$$

(4)

Response and attribution analysis

Response and attribution analyses were performed on the remainder parts of the STL-decomposed monthly time series. The response of terrestrial water storages to vegetation growth was first assessed using Pearson’s correlation analysis. For eliminating the influences relating to phenology and lag times, we used the maximum correlation based on one sliding window strategy. To achieve this, a parametric linear regression was established between the vegetation growth proxies and month-lagged terrestrial water storages.

$${{\rm{TW}}}_{t}={\eta }_{0}+{\eta }_{k}\times {{\rm{VP}}}_{t-k}$$

(5)

where ${{\rm{TW}}}_{t}$ is the terrestrial water storage at the specific time $t$, ${\eta }_{0}$ and ${\eta }_{k}$ are regression coefficients corresponding the k-month lag with respect to explanatory variables ${{\rm{VP}}}_{t}$. $k$ ranges from 0 to 24 (i.e., 0 represents no time effect, and 1–24 represents a one- to 24-month lag, respectively), and the lag month showing the highest determination coefficient is considered the best time lag.

We also used one nonparametric smoothing regression correlation, namely, the generalized additive model (GAM) that comprises a combination of generalized linear models and additive models, to depict the response of terrestrial water storages to vegetation growth. Although the GAM cannot provide the direction of the response, it has a powerful capacity for quantifying the magnitude of nonlinear and interactive correlations via a link function^74,75, as the following:

$$g({{\rm{TW}}}_{t})={\beta }_{0}+f({{\rm{VP}}}_{t})$$

(6)

where $g$ is the link function, ${\beta }_{0}$ is the model constant, and $f$ is a smoothing function for the explanatory variables ${{\rm{VP}}}_{t}$. The Akaike information criterion was used to select the appropriate model with respect to each explanatory variable. Here, we employed the adjusted coefficient of determination to evaluate the magnitude of the response.

The attributions of terrestrial water flux variations from hydrological cycle components were investigated first by checking the region-specific correlation between them using both GAM analysis and the sliding window strategy. Then, we then used partial least squares regression (PLSR) to calculate the contribution of changes in precipitation ($\triangle {\rm{P}}$), ET ($\triangle {\rm{ET}}$), and runoff ($\triangle {\rm{R}}$) to the changes in terrestrial water storages ($\triangle {\rm{TW}}$). One of the most intriguing aspects of PLSR is that the weights and regression coefficients of individual predictors in the most explanatory components can be used to infer correlations between the predictors and the response function.

$$\triangle {\rm{TW}}={\beta }_{0}+{\beta }_{1}\times \triangle {\rm{P}}+{\beta }_{2}\times \triangle {\rm{ET}}+{\beta }_{3}\times \triangle {\rm{R}}$$

(7)

where ${\beta }_{0}$, ${\beta }_{1}$, ${\beta }_{2}$, and ${\beta }_{3}$ refer to regression coefficients. Here, we used the PLSR-derived standardized coefficients to demonstrate the contributions.

Predicting terrestrial water storages

We predicted terrestrial water fluxes over drylands during 2015–2100 by applying a relationship established from the observations and the CMIP6 simulations under four greenhouse gas emission scenarios^2,76. In fact, terrestrial water storages are impacted not only by external climate drivers but also by anthropogenic forcing⁶⁴. Here, we assume that the anthropogenic effects related to human interventions such as government policy during the predicted period were identical to those during the historical period and that their changes were minimal. Machine learning approaches were used to delineate such projected relationships, and the multimodel ensemble means of GPP, LAI, precipitation (${\rm{P}}$), air temperature (${\rm{T}}$), ET, runoff (${\rm{R}}$), and terrestrial water storages were used as explanatory factors for delineating RZSM, SSM, and TWS. The plausibility of the selected variables was corroborated by variable importances (Supplementary Fig. 24).

To remove the bias between the observations and CMIP6 simulations, the monthly dataset was standardized to obtain anomaly values before applying the prediction procedure. To enhance confidence in the model simulations, this study selected four machine learning algorithms: random forest, support vector machine, neural network, and gradient boosting machine. The final model implemented for each grid is identified as the machine learning function (${\rm{ML}}$) that has the highest performance in the training period.

$${\rm{TW}} \sim {\rm{ML}}{({\rm{LAI}},{\rm{GPP}},{\rm{P}},{\rm{T}},{\rm{ET}},{\rm{R}},{\rm{TW}})}_{{\rm{CMIP}}6}$$

(8)

In practice, both the ensemble mean values of the CMIP6 simulations and the observed dataset were first decomposed into the trend, seasonal signal, and the remainder component using STL. A relationship was then established between each STL-based component of the CMIP6 simulations and the observed dataset for the period 1982–2014. This relationship can reveal the potential physical mechanisms underlying each principal component, although it might not indicate causality⁷⁷. Finally, we applied the established correlation to the CMIP6 simulations during 2015–2100 and added them to predict SM and TWS. Here, we used the ERA RZSM, GLEAM SSM, and GRACE TWSA as the observed terrestrial water storages for training.

The accuracy of the predicted SM and TWS anomaly was evaluated using reference values of ERA RZSM, GLEAM SSM, and GRACE TWSA during 2011–2014 when using the dataset during 1982–2010 for training. Pearson’s correlation coefficients and Nash–Sutcliffe efficiency indicated the applicability and reasonability of using ensemble machine learning models in predicting terrestrial water storages (Supplementary Figs. 24 and 25). In addition, to corroborate the machine learning-based predicted SM and TWS, we compared the output with the results based on a classic but extensively used linear mode, that is, the least absolute shrinkage and selection operator regression.