Evapotranspiration stress intensifies with enhanced sensitivity to soil moisture deficits in a rapidly greening China

Amidst drastic environmental changes, the intricate interplay and feedback mechanisms in the water–vegetation–atmosphere nexus experience alteration. Previous research primarily centres on the responses among variables within this system, with little known about whether and how these responses (sensitivities) change. Here, we employ the Evapotranspiration Stress Index (ESI) to represent the equilibrium of the nexus and develop a memory dynamic linear model based on Bayesian forward filtering. The model takes into account the carry-over effect in the “dry gets drier” self-amplify loop, allowing for a more effective estimation of the ESI time-varying sensitivity to associated influencing factors. Our analysis reveals that from 1950 to 2020, mainland China experienced a notable 4.74 % escalation in evapotranspiration stress. Surface soil moisture serves as a primary driver, whose sensitivity to ESI has surged by 9.49 % in the last decade compared to the early 2000s. Vapour pressure deficit (VPD) and leaf area index (LAI) also exerted a substantial role, with their sensitivities fluctuating approximately 22.91 % and −45.77 %, respectively. Moreover, the greening pace is linked to an increase in soil moisture sensitivity and a decrease in VPD sensitivity, suggesting that rapid greening may alter the ecological resilience against soil deficit and atmospheric drought. In comparison, the widely used moving window multiple linear regression (MLR) significantly overestimates sensitivity fluctuations, necessitating prudent interpretation of numerical estimates in related research findings. Our findings underscore the spatiotemporal variations in sensitivity, enriching the comprehension of ecosystem reactions to external factors, and offer essential insights for advancing greening endeavours.

1 Introduction

Over the last 2 decades, approximately one-third of land globally has experienced vegetation greening (Chen et al., 2019). Concurrently, 53 % to 64 % of these regions have undergone atmospheric desiccation (Yuan et al., 2019). However, research has also revealed a rise in available water resources (precipitation supply minus evapotranspiration loss) in ∼ 45 % proximate and downwind zones (Cui et al., 2022). Contrary to greening, satellite-based studies have indicated a 0.4-fold expansion in regions shifting from greening to browning (De Jong et al., 2013; Liu et al., 2023). Amidst the ongoing debates over greening or browning, and wetting or drying, it is unequivocal that the intricate interactions and feedback mechanisms among water, vegetation, and climate are altered.

Figure 1 Conceptual diagram of evapotranspiration stress and its self-amplify loop. (a) ESI represents the water–vegetation–atmosphere equilibrium, defined as the ratio of actual evapotranspiration (ET_a) to potential evapotranspiration (ET_p). When using the parameters of crops assumed in the PM equation to calculate ET_p (also known as ET₀, reference crop evapotranspiration in agriculture), ESI values greater than 1 may occur. (b) As regions experience more severe water deficiency or sparser vegetation cover, the balance tilts further to the right, and the ESI value decreases. A smaller ESI value indicates more severe drought conditions. (c) Schematic of the self-amplify feedback loop, where Gs represents stomatal conductance, GPP stands for gross primary productivity, and VPD denotes vapour pressure deficit. Blue lines depict the atmospheric feedback pathway, while green lines represent the terrestrial surface process pathway.

Defined as the ratio of actual to potential evapotranspiration, the Evapotranspirative Stress Index (ESI) can comprehensively reflect the equilibrium between water supply, vegetation status, and climatic conditions (Anderson et al., 2013, 2016) (Fig. 1a and b). Potential evapotranspiration (ET_p, also known as ET₀, reference crop evapotranspiration in agriculture) represents the evaporation under ideal circumstances – ample water, uniform vegetation (typically alfalfa) – and is driven solely by climate (Allen et al., 1998; Li et al., 2022b; Thornthwaite, 1948). Regions suffering from water scarcity or  sparse vegetation are subject to evapotranspirative stress (Fig. 1b), whereas ample water or vigorous  growth may result in actual evapotranspiration surpassing the hypothetical potential (Liu et al., 2022) (Fig. 1a). This occurs because ET_p calculations assume fixed canopy and surface properties for a hypothetical crop, whereas real vegetation may outperform these assumptions. Previous research has proven the high sensitivity of ESI to early-stage drought, as plants naturally curtail their water consumption before drought fully manifests (Nguyen et al., 2019; Otkin et al., 2018). Understanding ESI dynamics is thus critical for timely drought response, ecosystem stability, and sustaining crop yields.

The evolution of evapotranspiration (ET) stress includes a positive “dry gets drier” land–atmosphere feedback (Fu et al., 2022; Gentine et al., 2019; Seneviratne et al., 2010) (Fig. 1c). Regions subject to climate-induced ET stress are characterized by low soil moisture and weak vegetation vitality. Isohydric species, in response, close their stomata to conserve water, albeit at the cost of decreased photosynthesis and diminished carbon accumulation (Grossiord et al., 2020; Joshi et al., 2022). On the one side, this reduction in water vapour release into the atmosphere intensifies atmospheric dryness, further reinforcing climate constraints. On the other side, prolonged stress may lead to plant mortality due to hydraulic failure and carbon starvation, further decreasing vegetation cover (Kono et al., 2019; Mantova et al., 2022). Such a self-amplify loop highlights the necessity to consider temporal connections when exploring the ET stress dynamics. Research based on data analysis and field experiments has uncovered that such internal memory effects are particularly pronounced in water-limited areas (Liu et al., 2018; Richard et al., 2008; Xiao et al., 2024a). Statistically, the autocorrelation component can serve as a good representation of internal memory effects (Forzieri et al., 2022; Smith et al., 2022). For example, Kusch et al. (2022) employed a lagged normalized difference vegetation index (NDVI) as a representation of intrinsic-vegetation memory and compared its impact with external meteorological factors. An increasing autocorrelation coefficient signals declining ecosystem resilience or nearing critical thresholds (Liu et al., 2019; Scheffer et al., 2009). Unfortunately, much existing research has focused on external factors, neglecting the role of internal memory in ET stress dynamics (Feng et al., 2020; Fu et al., 2024b; Jung et al., 2010; Liu et al., 2020c; Peng et al., 2019). Therefore, integrating the co-regulation of external forcing factors with internal memory effects is essential for an in-depth investigation into the dynamics of ET stress and ecological responses.

Differences in regional evaporation stress patterns arise from different response to relevant factors. Soil moisture acts as the supply side of evaporation, while vapour pressure deficit (VPD), as the demand side, pulls water from soil pores and plant stomata into the atmosphere (Dong et al., 2020; Liu et al., 2020b). Increased ambient CO₂ can mitigate evaporation by enhancing water use efficiency (Li et al., 2023; Liu et al., 2017). Among external drivers, moisture and VPD are found to be particularly sensitive (Liu et al., 2020a; Sulman et al., 2016; Zhang et al., 2021a; Zhong et al., 2023), but vegetation structure, measured by the leaf area index (LAI), also influences evapotranspiration by regulating surface energy and land–atmosphere coupling (Fu et al., 2022). In addition to regional characteristic dependence, the differing sampling time spans are also important contributors to the inconsistent results. In an ever-changing world, the sensitivity of evaporation stress to environmental changes may also evolve (Hsu and Dirmeyer, 2023). However, research on the time-varying sensitivity is still in its infancy. Existing research investigates these temporal changes by segmenting time series, such as Zeng et al. (2022) using a 5-year moving window for NDVI sensitivity to precipitation or Li et al. (2022a) using machine learning to assess LAI sensitivity to soil water over 3-year blocks. However, these approaches only yield an average response for each segment, rather than the true “instantaneous or dynamic” sensitivity. Fortunately, the application of a dynamic linear model based on Bayesian forward filtering, as demonstrated by Zhang et al. (2022) in their recent research of global NDVI sensitivity to precipitation, offers a promising avenue to elucidate this problem.

Over the past 2 decades, China's mean air temperature has increased at a rate of 0.24 °C per decade, exceeding the global average of 0.20 °C per decade. Concurrent with climate warming and large-scale ecological restoration initiatives – including the Three-North Shelter Forest Program, Grain for Green Program, and Plain Greening Project (Fu et al., 2024a) – the regional leaf area index (LAI) has risen by 7.66 % (4.21 % per decade, p < 0.001), the highest rate globally (Chen et al., 2019). However, surface soil moisture has simultaneously declined at 0.08 % per decade, signalling hydrological trade-offs. These rapid environmental transformations, driven by both climate variability and anthropogenic activities, have established China as a critical natural laboratory for studying coupled hydrological–ecological–climatic interactions (Bai et al., 2020; Li et al., 2018; Zheng et al., 2022)

Here, a memory dynamic linear model (MDLM) integrates intrinsic lagged effects and concurrent extrinsic forcing to quantify both the response of the evapotranspiration stress index (ESI) to environmental changes and its temporal variability. Utilizing the MDLM, we calculate the time-varying sensitivity of ESI to key drivers at pixel scale in China over the past 2 decades, with 8 d temporal granularity. The robustness of the results is verified through a moving-window multiple linear regression model. Prior to model implementation, interpretable machine learning algorithms leveraging two classifiers are employed to screen key drivers of ESI across diverse climate–vegetation zones at 8 d, 16 d, and monthly timescales. This study aims to (a) untangle the dynamic trajectory of ET stress, (b) investigate the potential driving mechanisms of ESI, and (c) quantify the time-varying sensitivity of ESI to principal external drivers upon the incorporation of intrinsic memory effects. Our findings are expected to narrow the knowledge gap in temporal variability of evapotranspiration stress response, enhance our comprehension of the complex water–vegetation–climate interactions, and inform potential improvements for Earth system models.

2 Materials and methods

2.1 Datasets

The Evapotranspiration Stress Index (ESI) (Eq. 1) is the ratio of ET_a to ET_p. Although ET_p lacks consensus in definitions and estimation methods, the Penman–Monteith (PM) equation is adopted here due to its consistency with the ESI framework. The ERA5-derived ET_p product satisfies this criterion, delivering integrated radiation, surface temperature, and wind speed data required by the PM equation. To avoid potential overestimation of ET_p, we implement the Yang et al. (2019)-modified PM formula, which incorporates the CO₂ water-saving effect and is widely adopted in evapotranspiration studies (Lian et al., 2021; Zhang et al., 2023). CO₂ data were obtained from CarbonTracker (CT2022) for the period 2001–2020 at a 3° × 2° spatial resolution. Appendix A details the calculation procedures of ESI. Additionally, we assess the efficacy of three other mainstream evaporation products in representing ESI spatial patterns across China, benchmarking these against observational data from 26 flux stations.

$$\begin{array}{}\text{(1)}& \text{ESI}={\displaystyle \frac{{\mathrm{ET}}_{\mathrm{a}}}{{\mathrm{ET}}_{\mathrm{p}}}}\end{array}$$

In addition to the meteorological and CO₂ data utilized in ET_p calculations (Eq. A4), total precipitation (Prec), atmospheric pressure (P), and soil volumetric water content (svm) across different soil layers from ERA5 are also considered as influencing factors for the ESI. The root zone of vegetation serves as a critical interface for water–plant–atmosphere interactions (Gao et al., 2014, 2024; Wang-Erlandsson et al., 2016). Utilizing root zone depth estimations based on an 80-year drought return period from Stocker et al. (2023), we computed the spatial distribution of root-zone soil moisture (Rsvm) across China through a weighted averaging algorithm, as detailed in Appendix B.

Vegetation factors, key drivers of the Evaporative Stress Index (ESI), were incorporated through the Global Land Surface Satellite (GLASS) datasets, including the Leaf Area Index (LAI; version V60) and Fractional Vegetation Cover (FVC; version V40) (Liang et al., 2013, 2014, 2021). These datasets feature an 8 d temporal resolution and 0.05° spatial resolution spanning 2001–2020. Additionally, Moderate Resolution Imaging Spectroradiometer (MODIS) datasets, including the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) (MOD13C1) with a 16 d temporal resolution, were employed. Besides, to mitigate saturation effects and enhance robustness to noise, kernel density estimation was applied to smooth the NDVI data, creating kernel NDVI (kNDVI) (Camps-Valls et al., 2021) (Eq. 2).

$$\begin{array}{}\text{(2)}& \text{kNDVI}=\tanh \left({\text{NDVI}}^{\mathrm{2}}\right)\end{array}$$

Grids with NDVI values below 0.1 are omitted to eliminate barren, rocky, sandy or snow areas. The Aridity Index (AI) is calculated as a multiyear average (2001–2020) by the ratio of annual Prec to ET_p using the ERA5 dataset. Lower AI values indicate drier regions, classifying the study area into sub-humid (0.5 ≤ AI < 0.65), semi-arid (0.2 ≤ AI < 0.5), and arid (0.05 ≤ AI < 0.2), with regions exceeding an AI of 0.65 being categorized as non-dryland areas (Spinoni et al., 2015). To mitigate the effects of land-use change, the annual China Land Cover Dataset (Yang and Huang, 2021), featuring a 30 m grid spacing for the years 2001–2020, is employed. Only grids with consistent land cover types (cropland, forest, or grassland) throughout the study duration were selected for analysis, as depicted in Fig. A1.

The datasets and pre-processing procedures are summarized in Table S1 in the Supplement. All variables were interpolated to a spatial resolution of 0.1° to match the ERA5 dataset and were composited over three distinct temporal scales: 8 d, 16 d, and monthly. The selection of datasets was guided by criteria that included high accuracy, fine spatial resolution, and a broad temporal coverage.

2.2 Identification of the key drivers

Hydrothermal conditions and vegetation variables can each be represented by multiple proxies. For vegetation health, indicators like FVC, LAI, NDVI, and kNDVI are commonly used. In contrast, the energy conditions considered in this study are limited to T and Rn. A disparity in the number of variables may introduce bias, potentially favouring some feature groups over others, which may distort feature importance assessments (Liu et al., 2024). More critically, a high degree of correlation between features – known as multicollinearity – may cause overfitting, where the model becomes excessively tailored to the training data, compromising its ability to generalize.

Figure 2 Illustration of time-series decomposition and procedures for sensitivity analysis. (a) Conceptual diagram of the time-series structure, where the blue line denotes the intra-annual sequence, the green line signifies the interannual trend, and the brown line illustrates the seasonal cycle pattern. (b) The output of the MDLM. (c) The results from the 5-year moving window multiple linear regression model, in which Y represents the anomaly of ESI, X denotes the driving factors, θ symbolizes the sensitivity of Y to changes in X, and the subscript denotes the specific time point.

To identify the principal factors influencing ET stress, we use Explainable machine learning (SHapley Additive exPlanations), specifically XGBoost and Random Forest (RF) regressors. These models provide a robust approach for feature selection, with SHAP values used to evaluate the importance of each driver, where higher SHAP values indicate greater influence (Lundberg and Lee, 2017). Both XGBoost and RF reduce overfitting risks due to their ensemble nature (Breiman, 2001; Chen and Guestrin, 2016). To optimize the regression model, we use grid search techniques and 5-fold cross-validation for parameter refinement. Furthermore, regularization and the “early stopping” mechanism help manage overfitting and control model complexity. Table S2 in the Supplement presents the performance metrics of two models. Despite the above efforts to prevent overfitting, the relatively low validation accuracy still reveals the SHAP model's susceptibility to misallocating feature importance among highly correlated variables. To address this, we generated multiple training subsets through categorical divisions of temporal scales (8 d, 16 d, monthly) and spatial partitions (drought gradient and land cover types: cropland, forest, grassland). This approach enabled us to obtain diverse importance rankings across China's regions and derive statistically robust importance hierarchies through distribution analysis (Figs. 4 and S3 in the Supplement).

Moreover, we employ partial dependence plots (PDPs) to elucidate the effects of key drivers on the ESI, which show how the model's predictions change in response to variations in individual drivers, while keeping other variables constant (Friedman, 2001; Štrumbelj and Kononenko, 2014). To enhance the robustness, we performed 30 iterations with varying random seeds and averaged the results.

Persistent common trends and analogous seasonal patterns can result in overfitting, thereby concealing the actual causal relationships (Li et al., 2022a) (Fig. 2a). In our SHAP model, we input deseasonalized and detrended anomalies of each variable. Seasonality is mitigated by subtracting the climatological mean, and long-term trends are eliminated through the application of a locally weighted smoothing filter (LOWESS) with a span of 0.3 (Cleveland, 1979).

2.3 Exploration of the time-varying sensitivity

The sensitivity of evaporative stress to key drivers (identified in Sect. 2.2) proxies the ESI's response to driver variations. Given self-amplify mechanisms, the previous state of evaporative stress is modelled by a lag-5 autocorrelation. The mathematical expression for sensitivity is as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathit{\delta }{\text{ESI}}_t=& {\mathit{\theta }}_{\text{ESI1}}\mathit{\delta }{\text{ESI}}_{t-\mathrm{1}}\sim {\mathit{\theta }}_{\text{ESI5}}\mathit{\delta }{\text{ESI}}_{t-\mathrm{5}}+{\mathit{\theta }}_{\text{Hydr}}\mathit{\delta }{\text{Hydr}}_t\\ & +{\mathit{\theta }}_{\text{Ener}}\mathit{\delta }{\text{Ener}}_t+{\mathit{\theta }}_{\text{Clim}}\mathit{\delta }{\text{Clim}}_t+{\mathit{\theta }}_{\text{Vege}}\mathit{\delta }{\text{Vege}}_t+\mathit{\epsilon },\end{array}\end{array}$$

where δHydr_t, δEner_t, δClim_t, and δVege_t are anomalies in hydrological, energy, climatic, and vegetation conditions, respectively, with subscripts indicating that the anomalies are calculated from the current time (t). The corresponding θ represents the sensitivity of evaporative stress investigated in this study, and ε is the error term. ${\mathit{\theta }}_{\text{ESI1}}\mathit{\delta }{\text{ESI}}_{t-\mathrm{1}}\sim {\mathit{\theta }}_{\text{ESI5}}\mathit{\delta }{\text{ESI}}_{t-\mathrm{5}}$ denote the stress states from the previous 1 to 5 stages. Considering the lagging effect, we constructed a memory dynamic linear model based on Bayesian forward filtering (MDLM) to analyse the time-varying sensitivity of ESI. This model draws upon the foundational work of Liu et al. (2019) and Zhang et al. (2021b), which primarily examines the vegetation's response to antecedent growth states and prevailing climatic conditions.

Dynamic linear models serve as a sophisticated statistical methodology for modelling time-series signals by incorporating various contributing factors (Prado et al., 2021). For each pixel, we developed a MLDM that integrates long-term evaporation stress trends, seasonal cycles, the effects of internal memory from the previous five stages, and external drivers. DLM decomposes the ESI series into three components: local mean and trend, seasonal pattern, and regression part. The MDLM, utilizing Bayesian forward filtering, captures the dynamic relationship between ESI's regression components and driving factors at each time step, providing an “instantaneous” response. Additionally, this approach allows sensitivity analysis in both stationary and non-stationary scenarios (Liu et al., 2019; Zhang et al., 2021b), with model precision increasing as the sample size grows. Therefore, we calibrated MLDM using an 8 d resolution, the finest temporal scale available, with data allocation comprising 80 % for training and 20 % for validation. Further details are available in Text S1 in the Supplement.

We also employed multiple linear regression (MLR) with a 5-year moving window to calculate a rough sensitivity sequence for comparison. To remove seasonality and isolate the relationship between ESI anomalies and drivers, we subtracted the long-term average over an 8 d interval and applied the LOWESS method to eliminate local trends. Using Eq. (3), we established a MLR model and analysed the regression coefficients (sensitivity coefficients) for each window. It is crucial to acknowledge that the sensitivity derived from the MLR model represents an average response across a designated window (Fig. 2c), providing only a rough sequence of sensitivity to confirm the overarching trend throughout the study period.

The regression coefficients obtained by the two methods are partial regression coefficients, reflecting the magnitude and direction of the relationship between the variable and the ESI while holding other independent variables constant to remove the obfuscating effects of collinearity (Toyoda, 2024). The precise elucidation of real-world phenomena's intricacies necessitates incorporating a broad spectrum of representative drivers. To prevent overfitting due to multicollinearity, we evaluated the variance inflation factor (VIF) for the independent variables before modelling (Belsley, 1991) (Table S4 in Supplement).

2.4 Other statistical analysis

Besides, we employed Sen's slope estimator (Sen, 1968) to quantify the temporal trend and the Mann–Kendall (MK) test to ascertain the statistical significance (Kendall, 1949; Mann, 1945). The slope > 0 signifies a positive trend, and a slope < 0 denotes a negative trend over time.

This study comprehensively investigates the dynamic trajectory of evaporative stress by comparing the ESI calculated using the conventional PM formula from 1950 to 2020 with the ESI derived from the PM formula incorporating CO₂ water-saving effects from 2001 to 2020, enabling a detailed examination of its variation characteristics. Subsequent analyses are conducted using the latter due to limitations in the influencing factors' data series. The SHAP model is applied at 8 d, 16 d, and monthly scales to obtain temporally robust importance rankings. Furthermore, to calculate sensitivity, the MDLM operates at an 8 d temporal resolution, as a larger sample size enhances parameter accuracy. The temporal evolution patterns of sensitivity coefficients are explored at an interannual scale.

Figure 3 Dynamic trajectory of ESI across China from 1950 to 2020. (a) Time series and linear regression of the mean regional ESI, with the traditional PM-equation-derived ESI shown in dark green and the CO₂ water-conservation-adjusted ESI in light green. (b) Intra-annual distribution of ESI, with each line corresponding to a different year and darker hues indicating more recent years. (c) ESI trend over the preceding 7 decades, quantified by Sen's slope, where lighter shaded regions denote statistically significant trends at the 0.05 level. The inset graphically summarizes the significant (darker shades) and non-significant (lighter shades) changes of ESI across various vegetation types. Panel (d) and (e) illustrate the distribution of ESI trends in both arid and humid zones, as well as across diverse land-use categories.

3 Results

3.1 Evapotranspiration stress is intensifying in China

Mainland China demonstrated a pronounced decline in the ESI from 1950 to 2020 (total change of 4.74 %, decreasing by 0.54 % d⁻¹, p < 0.001). While the CO₂ water-saving effect partially offset the decline, the downward trend remains largely parallel to the long-term series, indicating a continued exacerbation of evaporative stress (−0.70 % per decade, p = 0.31). Annual tracking over 70 years enhances statistical significance due to the expanded sample size (Fig. 3a). Collectively, the ESI for each month is on a downward trend, with stability in the early growing season (0.58 ± 0.08). Commencing in April, the onset of the growing season prompts a marked increase in ESI and culminates in a peak between 0.9 and 1.0 by August (Fig. 3b).

Spatially, over one-quarter (28.01 %) of the regions nationwide have registered a significant reduction in ESI (p < 0.05), concentrating in the North China Plain and the Northeast, where the rate surpassing 2.40 % per decade. Concurrently, stress in the Northwest and the middle to lower Yangtze River basin is easing (Fig. 3c). ESI declines more sharply in dryland (total change of 7.63 %, decreasing by 0.91 % per decade, p < 0.001) than in non-dryland lands (total change of 2.92 %, decreasing by 0.32 % per decade, p < 0.01) (Figs. 3d and S2 in the Supplement). This decline is observed universally across various land-use categories, with cropland seeing the greatest reduction and the highest percentage of areas (−1.14 % per decade, p < 0.001; 7.59 × 10⁵ km², 50.24 % of the total cropland area) (Figs. 3c, e, S1, and S2 in the Supplement). Similar dynamics are observed in ESI factoring in the water conservation benefits attributed to CO₂. Even so, the decline in ESI for drought-prone and agriculture-intensified regions over the past 2 decades has outpaced that of the preceding 70 years (−1.63 vs −0.91 % per decade in dryland, −1.60 vs −1.14 % per decade in cropland), highlighting the increasing water stress in these areas.

Figure 4 Identification of the key external drivers of evapotranspiration stress index (ESI). (a) Matrix plot employs bubble size to represent the absolute values of SHAP, with a colour gradient from blue to red denoting negative and positive effects, respectively. Numerical annotations within each bubble denote the rank of importance for the variables, with higher values signifying greater influence. These ranks are obtained from an explainable machine learning utilizing XGBoost regressors. (b) Box plot provides a statistical representation of the variable rankings across varying underlying surface conditions (corresponding to the row it's in), where the black line denotes the mean ranking and the red dot signifies the median ranking. The analysis incorporates vegetation factors such as FVC, EVI, NDVI, kNDVI, and LAI; meteorological factors including atmospheric pressure (P), carbon dioxide (CO₂) concentration, wind speed (U), and VPD; energy factors encompass air temperature (T) and net radiation (Rn); and water supply factors, which account for precipitation (Prec) and the soil volumetric water content (svm) across various soil layer depths, with subscripts specifying the depth range, with 0–7 indicating the soil layer from 0 to 7 cm, and Rsvm representing root-zone soil moisture. Inputs to the model are the anomalies of all aforementioned variables.

Figure 5 Importance and partial dependence of external driving factors influencing ESI. (a) Relative importance of hydrologic, energy, climatic, and vegetation factors in explaining ESI, with each category encompassing multiple indicators. The significance of these indicators is visually encoded by the colour scheme of the bars matching their labels on the left. (b) Pie chart illustrates the proportion of the four categories described in (a). (c–j) Partial dependence graphs correspond to the following variables, in order: P, U, CO₂ concentration, VPD, T, Rn, LAI, and soil moisture content for the top 0–7 cm. The thickened lines indicate the average effects, and the light lines around them indicate 30 random incidents from the datasets. This figure is based on an 8 d data series and generated through a SHAP model employing a XGBoost regressor, which processes the anomalies of the variables as inputs.

3.2 Soil moisture as the primary external driver

Both RF and XGBoost models pinpoint shallow soil moisture (0–7 cm) as the paramount driver, indicating that its scarcity is most likely to trigger evaporative stress, especially in water-limited regions, consistent with Liebig's law of the minimum (Danger et al., 2008; Tang and Riley, 2021). Among vegetation proxies, LAI performs the best, likely attributable to its strong correlation with stomatal conductance (Fig. 4b). Similar importance ranking distributions were observed using RF regressor (Fig. S3b). Comparatively, the XGBoost regressor exhibited enhanced stability across multiple subsets, with minimal divergence between training and validation metrics (Table S2), prompting its selection for subsequent importance quantification analysis. Cross-validation reveals that, across an 8 d interval, water, energy, climate, and vegetation factors on ESI variability decline sequentially, contributing 42 %, 24 %, 21 %, and 11 % to ESI, respectively (Fig. 5b). Within each category, the drivers exerting the greatest influence are svm0–7 (0.19), Rn (0.18), VPD (0.14), and LAI (0.03) (Fig. 5a). Monthly-scale analysis reveals amplified hydrological regulation (water contribution = 0.56), with concurrent suppression of vegetation and energy controls (Figs. 5b, S4, and S5 in the Supplement), demonstrating scale-dependent hierarchy in ecosystem stress controls.

Following collinear variable filtering, we identified svm0–7, Rn, T, VPD, U, P, CO₂, and LAI as critical external drivers of evaporative stress (parameter set R0). Given precipitation's secondary ranking to svm0–7, we established scenario R1 substituting precipitation to examine hydrological driver divergence (Appendix C). SHAP model calibration across temporal scales revealed ESI response patterns for both variable sets. Reduced multicollinearity enhanced model stability (Table S3 in the Supplement). ESI generally exhibits a nonlinear relationship with several critical factors (Fig. 5c–j), decreasing with a higher VPD, T, and Rn but increasing with U, LAI and svm0–7. Similar patterns are observed across temporal resolutions, reinforcing the robustness of our findings (Figs. S4 and S5). Moreover, the non-monotonic interactions between P and CO₂ with ESI highlight the imperative for further investigation into how these factors exacerbate or mitigate evaporative stress across various numerical ranges. However, except for Rn, increasing temporal scale generally decreases the linearity and monotonicity of PDPs (Fig. 5c–j), likely due to data averaging or noise distortion obscuring relationships (Figs. 5, S4, and S5).

3.3 Undergoing dramatic changes in sensitivity of key drivers

To quantify the influence of these variables on the ESI, denoted by sensitivity parameter θ (Table S5 in the Supplement), factorial simulations are conducted using the MDLM and MLR model. The findings demonstrate solid stability and powerful interpretability after successfully passing the multicollinearity test (Table S4).

Figure 6 Sensitivity of the ESI to external drivers and its dynamic trajectory. (a) Distributions of R², RMSE, and MAE for MLR and MLDM models across all grid cells in training and validation sets, with triangles marking median positions. (b–i) θ time series derived from MLDM (red) and MLR (blue). Dashed lines indicate linear or quadratic fit trends, with shaded regions representing 95 % confidence intervals of the fits. Text annotations display Sen's slope estimates with statistical significance, along with R² values and fit significance; quadratic fits were applied when linear trends were insignificant. Box plots illustrate θ distributions: boxes denote interquartile ranges (25th–75th percentiles), horizontal lines represent medians, dots indicate means, whiskers span 5th–95th percentiles, and text annotations specify t test significance of differences between model results.

Soil moisture stands out as the paramount and most sensitive factor, demonstrating a significant upward trend in both models, albeit with quantitative differences (Fig. 6b). MDLM results indicate θ_svm increased by 9.49 % over the past 2 decades (0.03 ± 0.01 yr⁻¹, p < 0.001), peaking at a 13 % rise from the previous decade in 2016. In contrast, the MLR model, with a slightly lower mean θ_svm (2.36 vs. 2.98), exhibited a stronger growth trend (0.05 ± 0.01 yr⁻¹, p < 0.001), peaking earlier in 2014. This shift in peak timing is an inherent consequence of the moving average process. Similarly, Tang et al. (2025) reported analogous evolutionary patterns and attributed this phenomenon primarily to climatic drivers. Among other pivotal drivers, θ_VPD and θ_LAI exhibit considerable numerical values and trends (Figs. 6 and S6 in the Supplement). On average, a 1 kPa rise in VPD correlates with an ∼ 16.84 % rise in evaporative stress according to the MDLM, slightly lower than the 19.74 % increase estimated by the MLR. The MLR also predicts a higher annual growth rate (0.37 ± 0.07 vs 0.24 ± 0.05 % yr⁻¹, p < 0.001) (Fig. 6c). Regarding LAI, the MLR projects a greater decline, with a cumulative reduction of 82.16 % (compared to 45.77 % in the MDLM) and an annual rate of −0.19 ± 0.02 % yr⁻¹ (vs.-0.16 ± 0.03 % yr⁻¹ in the MDLM, p < 0.05) (Fig. 6d). This suggests that the perceived decline in vegetation greening efficacy for drought relief, as assessed by conventional MLR methods, might be overstated. The consensus result from Parameter Set R1 corroborates the robustness of this observation (Fig. S6d).

Figure 7 Spatiotemporal heterogeneity in sensitivity of ESI to svm0–7, VPD, and LAI. (a–c) Spatial distributions of temporal trends (estimated by Sen's slope) for θ_svm (sensitivity of ESI to svm0–7), θ_VPD, and θ_LAI, respectively. Orange indicates increasing trends, and blue denotes decreasing trends, with hatched patterns marking significance at p < 0.05 (Mann–Kendall test). Insets display spatial distributions of multi-year mean θ. (d–f) Discrepancies in Sen's slope estimates and multi-year mean values between MLR- and MDLM-derived θ (a–c). Pie charts quantify the proportion of land area exhibiting corresponding trend patterns. Grey indicates non-vegetated areas.

Figure 8 Sensitivity comparisons between MDLM and DLM models across regional land types. (a–c) Distributions of sensitivity parameters θ (ESI to svm, VPD, and LAI) in dryland, non-dryland, forest, cropland, and grassland. Box plots show interquartile ranges (25th–75th percentiles, boxes), medians (central lines), means (dots), and data ranges (5th–95th percentiles, whiskers). Corresponding probability density curves are plotted alongside. Top asterisks indicate t test significance between different groups: three asterisks (^∗) denote p < 0.001, N.S. = non-significant differences.

Spatially, θ_svm is predominantly higher in the southeast and lower in the northwest, while θ_Prec exhibits a complementary pattern (Figs. 7a and S10a in the Supplement). This indicates that in humid regions, ESI is primarily affected by soil moisture regulation, whereas in arid regions, short-term precipitation can substantially alleviate evaporative stress. The specific mechanisms are discussed in Sect. 4.2. Additionally, θ_svm also varies by land-use type, with forest showing the strongest sensitivity and most drastic changes, followed by cropland and grassland (MLDM: 4.27 > 3.60 > 1.17, Fig. 8a and d). A plausible explanation is that regions endowed with optimal water conditions and robust ecological structure exhibit increased sensitivity due to the magnified impacts resulting from the disequilibrium of hydrology, climate, and vegetation (Forzieri et al., 2022). Notably, sensitive hotspots were distributed near the “Hu Line” and in the lower Yangtze River basin (Fig. 7a). Given the recent shifts in forest composition (large-scale planting of single tree species) coupled with the increasing extreme climatic events (Ruan et al., 2022; Yin et al., 2022), these areas warrant special attention. Furthermore, sizeable areas (MDLM: 4.00 × 10⁶ km², approximately 48.3 % of China mainland) experienced a rise in θ_svm (Sen's slope > 0, Fig. 7a), with a similar pattern observed in MLR (Fig. S9a in the Supplement). However, MLR underestimates θ_svm in nearly all vegetated areas (∼ 93.6 %), while markedly overestimating its increasing trend in ∼ 58.5 % of these areas, particularly in the south (Figs. 7d and S9a).

Similarly, MLR shows a comparable spatial pattern of underestimation in values and overestimation in trends when estimating θ_VPD (Figs. 7b, e and S9b). It exhibits notable regional disparities, with dryland encountering a 1.64-fold greater increase in ET stress per unit of VPD compared to wet regions (MDLM: −23.02 % vs −14.07 %). θ_VPD is decreasing in arid regions (−0.89 % per decade, p < 0.01), in contrast to the increasing trend in humid regions (3.36 % per decade, p < 0.001); therefore the discrepancy between dryland and non-dryland is anticipated to widen (Fig. 8b and e). Notably, MLR indicates an even greater trend divergence between dryland and non-dryland regions (−2.41 % vs. 6.58 % per decade, p < 0.001, Fig. 8e). Interestingly, the response of evaporative stress to LAI (θ_LAI) exhibits in both positive and negative ways. Approximately 29.6 % (2.81 × 10⁶ km²), counterintuitively, experienced an increase in LAI coupled with intensified ET stress (θ_LAI < 0), mainly in the Loess Plateau, North China Plain, and lower Yangtze River region (Fig. 7c). In these areas, rapid greening and intensive agricultural expansion are linked to heightened water deficits (Feng et al., 2016; Lu et al., 2018; Wang et al., 2024). In the lower Yangtze River region, increased LAI is associated with reduced evaporation (Zhang et al., 2021c), possibly due to shading effects that lower surface radiation, soil temperature, and soil evaporation (Forzieri et al., 2020). Regions where increased LAI alleviated stress are predominantly near the “Hu Line” (θ_LAI > 0, 4.01 × 106 km², ∼ 42.27 %), where sensitivity trends showing a significant decline (Sen's slope < 0, Fig. 7c). This region is a focal point for China's afforestation efforts (Zhan et al., 2023), implying that the benefits of sustained greening in reducing ET stress might be lessening. Additionally, artificial agronomic regulation renders minimal and stable θ_LAI in cropland (MLDM: θ_LAI = 1.11 %, Sen's slope = −0.48 % per decade, Fig. 8c and f). Figure 8 further illustrates that the values and trends of θ_svm, θ_VPD, and θ_LAI, as estimated by MLDM and MLR, display statistically significant differences (p < 0.001) across most climate zones and vegetation types, underscoring substantial methodological discrepancies. Within the same land-use categories, the MLR model consistently demonstrates more dispersed sensitivity trend distributions and higher intra-group variability (Fig. 8).

Figure 9 Sensitivity and its trend of ESI to svm0–7, VPD, and LAI along greening and aridity gradient. (a) The grouping statistics of θ_svm across different bins, categorized by differing degrees of aridity and trends in LAI. Median value is shown for each grid and the number of grid points within each bin is shown in Fig. S12 in the Supplement. Similar to (a), (b) pertains to the grouping statistics of the trends of θ_svm. (c) The conceptual diagram delineates the covariation between LAI trends and θ_svm, including their respective trends, in regions with an Aridity Index (AI) below 1.5. It offers a conceptual interpretation of the patterns depicted in (a) and (b) for these areas. Note that an increase in ESI indicates a mitigation of evaporative stress. Panels (d–f) and (g–h) are analogous to (a–c); however, their main focus is VPD and LAI, respectively. In panels (c–i), the upper right corner indicates the aridity of regions adhering to the specified pattern. A red line with a “+” sign denotes direct proportionality, whereas a “−” sign accompanied by a blue line indicates inverse proportionality. A grey dashed line represents an ambiguous relationship, and a black dashed line highlights the sensitivity parameter θ and its corresponding trend.

3.4 Sensitivity of soil moisture increases with enhanced greening

Given the prominent greening tendency in China, we further explored the connection between greening and the sensitivity of dramatically changing variables (Fig. 9). In extremely humid regions (AI > 1.5), faster greening trends are associated with lower θ_svm values (Fig. 9a). Conversely, in less humid regions (AI < 1.5), particularly in semi-arid and semi-humid zones (0.2 < AI < 0.65), both θ_svm and its trend increase alongside the gradient of greening trends, suggesting higher sensitivity and greater acceleration (Fig. 9a and b). Specifically, when 0.2 < AI < 0.5, θ_svm rises progressively from 1.30 to 2.20, while the associated trend simultaneously escalates from ∼ 0 to 0.2.

Unlike θ_svm, θ_VPD and θ_LAI (and their trends) do not exhibit a co-varying gradient pattern along the greening trends (Fig. 9d–i). Remarkably, drier regions with slower pace of greening display greater absolute θ_VPD values (|θ_VPD|, Fig. 9d), suggesting that heightened VPD amplifies ET stress per unit increase. However, in rapid greening areas, |θ_VPD| tends to be diminished (Fig. 9d), likely because denser vegetation enhances transpiration, increasing atmospheric moisture and partially offsetting moisture deficit stress (Fig. 9d and f). These findings imply that vegetation greening may, to some extent, improve the ecosystem's resilience to atmospheric drought. Concerning θ_LAI, a faster LAI increase corresponds to a steeper decline in θ_LAI, reducing the efficacy of LAI in mitigating ET pressure (Fig. 9h and i). Consequently, greening-induced water stress mitigation may be less pronounced than anticipated.

4 Discussion

4.1 Mechanistic advancements of the MDLM

Previous research has largely focused on the long-term average response of individual variables in the hydrological–climatic–vegetation nexus. Our study employs the ESI as a metric to assess the dynamic equilibrium of the system and introduce the MDLM to capture its time-varying response (Fig. 1). The MDLM, characterized by its mechanism-driven architecture and operational simplicity, can inherently detrend and deseasonalize data, eliminating the need for anomaly pre-processing. It shows a 44.26 % improvement in simulation accuracy over the moving-window MLR (R² increasing from 0.61 to 0.88, Figs. 6a, and S6–S8), consistent with high fitting accuracies (> 0.90) reported in other studies (Zhang et al., 2021b). This improvement is primarily attributed to Bayesian forward filtering embedded in the MDLM for dynamic parameter estimation and its capability to approximate nonlinear dependencies as quasi-linear relationships at discrete temporal nodes (Fig. 3). Consequently, it provides time-specific analyses, unlike the window-averaged sensitivity produced by MLR (Liu et al., 2019; Zhang et al., 2022). Pixel-scale model comparisons indicate that the MDLM outperforms the MLR across the entire study area (Figs. S7 and S8 in the Supplement), demonstrating its high flexibility and broad applicability to diverse climatic and terrain conditions. More importantly, comparing the θ sequences for key variables (svm, VPD, and LAI) reveals that the MLR significantly overestimates their temporal trends (Figs. 6–8, S6, S9, and S10). These MLR processes – truncating sequence, smoothing, and simplifying nonlinear dynamics – not only lead to signal loss but also cause substantial differences in θ evolution trajectories (Figs. 6g, h and S6g). These findings highlight the need to re-evaluate the suitability of traditional moving average methods for determining inter-variable coupling strength or temporal sensitivity trends.

4.2 Enhanced evapotranspiration stress and its moisture-driven mechanisms

Our findings indicate that despite the water-saving benefits of elevated atmospheric CO₂ levels, various regions in China have distinctly experienced increased evaporative stress (Figs. 3a, b and S1). This phenomenon is in concordance with several recent detrimental environmental shifts, including land surface desiccation, rising evaporation demand, and groundwater depletion (Jasechko et al., 2024; Qing et al., 2023; Yuan et al., 2019). The pronounced increase in ET stress in dryland and cropland underscores the vulnerability of these ecosystems (Figs. 3c–e, S1, and S2a, c), primarily attributed to the considerable desiccation of soil moisture (Fig. S11a in the Supplement). Soil moisture, as the exclusive direct water supply source for terrestrial evaporation processes – whether from land surface or plant transpiration – exerts a paramount role (Smith and Boers, 2023; Zhao et al., 2023) (Figs. 4–6 and S3–S6). Over the past 20 years, increased soil water sensitivity and variability has led to greater ET stress intensification (Figs. 6b, 7a, and S11a). Sporadic studies on time-varying sensitivity also reveal that water supply increasingly impacts ecosystems (Hu et al., 2023; Li et al., 2022a; Zeng et al., 2022; Zhang et al., 2022). Our simulations, using precipitation as a water proxy for moisture, also corroborated these findings (Figs. S6b and S10). Hence, addressing the threats posed by soil moisture deficits in arid cropland is crucial for ensuring water and food security.

Intriguingly, our research identifies a complementary spatial pattern between θ_svm (southern high–northern low) and θ_Prec across China (Figs. 7a and S10a). In northern arid/semi-arid regions with low initial soil moisture, precipitation acts as an “instantaneous input” that rapidly replenishes soil moisture and substantially alleviates evaporative stress, exhibiting higher sensitivity. Conversely, in southern humid regions where soils maintain near saturation, surpassing the pre-existing hydro-equilibrium threshold triggers disproportionately severe ecosystem perturbations. Consequently, ESI in these regions demonstrates heightened dependence and sensitivity to soil moisture dynamics. Analogously, despite a relatively moderate increase in forest ET stress (Figs. 3c–e and S2), soil moisture sensitivity emerges both highly pronounced and rapidly enhanced recently (Fig. 8a). In regions with favourable hydrological conditions and stable ecological structures, the collapse of the hydrological–climatic–vegetation system can precipitate profoundly adverse effects (Forzieri et al., 2022). Furthermore, in the face of increasingly intense and severe extreme climatic events, complex ecosystems require closer monitoring to formulate proactive strategies.

4.3 Ecological implications of “greening but drying” feedbacks in a changing climate

Additionally, the pace of greening may influence ecosystems' ability to withstand different types of water stress. On the one hand, in semi-arid and semi-humid regions, the sensitivity of ESI to soil moisture and its variability are proportional to the greening trend (Fig. 9a–c). Therefore, once a soil drought occurs, the adverse impact on ET stress increases, suggesting that overly rapid greening can diminish an ecosystem's resilience to soil drought. Specifically, when afforestation surpasses the carrying capacity dictated by regional hydroclimatic conditions, a slight soil deficit can trigger the self-amplify cycle and result in vegetation degradation and mortality, known as overshoot drought (Zhang et al., 2021b). From a long-term lens, reliance on irrigation during the early stages of artificial vegetation establishment, aimed to ensure plant survival, can hamper the formation of deep root systems (Moreno-Mateos et al., 2020; Xiao et al., 2024b). This practice weakens the plant–groundwater linkage in drylands, subsequently diminishing their resilience (Wang et al., 2023). On the other hand, VPD sensitivity and the greening trend exhibit an inverse gradient (Fig. 9d), particularly in arid lands with high θ_VPD values Fig. 8b), indicating that greening in drylands has improved the ability to cope with atmospheric drought. Two mechanisms may explain this phenomenon:

1. In certain arid zones, deep-rooted vegetation accesses deep soil moisture during droughts, releasing it into the atmosphere via transpiration, which increases humidity but may also worsen soil aridity (Liang et al., 2024; Sun et al., 2021).

2. Greening reduces surface albedo, limiting surface temperature rise and mitigating atmospheric drought through indirect thermal regulation (Zhang et al., 2024).

Concerning LAI, a faster greening trend correlates with a diminishing contribution of LAI growth in alleviating water stress (Figs. 9d, 8h, and S13h in the Supplement). In certain regions, enhanced LAI may exacerbate water stress (Fig. 7c). Gleason et al. (2017) found that high vegetation density may exacerbate inter-species competition for water, leading to increased water stress for individual plants or specific communities, providing a cogent explanation for our findings. Collectively, although greening has traditionally been considered a positive environmental adjustment, we advocate for a more critical and dialectical understanding. In adherence to sustainable development imperatives, policymakers and decision-makers are tasked with the precise identification of regions amenable to greening initiatives and the establishment of quantifiable benchmarks that will safeguard ecological equilibrium and foster enduring sustainability.

4.4 Research limitations and prospects

Although this study provides valuable insights into the interplay between regional greening and complex ecological couplings in China, it remains subject to certain limitations in methodology and data. First, the observational sample size from flux stations is insufficient, as it has been employed solely for the screening of reanalysis products rather than dynamic factorial simulations. Second, deviations in Parameter Set R0 and R1 indicate that other local factors such as soil texture, vegetation acclimation, plant demographic rates, and vegetation–groundwater dependency may exert considerable influence (Abel et al., 2024; Fu et al., 2024b; Patel et al., 2021). More in-depth studies based on long-term ecophysiological observations and purpose-built field experiments are needed to further unravel the complex mechanisms of VPD, soil moisture, and LAI in land–atmosphere interactions.

Finally, what are the key implications given the glimpse of intensifying evapotranspiration stress coinciding with its increasing sensitivity to soil water deficits in a greening China? Soil water variability would lead to greater changes in the ecosystems undergoing rapid greening. That is, stronger drought effects can be expected when soil water is anomalously low. Concurrently, the increase in LAI may not alleviate evaporative stress as much as expected. Combined with the increasing frequency and severity of extreme climate events, a “greening but drying” trend may thus be more prevalent in the future, potentially raising the risk of ecosystem imbalance.

5 Conclusions

Our study introduces a memory-driven dynamic linear model (MDLM) that integrates the “dry gets drier” legacy effects to evaluate time-specific sensitivities of evapotranspiration stress drivers in China under rapid vegetation greening. Using the Evapotranspiration Stress Index (ESI) as a proxy for ecosystem water–atmosphere–vegetation equilibrium, we documented a 4.74 % intensification of ET stress across mainland China from 1950 to 2020. Soil moisture sensitivity was the dominant driver, rising by 9.49 % during 2001–2020. Enhanced vegetation greening exhibited stronger coupling with elevated soil moisture sensitivity but reduced VPD sensitivity, reflecting greater susceptibility to soil drought but enhanced resilience to atmospheric aridity. Compared to the conventional moving-window multiple linear regression (MLR) method, the MDLM framework increased the coefficient of determination (R²) by 44.26 %, offering more accurate sensitivity estimates. Methodologically, this work underscores the need to critically re-evaluate trend overestimation inherent in traditional MLR methods. Our findings advance mechanistic understanding of complex regional ecosystem dynamics and offer guidance for steering greening strategies toward a more stable equilibrium among water, atmosphere, and vegetation.

Appendix A: Calculation and validation of the ESI

We utilize the ratio of actual evapotranspiration (ET_a) to potential evapotranspiration (ET_p) to calculate the Evapotranspiration Stress Index (ESI). Before investigating the dynamic changes of ESI, we compared three widely used datasets: the Global Land Evaporation Amsterdam Model v3.7a (daily, 0.25°), the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) MOD16A2GF Version 6.1 (8 d, 500 m), and the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5-Land (hourly, 0.1°). The results are subsequently validated against observational data from 26 eddy covariance (EC) stations (Fig. A1), sourced from various observation networks (station details in Table S1 in the Supplement). Each site has at least 1 year of continuous observations and energy balance residuals < 35 W m⁻² (calculated using Eq. A1). Outliers are removed using a 3-standard-deviation method, and rainy day observations are excluded. Subsequently, data imputation is conducted using the IterativeImputer tool in Python. The latent heat measured at the EC station is converted to ET_a (Eq. A2), and ET_p is calculated using the traditional Penman–Monteith method (Eq. A3).

$$\begin{array}{}\text{(A1)}& \text{Residual}=R_{\mathrm{n}}-G-\text{LE}-H,\end{array}$$

where R_n represents net radiation; G denotes soil heat flux, which is zero on a daily scale; LE refers to daily latent heat; and H represents sensible heat flux, with all units expressed in W m⁻².

$$\begin{array}{}\text{(A2)}& {\mathrm{ET}}_{\mathrm{a}}={\displaystyle \frac{\text{LE}}{\mathit{\lambda }}}\end{array}$$

In Eq. (A2), the unit of ET_a is mm d⁻¹, with the constant value of λ at 0 °C taken as 28.94.

$$\begin{array}{}\text{(A3)}& {\mathrm{ET}}_{\mathrm{p}}={\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}U(e_{\mathrm{s}}-e_{\mathrm{a}})}{\mathrm{\Delta }+\mathit{\gamma }(\mathrm{1}+\mathrm{0.34}U)}},\end{array}$$

where the unit of ET_p is mm d⁻¹; Δ is the slope of the saturation vapour pressure versus temperature curve (kPa °C⁻¹); γ is the psychrometric constant (kPa °C⁻¹); U is wind speed at 2 m (m s⁻¹); T is surface air temperature (°C); and e_s and e_a are saturated and actual vapour pressure (kPa) and whose difference serves as a crucial indicator of atmospheric moisture deficit, known as vapour pressure deficit (VPD). At the grid scale, all climatic variables used to calculate ET_p, including Rn, T, and U, are obtained from the ERA5 dataset. VPD is calculated based on the Clausius–Clapeyron relation using 2 m temperature and dew point temperature (Held and Soden, 2006; Zhong et al., 2023).

To avoid overestimating ET_p, we employed the formula developed by Yang et al. (2019) to calculate grid ET_p, accounting for the CO₂ water-saving effect (Eq. A4). This formula has been widely recognized and applied in recent studies, providing a solid theoretical basis for accurate evapotranspiration estimates (Lian et al., 2021; Zhang et al., 2023). All meteorological data used in the formula came from ERA5, while monthly CO₂ data from CarbonTracker (CT2022) covered 2001–2020 at a 3° × 2° resolution. After considering the water-saving effect of CO₂, the mean, maximum, and minimum values of ET_p in China decreased by 9.48, 30.85, and 19.42 mm, respectively, while the pattern is mirrored (Fig. A2a, b, and d). The conventional ET_p calculation method overestimated values across 99.6 % of mainland China, particularly in northern arid regions (Fig. A2c).

$$\begin{array}{}\text{(A4)}& {\mathrm{ET}}_{\mathrm{p}}={\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}U(e_{\mathrm{s}}-e_{\mathrm{a}})}{\mathrm{\Delta }+\mathit{\gamma }\mathit{\{}\mathrm{1}+U[\mathrm{0.34}+\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{4}}(\left[{\mathrm{CO}}_{\mathrm{2}}\right]-\mathrm{300})\left]\mathit{\right\}}}},\end{array}$$

where 2.4 × 10${}^{-\mathrm{4}}(\left[{\mathrm{CO}}_{\mathrm{2}}\right]-\mathrm{300})$ denotes the effect of atmospheric CO₂ concentration (ppm) on surface stomatal resistance.

Figure A3 compares the ESI derived from three datasets, revealing substantial discrepancies attributable to methodological differences. For ET_p estimation, MODIS and ERA5 apply the PM formula calibrated for idealized vegetation types, consistent with our research hypotheses in Fig. 1. In contrast, GLEAM integrates actual vegetation heterogeneity. The ET_a < ET_p constraint in MODIS and GLEAM limits their ability to capture water surplus conditions (Fig. A3c and d), whereas ERA5 – free of this restriction – produces ESI spatial patterns aligned with eddy covariance (EC) flux measurements (Fig. A3a and b). Incorporating CO₂ water-saving effects into ET_p calculations identifies an expansion of non-water-limited regions, enhancing aridity-humidity characterization ontological congruence (Fig. A3a, b, and e). Note that our research is dedicated to analysing the temporal dynamics of the ESI, particularly in relation to the intensification of drought conditions. The conceptual diagram (Fig. 1) demonstrates that lower ESI values signify increasingly severe imbalances and stress levels. The study does not address how often the threshold of 1 is crossed, and it does not explore the physical significance of this threshold.

Figure A1 Location of the study area and distribution of multi-year average ESI. The extent of the study area, employing a tri-chromatic map to differentiate between cropland, forest, and grassland; dots denote the positions of EC flux stations, with colours representing the International Geosphere-Biosphere Programme (IGBP) classifications and shapes indicating observational sources. The inset displays time-series plots of regional average air temperature (red), surface soil moisture (blue), and LAI (green) over the past 20 years, with a dashed black line denoting the linear trend. The red slanted line denotes the position of the “Hu Line” (also known as the Hu Huanyong Line), with its neighbouring transitional area defined as the “Hu Zone”, fundamental geographical demarcation in China that characteristically differentiates regional hydrological regimes, vegetation coverage gradients, and demographic distribution patterns (Hu, 1935; Li et al., 2024).

Figure A2 Comparison of the spatial distribution of the multi-year average potential evapotranspiration (ET_p). (a) Calculated by the traditional PM equation and (b) the modified PM equation that incorporates the CO₂ water-saving effect. The maximum value, minimum value, and average value are marked in the lower left corner of the figure. (c) The spatial distribution of the difference between ET_p and ET_pPMCO₂. (d) The probability density curves of the two.

Figure A3 Spatial distribution of the multi-year average of ESI calculated from ET_a and ET_p products from different datasets. (a) was calculated using flux tower observational data, (b) derived from ERA5, (c) originated from GLEAM, (d) obtained from MODIS, and (e) sourced from ERA5 but incorporating the CO₂ water-saving effect in ET_p calculation.

Appendix B: Spatiotemporal patterns of Rsvm in China

Previous studies typically generalized the 0–100 cm soil layer as the root zone. However, since the root zone is a highly dynamic interface for water–vegetation–atmosphere interactions, employing more realistic root zone soil moisture data as an indicator for evapotranspiration stress is of critical importance (Gao et al., 2014, 2024; Wang-Erlandsson et al., 2016). This research utilizes root depth data based on an 80-year drought recurrence interval, estimated by Stocker et al. (2023) through the mass accumulation curve methodology. We calculated actual root zone soil moisture (Rsvm) using the weighted averaging method with ERA5 layered soil moisture data (2001–2020). For cases where root depth exceeded 289 cm (the maximum depth covered by ERA5 data), the weighted average soil water content from 0–289 cm was adopted as the Rsvm for this pixel. The computational procedure is outlined below.

The ERA5 soil moisture dataset is stratified into four layers (Layer 1–4), corresponding to soil depths of 0–7, 7–28, 28–100, and 100–289 cm. For each pixel, the contribution of each ERA5 layer to the root zone was calculated based on root depth parameters.

Layer 1 (0–7 cm):

$$\begin{array}{}\text{(B1)}& d_{\mathrm{1}}=min(\mathrm{7},\text{root\hspace{0.17em}\hspace{0.17em}depth}).\end{array}$$

Layer 2 (7–28 cm), activated only when root depth > 7 cm:

$$\begin{array}{}\text{(B2)}& d_{\mathrm{2}}=max[\mathrm{0},min(\mathrm{28},\text{root\hspace{0.17em}\hspace{0.17em}depth})-\mathrm{7}].\end{array}$$

Layer 3 (28–100 cm), activated only when root depth > 28 cm:

$$\begin{array}{}\text{(B3)}& d_{\mathrm{3}}=max[\mathrm{0},min(\mathrm{100},\text{root\hspace{0.17em}\hspace{0.17em}depth})-\mathrm{28}].\end{array}$$

Layer 4 (100–289 cm), activated only when root depth > 100 cm:

$$\begin{array}{}\text{(B4)}& d_{\mathrm{2}}=max[\mathrm{0},min(\mathrm{289},\text{root\hspace{0.17em}\hspace{0.17em}depth})-\mathrm{100}].\end{array}$$

The Rsvm is

$$\begin{array}{}\text{(B5)}& \mathrm{\text{Rsvm}}={\displaystyle \frac{d_{\mathrm{1}}\times {\text{svm}}_{\mathrm{1}}+d_{\mathrm{2}}\times {\text{svm}}_{\mathrm{2}}+d_{\mathrm{3}}\times {\text{svm}}_{\mathrm{3}}+d_{\mathrm{4}}\times {\text{svm}}_{\mathrm{4}}}{\text{root\hspace{0.17em}\hspace{0.17em}depth}}},\end{array}$$

where root  depth denotes actual root depth (cm), d_1–4 represent weighted thicknesses of individual layers (cm), svm_1–4 indicates soil moisture values per layer, and Rsvm corresponds to root zone soil moisture (m³ m⁻³).

Figure B1 investigates the spatiotemporal patterns of root zone soil moisture across China. The regional mean root depth is 2.33 m, with northern arid/semi-arid regions dominated by shallow-rooted vegetation (e.g., grasslands and shrubs) exhibiting lower root zone moisture values (< 0.30 m³ m⁻³). Deeper root systems (> 15 m) occur in northwestern arid zones and southern karst regions, where vegetation taps into groundwater via deep roots (Fig. B1a). The 20-year mean Rsvm is 0.33 m³ m⁻³, displaying a clear south-to-north gradient, which corresponds to precipitation patterns (Fig. B1b–d). Notably, over 50 % of vegetated areas exhibit declining Rsvm trends, particularly in northwestern China and the North China Plain, where rates exceed 0.002 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. B1c). These regions have mean annual temperatures of 18–26 °C and annual precipitation below 900 mm (Fig. B1d). These soil-drying trends likely arise from climatic stressors and unsustainable groundwater extraction, calling for urgent policy interventions to protect agricultural productivity and ecosystem stability. Conversely, northeastern and southeastern China demonstrate rising Rsvm trends, where annual precipitation remains abundant. This hydrological shift requires systematic assessment of potential reductions in flood mitigation capacity, underscoring the need for enhanced monitoring protocols in watershed management to address evolving hydrological risks.

Figure B1 Spatiotemporal patterns of root-zone soil moisture (Rsvm) in China. (a) Spatial distribution of root depth. (b) Multi-year mean Rsvm (2001–2020). The inset illustrates the probability density distribution with red and blue markers denoting the mean and median values, respectively. (c) Sen's slope of Rsvm during 2001–2020, where shaded areas indicate regions with statistically significant trends (p < 0.05). The inset displays areal proportions of Rsvm changes: orange represents decreasing trends, blue increasing trends, and grey indicates no-data regions. (d) Climatological phase space diagram corresponding to (a)–(c), with the x axis showing multi-year mean air temperature and the y axis representing precipitation. Colour gradations within bins reflect median values, while dots mark regions exhibiting statistically significant Rsvm trends (p < 0.05).

Appendix C: Comparative impacts of precipitation and surface soil moisture on ESI

Since precipitation ranked second in importance among parameters (after svm0–7), we developed a comparative Parameter Set R1 using precipitation anomalies to investigate the distinct influences of precipitation and soil moisture on the ESI. SHAP analysis across multiple temporal scales in China revealed contrasting feature importance patterns between parameter sets (see Table S3 for model performance metrics): R1 displayed a broader distributions of contribution indicating more balanced interactions among parameter, while R0 prioritized dominant parameters (svm0–7, Rn). Both configurations consistently showed Rn, VPD, and T negatively regulating ESI, contrasting with positive regulation by LAI and U (Fig. C1). Notably, Prec and svm0–7 synergistically enhanced ESI through moisture supply mechanisms (Figs. C1, C2, 5j, S4j, and S5j), consistently ranking as top predictors. In R1, Prec showed persistent secondary influence to VPD, whereas svm0–7 dominated feature importance at both 16 d and monthly scales (Fig. C1b and c). Unlike the linear relationship of soil moisture, precipitation anomalies induced a biphasic “plateau-ramp” ESI response, requiring anomalies exceeding −5 to 0 mm for stress mitigation. Conversely, svm0–7 exhibited acute sensitivity with near-vertical monthly response gradients (Fig. S5j), indicating immediate drought relief from minimal moisture replenishment (−0.02 to −0.01 m³ m⁻³) during extreme deficits. Our deseasoned-detrended anomaly analysis inherently accommodates negative values. Subsequent research should prioritize identifying precise response thresholds and turning points.

Figure C1 Multi-temporal beeswarm plots between Parameter Set R0 (a–c) and R1 (d–f). Parameter set R0 includes svm0–7, Rn, T, VPD, U, P, CO₂, LAI, with R1 substituting svm0–7 with Prec. Red denotes high feature values with substantial impacts on ESI. High feature values distributed along the SHAP positive half-axis signify positive contributions to ESI. All features are arranged in descending order of contribution magnitude, where wider x-axis distributions reflect stronger model impacts.

Figure C2 Precipitation PDPs under Parameter Set R1 across temporal scales, with soil moisture counterparts shown in Figs. 5j, S4j, and S5j.