Declining Sensitivity and Increasing Resistance Time of Ecosystem Water Use Efficiency to Meteorological Drought

Drought is a dominant factor influencing terrestrial ecosystem water-use efficiency (WUE). However, the coupling relationship between WUE and drought remains insufficiently understood. Currently, the coupling relationship is primarily assessed using correlation coefficients or linear regression slopes. However, the optimal drought timescale at which WUE responds to drought has largely been overlooked. Therefore, this study investigated the spatiotemporal patterns of the WUE – meteorological drought coupling relationship across global terrestrial ecosystems from 1982 to 2018 with satellite- derived and model-simulated WUE, together with the Standardized Precipitation-Evapotranspiration Index (SPEI), and explored the potential causal mechanisms. Within the framework of WUE-SPEI coupling, the maximum correlation coefficient between WUE and SPEI represents the sensitivity of WUE to meteorological drought (R_max), whereas the corresponding optimal drought timescale represents its resistance time (T_opt). The results indicated that the sensitivity of WUE to meteorological drought decreased at a rate of −0.0003 yr⁻¹ (p< 0.01), while the resistance time increased at a rate of 0.0155 month yr⁻¹ (p< 0.01), indicating a weakening of the coupling between WUE and meteorological drought. Attribution analysis indicated that CO₂ fertilization was the primary factor contributing to the weakening of the coupling relationship. Surface soil moisture was the most critical hydrometeorological driver, exhibiting nearly opposite effects and significant threshold effects on R_max and T_opt. Peter & Clark Momentary Conditional Independence (PCMCI+) algorithm was further employed to construct a causality diagnosis framework for identifying the relationships between WUE-drought coupling and temperature, precipitation, radiation, wind speed, vapor pressure deficit, and surface and root-zone soil moisture. The results showed that the decrease in the R_max had direct negative causal effects on precipitation, temperature, and radiation. In contrast, the increase in the T_opt was primarily driven by a negative causal effect of radiation. This study highlights the weakened coupling between WUE and meteorological drought, suggesting that vegetation's carbon-water trade-off is evolving toward drought adaptation, which is crucial for understanding the adaptive strategies of vegetation in response to climate change.

1 Introduction

Ecosystem water use efficiency (WUE) is defined as the ratio of carbon assimilation to ecosystem water evapotranspiration (Beer et al., 2009). High WUE indicates that an ecosystem can more effectively carry out photosynthesis and growth under limited water resources (Wang et al., 2022). Therefore, the WUE is widely used to characterize the trade-off between carbon uptake and water loss within ecosystems (Xue et al., 2022). Meteorological drought, as characterized by the Standardized Precipitation Evapotranspiration Index (SPEI), is one of the most common and destructive natural hazards (Li et al., 2026). The SPEI combines sensitivity to atmospheric evaporative demand, as captured by the Palmer Drought Severity Index (PDSI), with the multi-scalar characteristics of the Standardized Precipitation Index (SPI) (Vicente-Serrano et al., 2010). Therefore, the SPEI is widely regarded as an effective index for capturing the impacts of drought on terrestrial ecosystem carbon-water cycles (Vicente-Serrano et al., 2012, 2013). Against the backdrop of global climate change, the frequency, duration, and intensity of meteorological drought are projected to continue increasing in the future (Mokhtar et al., 2021; Trenberth et al., 2014; Chen et al., 2023), which will have increasingly profound effects on terrestrial ecosystems. Therefore, understanding the impacts of meteorological drought on terrestrial carbon-water dynamics is of great significance for ecosystem sustainability (Wang et al., 2026a).

Previous studies have shown that WUE responses to drought exhibit significant spatiotemporal heterogeneity. Across aridity gradients, WUE exhibits a negative response to drought in arid ecosystems, whereas it shows a positive response in humid ecosystems (Huang et al., 2017). From a temporal perspective of drought development, WUE decreases during summer droughts but increases during autumn droughts (Ma et al., 2019). Furthermore, WUE demonstrates a two-phase relationship with drought intensity: increasing under moderate drought but declining under severe drought (Lu and Zhuang, 2010). The timing of drought occurrence also plays a crucial role in regulating both the direction and magnitude of WUE responses, as ecosystem physiological and physical processes exhibit varying sensitivities to water stress across different growth stages (Huang et al., 2021; Wang et al., 2021).

However, these studies remain limited, as they have largely been conducted at specific or fixed drought timescales, thereby neglecting the cumulative and lagged effects of drought on vegetation dynamics. In fact, plants can uptake water stored from past precipitation in deeper unsaturated soil layers or groundwater to maintain carbon–water balance (Zhang et al., 2026). In other words, drought not only impacts vegetation growth synchronously but also exhibits lagged and cumulative effects, where past drought conditions influence current vegetation growth (Huang et al., 2018; Kannenberg et al., 2020). Vicente-Serrano et al. (2013) emphasized that vegetation responses to drought are inherently multi-timescale in nature. Longer timescales imply stronger memory effects, which may buffer the impacts of recent drought events and thereby reduce vegetation drought sensitivity (Seddon et al., 2016). In contrast, shorter timescales indicate more rapid vegetation responses to moisture deficits, reflecting higher drought sensitivity (Jiao et al., 2021). Therefore, neglecting the role of timescales in assessing WUE-drought coupling may hinder our understanding of how WUE responds to drought. Although increasing attention has been paid to the optimal timescale of ecosystem responses to drought (Guo et al., 2026; Yuan et al., 2024; Xing et al., 2026), most existing studies identify a single fixed optimal drought timescale over the entire study period, thereby primarily revealing its spatial distribution patterns. For example, Yuan et al. (2024) used Spearman correlation analysis to reveal the lagged and legacy effects of meteorological drought on the Normalized Difference Vegetation Index (NDVI) in northern China during 1982–2022. However, the optimal drought timescale is not static; under climate change, its temporal evolution can indirectly indicate changes in ecosystem resistance and sensitivity to drought. Thus, it remains unclear whether the response duration of WUE to drought has lengthened or shortened, whether such temporal changes indicate a weakening or strengthening of WUE-drought coupling, and which environmental or physiological factors dominate these changes. Therefore, incorporating time-varying optimal drought timescales into analyses of WUE-drought coupling is essential for a more comprehensive and accurate understanding of vegetation carbon-water trade-offs under drought stress.

Therefore, we calculated the coupling relationship between WUE and SPEI across time scales ranging from 1 to 24 months, including the maximum correlation coefficient (R_max) and optimal drought timescale (T_opt). Specifically, the R_max represents the sensitivity of WUE to meteorological drought, whereas the T_opt represents its resistance time. Subsequently, we utilized Dynamic Global Vegetation Models from TRENDY project to quantify the contributions of CO₂, climate change (CLI), and land-use change (LCC) to this coupling relationship. To further explore the impact of climate change, we employed the eXtreme Gradient Boosting (XGBoost) algorithm combined with SHapley Additive Explanations (SHAP) to identify both the relative importance and modes of influence of key hydrothermal factors in shaping the coupling relationship between WUE and drought. Finally, we employed Peter-Clark Momentary Conditional Independence Plus (PCMCI+) to uncover the complex causal network between hydrothermal factors, R_max, and T_opt. These findings are expected to enhance our understanding of the vulnerability of terrestrial ecosystems to drought and provides valuable insights for supporting ecosystem sustainability under climate change.

2 Data sources and processing

2.1 GPP and ET data

Three widely used GPP datasets were employed in this study: FLUXCOM GPP, GLASS GPP, and NIRv GPP. FLUXCOM GPP, which is derived from a machine learning-based integration of eddy covariance measurements and remote sensing data, offers global coverage with monthly temporal resolution and a spatial resolution of 0.5°. The GLASS GPP product, generated using a light use efficiency model in conjunction with AVHRR reflectance data, provides 8 d composite estimates with a spatial resolution of 0.5° from 1982–2018. To aggregate the 8 d time-scale GPP data into monthly values, the maximum value method was utilized. NIRv GPP, a recently developed satellite-derived index based on near-infrared reflectance of vegetation, features a monthly temporal resolution and a 0.05° spatial resolution, and exhibits a strong correlation with tower-based GPP measurements. These datasets collectively represent distinct methodological approaches to quantifying terrestrial carbon uptake while maintaining complementary spatiotemporal characteristics suitable for cross-comparison analysis.

Three widely used ET datasets were analyzed in this study: ERA5 ET, GLASS ET, and GLEAM ET. ERA5 ET, produced by the European Centre for Medium-Range Weather Forecasts through atmospheric reanalysis modeling, provides monthly temporal resolution at a 0.25° spatial grid, incorporating land-atmosphere interaction processes. The AVHRR-based GLASS ET data employs the Bayesian Model Averaging method, which merges five process-based ET algorithms to improve ET estimation. It provides ET data spanning from 1981 to 2022 with an 8 d temporal resolution and a 0.5° spatial resolution. In GLEAM, multiplicative evaporative stress factors are applied to convert the estimated potential evapotranspiration values of three land components – bare soil, high-canopy, and short-canopy – into bare soil evaporation (Eb) and transpiration (Et). Interception loss (Ei) is calculated separately using an analytical model driven by precipitation and vegetation characteristics. The actual evapotranspiration is estimated as the sum of these three components. GLEAM ET provides data with a monthly temporal resolution and a 0.25° spatial resolution. These datasets collectively represent diverse retrieval approaches (reanalysis, hybrid model, and observation-driven algorithm) while maintaining spatially and temporally complementary resolutions, enabling robust intercomparison of terrestrial water flux patterns across multiple scales.

2.2 FLUXNET data

To validate the reliability of WUE data observed from multi-source satellite products, this study collected flux data from FLUXNET 2015 (https://fluxnet.org/, last access: 14 January 2025). To ensure the long-term reliability of WUE observations, we selected sites with continuous observation records of more than 5 years and with complete and valid records were selected. Based on this criterion, 85 flux stations were ultimately selected for analysis at the monthly scale (Fig. S1, Table S1 in the Supplement). Among these sites, only high-quality gap-filled data (QC > 0.8) were retained (Guan et al., 2025). Additionally, to ensure the accuracy of ET estimates, the energy balance closure criterion was applied: when the energy imbalance (net radiation minus the sum of latent heat flux, sensible heat flux, and soil heat flux) exceeded 50 % of the available energy, the latent heat flux (LE) was treated as missing (Knauer et al., 2018). GPP was estimated using the nighttime respiration allocation method (“GPP_NT_VUT_REF”) to ensure consistency across sites.

2.3 TRENDY v12 multi-model simulated GPP and ET data

The latest version TRENDY v12 (Trends in Net Land-Atmosphere Carbon Exchange) is a collaborative initiative that integrates multiple Dynamic Global Vegetation Models (DGVMs) (are available on request to Stephen Sitch (s.a.sitch@exeter.ac.uk) and Pierre Friedlingstein (p.friedlingstein@exeter.ac.uk)), aimed at quantifying global carbon budgets using forcing data on carbon dioxide concentrations, climate variables, and land use changes (Friedlingstein et al., 2023). All models utilize the same forcing data, with historical climate fields derived from the CRUv.4.07 and CRU-JRA55 datasets, and global atmospheric CO₂ concentrations obtained from a combination of ice core records and atmospheric observations. Each DGVM model runs four simulation scenarios: S0 (no changes in CO₂, climate, or land cover), S1 (changes in CO₂, static climate conditions, and static land use), S2 (changes in CO₂ and climate conditions, static land use), and S3 (changes in CO₂, climate conditions, and land use).

2.4 SPEI data

The SPEI is commonly used to assess climate drought conditions. By incorporating regional precipitation, potential evapotranspiration, and temperature, the SPEI quantifies the hydrological balance over a specific period and provides an objective method for comparing drought severity across different regions and time periods. It is of considerable importance for monitoring and predicting the impact of climate change on water resources and ecosystems. To ensure consistency with the climate data version used in the TRENDY model (CRU TS 4.07), we selected the SPEI version 2.9 dataset, which is derived from CRU TS 4.07 (https://spei.csic.es/spei_database, last access: 15 January 2025). This dataset has a spatial resolution of 0.5° × 0.5° and covers timescales ranging from 1 to 24 months (Beguería et al., 2014).

2.5 Hydrometeorological factors data

Since CRU climate data (CRUts4.07 and CRU-JRA55) is widely used as climate forcing in DGVMs, we adopt the same climate factor data in this study to ensure consistency with these models. Specifically, we collected monthly mean temperature (Temp), precipitation (Pre), and actual vapor pressure data (V_ap) from the CRUts4.07 climate product, as well as downward shortwave radiation (Rad) and wind speed (WS) data from the CRU Japanese Reanalysis (CRU JAR) (https://crudata.uea.ac.uk/, last access: 15 January 2025). Soil moisture data for the surface layer (0–10 cm, SMsurf) and root zone (10–100 cm, SMroot) were obtained from the Global Land Evaporation Amsterdam Model (GLEAM) (https://www.gleam.eu/, last access: 3 October 2024). Additionally, the vapor pressure deficit (VPD) was calculated based on monthly average temperature and actual vapor pressure using the following equation (Wang et al., 2025):

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}\lg E_w=C_{\mathrm{1}}& \times \left(\mathrm{1}-{\displaystyle \frac{T_{\mathrm{1}}}{T}}\right)+C_{\mathrm{2}}\\ & \times \lg \left({\displaystyle \frac{T}{T_{\mathrm{1}}}}\right)+C_{\mathrm{3}}\\ & \times \left[\mathrm{1}-{\mathrm{10}}^{C_{\mathrm{6}}\times \left(\frac{T}{T_{\mathrm{1}}}-\mathrm{1}\right)}\right]+C_{\mathrm{4}}\\ & \times \left[{\mathrm{10}}^{C_{\mathrm{7}}\times \left(\mathrm{1}-\frac{T_{\mathrm{1}}}{T}\right)}-\mathrm{1}\right]+C_{\mathrm{5}}\end{array}\text{(2)}& {\displaystyle }\mathrm{VPD}=E_w-V_{\mathrm{ap}}\end{array}$$

Where, E_w represents the saturated vapor pressure. and V_ap represents the actual vapor pressure. The constants used are as follows: $C_{\mathrm{1}}=\mathrm{0.1079574}\times {\mathrm{10}}^{\mathrm{2}}$, $C_{\mathrm{2}}=-\mathrm{0.5028}\times \mathrm{10}$, $C_{\mathrm{3}}=\mathrm{0.150475}\times {\mathrm{10}}^{\mathrm{3}}$, $C_{\mathrm{4}}=\mathrm{0.42873}\times {\mathrm{10}}^{-\mathrm{3}}$, C₅=0.78614, $C_{\mathrm{6}}=-\mathrm{0.82969}\times \mathrm{10}$, $C_{\mathrm{7}}=\mathrm{0.476955}\times \mathrm{10}$. Additionally, T₁=273.16 K (the triple point temperature of water), and T= 276.15 +t (K), where t is the air temperature in Celsius.

2.6 Land use types data

This study utilizes the annual HLDA+ Global Land Use Change dataset (https://doi.org/10.1594/PANGAEA.921846, Winkler et al., 2020), which includes six general land use/cover categories: urban areas, croplands, pastures/rangelands, Forest, unmanaged grasslands/shrublands, and sparse/bare soils. We excluded pixels with land use types of “urban areas” and “sparse/bare soils” to define permanent vegetation areas (Winkler et al., 2021), and all subsequent studies were conducted within these permanent vegetation areas (Fig. S4).

2.7 Aridity index data

Using the third edition of the Global Aridity Index and Potential Evapotranspiration dataset (Global-AI-PET-v3, https://doi.org/10.6084/m9.figshare.7504448.v5), the globe is classified into four aridity gradients based on the following thresholds: AI < 0.2 as arid (AR), 0.2 ≤ AI < 0.5 as semi-arid (SAR), 0.5 ≤ AI < 0.65 as dry sub-humid (DSH), and AI > 0.65 as humid (HU) (Fig. S3).

2.8 Pre-processing

All datasets were resampled to a spatial resolution of 0.25° using bilinear interpolation (Li et al., 2022b). Moreover, seasonal cycles and long-term trends in WUE and SPEI data may lead to spurious correlations (Boulton et al., 2022). Therefore, it is essential to perform deseasonalization and detrending on WUE and SPEI prior to subsequent calculations (Smith and Boers, 2023). This study employed the Seasonal and Trend decomposition using Loess (STL) method to decompose WUE and SPEI time series of each grid cell into the overall trend, seasonal component, and residual component, as implemented using the stl() function in the “stats” package in R (v4.2.1). The STL residual component, which represents the deseasonalized and detrended WUE and SPEI time series, were utilized for further analysis (Wang et al., 2023).

3 Methods

3.1 Definition of WUE

The grid-scale ecosystem WUE was calculated as the ratio of GPP to ET according to Eq. (3):

$$\begin{array}{}\text{(3)}& \mathrm{WUE}={\displaystyle \frac{\mathrm{GPP}}{\mathrm{ET}}}\end{array}$$

Where, the unit of WUE is g C m⁻² month⁻¹, the unit of GPP is g C m⁻² month⁻¹, and the unit of ET is mm⁻¹ month⁻¹. In this study, based on these three independent sets of GPP products and three sets of ET products, we first computed WUE for each GPP-ET product combination. Then, we integrated the results of the nine WUE combinations using an arithmetic mean, ultimately generating a multi-product averaged WUE. This approach effectively reduces the impact of systematic errors from individual products, thereby enhancing the robustness of the estimated results (Wu et al., 2025).

Representative long-term scale sites from different land use types were selected to validate the satellite products. In the calculation of site-scale ecosystem WUE, GPP can be directly obtained, while ET is computed using Eq. (4):

$$\begin{array}{}\text{(4)}& \mathrm{ET}={\displaystyle \frac{\mathrm{LE}}{\mathrm{2.501}-\mathrm{2.361}\times {\mathrm{10}}^{-\mathrm{3}}\times T}}\end{array}$$

where ET represents actual evapotranspiration (mm month⁻¹), LE (W m⁻²) represents latent heat flux, and T (°C) is the air temperature. The validation results indicated a high correlation between satellite-based WUE and the flux-based WUE (R= 0.62, Fig. S2), implying that satellite-based WUE results were highly reliable.

3.2 Calculation of coupling relationship between WUE and drought

The coupling relationship was defined as the maximum correlation coefficient (R_max) and optimal lag time (T_opt) between WUE and SPEI. In general, WUE decreases with increasing drought severity, whereas it tends to be higher under mild or non-drought conditions (Li et al., 2025), indicating an overall positive correlation between WUE and SPEI. Accordingly, R_max was defined as the maximum (non-absolute) correlation coefficient between SPEI at different time scales and WUE, representing the maximum sensitivity of WUE to SPEI (Vicente-Serrano et al., 2013). SPEI typically exhibits a lagged effect on WUE, whereby past drought conditions influence current WUE (Huang et al., 2018; Kannenberg et al., 2020). Therefore, T_opt was defined as the SPEI time scale corresponding to R_max, representing the resistance time of WUE to SPEI (Li et al., 2024).

3.3 Trend analysis

An 18-year moving window was employed to investigate the temporal variations in R_max and T_opt over the past four decades. Then, we utilized a combination of Theil-Sen slope estimation and the Mann-Kendall (MK) test to identify and analyze trends in long-term time series data. The Theil-Sen method was applied to calculate robust linear trends, with p-values that remain resilient to the influence of outliers (Gocic and Trajkovic, 2013). At the same time, the non-parametric MK test was employed to assess the significance of monotonic trends by evaluating their slope values (Ma et al., 2020). To ensure the robustness of the identified temporal trends, we further examined the trends of R_max and T_opt using 16-year and 20-year moving windows.

3.4 Attribution analysis and causality diagnosis of WUE-drought coupling relationship

We compared the coupling relationships derived from the simulation results of all DGVMs under the S3 scenario with those calculated from remote sensing observations, and ultimately retained only the seven models that were consistent with the remote sensing results: E3SM, EDv3, JSBACH, JULES, LPJ-GUESS, SDGVM, and VISIT (Zeng et al., 2022). We not only computed WUE for each individual model but also calculated the multi-model mean WUE under the S1, S2, and S3 scenarios (Fig. S6). Furthermore, we derived R_max and T_opt to examine the relationship between WUE and SPEI. Considering that the coupling relationship under the S1 scenario can only represent the effects of CO₂ changes, we also computed the coupling relationships for the (S2–S1) and (S3–S2) scenarios to quantify the contributions of CLI and LCC to the changes in R_max and T_opt.

Machine learning models are capable of effectively capturing the nonlinear relationships between multidimensional predictors and target variables (Yan et al., 2025b). With the advancement of interpretability techniques such as SHapley Additive Explanations (SHAP), machine learning has gradually evolved from a black-box paradigm into an explainable artificial intelligence framework. In this study, “threshold effects” refer to nonlinear transition points at which the marginal contribution of a driving factor to R_max or T_opt changes markedly. Specifically, a threshold was interpreted as the point where the SHAP value showed a transition between positive and negative contributions (Chen et al., 2026). To avoid overinterpreting unstable patterns, only thresholds that appeared within the main distribution range of the predictor and showed consistent response patterns were discussed. Therefore, to complement our understanding of CLI, this study selected a series of hydrometeorological factors – including VPD, SMsurf, SMroot, Temp, Pre, WS, and Rad – and employed the XGBoost algorithm to quantify their respective influences on R_max and T_opt (Piao et al., 2020). To increase the sample size, all time-series data from each grid cell were extracted. The dataset was divided into training and test subsets at a ratio of 80:20. Bayesian optimization combined with five-fold cross-validation was used for hyperparameter tuning. For the R_max model, the learning rate was set to 0.06, the number of trees to 1431, and the maximum depth to 4. For the T_opt model, the learning rate was set to 0.07, the number of trees to 1287, and the maximum depth to 4. The model performance was evaluated using the mean coefficient of determination (R²) and root mean square error (RMSE) obtained through five-fold cross-validation (Fig. S19).

This study further employed the improved Peter-Clark Momentary Conditional Independence plus algorithm (PCMCI+) (Runge et al., 2019; Runge, 2020) to investigate the causal mechanisms between R_max and T_opt, and hydrometeorological factors. These factors included VPD, SMsurf, SMroot, Temp, Pre, WS, and Rad. The time series of all variables were detrended and converted into anomalies to satisfy the stationarity requirements of time-series analysis. The causal diagnostic analysis was then implemented using the PCMCI+ algorithm in the Python package Tigramite (https://github.com/jakobrunge/tigramite, last access: 20 April 2025). PCMCI+ iteratively identifies directed causal associations among variables through conditional independence tests. To reduce uncertainty in causal inference, we constrained the network to physically plausible associations consistent with established mechanisms. Instead of the conventional ParCorr method, a Gaussian process-based nonlinear conditional independence test (GPDC) was applied, and p-values were computed through permutation testing (500 bootstrap iterations), with significance determined using an FDR-corrected threshold of p< 0.05. Compared with structural equation modeling, PCMCI+ does not require a prespecified model and can automatically infer causal structures from time-series data (Poppe Terán et al., 2023). A detailed methodological description of the PCMCI+ algorithm is provided by Runge (2020).

4 Results

4.1 Spatiotemporal Characteristics of R_max and T_opt

Based on remote sensing observations and the TRENDY multi-model ensemble, we computed the R_max and T_opt trends between ecosystem WUE and drought (Fig. 1). The results indicated that the R_max derived from both datasets declined significantly at a rate of −0.0003 yr⁻¹ (p< 0.01). Furthermore, the T_opt calculated from remote sensing observations and the TRENDY multi-model ensemble increased at rates of 0.0155 month yr⁻¹ and 0.0222 month yr⁻¹ (p< 0.01), respectively. Notably, the R_max derived from the TRENDY multi-model ensemble were approximately 0.02 lower than those from remote sensing observations, whereas the T_opt were about 0.5 months higher; however, both datasets exhibited nearly identical temporal trends.

Figure 1 Temporal trends of the R_max and T_opt between global WUE and SPEI calculated from (a) the remote sensing observations and (b) the TRENDY multi-model ensemble using an 18-year moving window.

The spatial patterns of the multi-year mean values and trends of R_max and T_opt were generally consistent, although the R_max derived from the TRENDY multi-model ensemble was smaller and the T_opt was larger than those obtained from remote sensing observations (Fig. 2). It was also noteworthy that the absolute values of the R_max and T_opt trends estimated from remote sensing observations were slightly larger than those derived from the TRENDY multi-model ensemble. Specifically, in the high-latitude regions of the Northern Hemisphere, the R_max generally showed a decreasing trend, with the exception of the western Siberia, and central North America. Notably, the decrease was most pronounced in the central Siberian. In contrast, the Southern Hemisphere exhibited more regions where the R_max have increased, particularly in areas such as the central-northern Africa, southern Africa, northern Australia, and the Amazon Basin, where the increasing rate was notably high. The T_opt exhibited an opposite trend to the R_max, showing a clear increase between 30 and 60° N, while the remaining regions primarily exhibited a decreasing trend.

Figure 2 Spatial distribution of the mean and change rates of R_max and T_opt calculated from the remote sensing observations and the TRENDY multi-model ensemble using an 18-year moving window. From left to right, the three columns represent remote sensing observations, the TRENDY multi-model ensemble and statistics across different vegetation types and aridity gradients. Points indicate statistically significant trends (p< 0.01).

From the perspectives of vegetation type and drought gradient, the mean values of R_max and T_opt showed relatively small differences, whereas their trends exhibited pronounced heterogeneity across these two dimensions (Fig. 2c). Specifically, the mean R_max varied little across aridity gradients but differed among vegetation types, with Shrub showing higher mean values than other vegetation types, particularly Shrub in HR, which reached 0.54. The mean characteristics of T_opt were consistent with those of R_max, with the highest mean T_opt observed in HR across aridity gradients and among shrubs across vegetation types. In terms of the rate of change, R_max showed an overall decreasing trend globally but increased in AR, with the highest growth rates observed in Forest and Shrub, at 0.0017 and 0.0021 yr⁻¹, respectively. Across vegetation types, R_max also tended to increase in Pasture. In contrast, the rate of change in T_opt showed marked regional differences across aridity gradients, with higher increasing rates in AR, reaching up to 0.1181 month yr⁻¹ in Pasture, while SH showed a decreasing trend at −0.0010 month yr⁻¹. Across vegetation types, T_opt increased more rapidly in Cropland and Pasture, whereas Forest displayed a decreasing trend, with the largest decline occurring in Forest within SAR at −0.0394 month yr⁻¹.

4.2 Separation of the Driving Factors of the changes in R_max and T_opt

Since the spatial distributions simulated by individual DGVM exhibited certain deviations (Figs. S9–S11), this study adopted TRENDY multi-model ensemble to separate the effects of CO₂ fertilization, CLI, and LCC on R_max and T_opt. The results from the TRENDY multi-model ensemble indicated that the negative effect of CO₂ on the R_max and the positive effect on T_opt were primarily concentrated in the high-latitude regions of the Northern Hemisphere. However, the opposite contributions were observed in the central-southern North America, central and eastern Siberia, as well as southern Central Asia and the Tibetan Plateau. In the Southern Hemisphere, CO₂ exerted a predominantly positive effect on the R_max and a negative effect on T_opt, revealing an overall opposite spatial pattern. In comparison to CO₂, the contributions of CLI and LCC were relatively smaller. The negative contribution of CLI to the R_max was mainly distributed in the southern Europe, western Siberia, and the Amazon Basin, while the negative contribution of LCC was primarily concentrated in the eastern Europe and central-northern Africa. Both factors showed an overall opposite spatial distribution pattern for their contributions to T_opt relative to the R_max (Fig. 3).

Figure 3 The contributions of CO₂, CLI, and LCC to the R_max and T_opt under an 18-year moving window. (The left and right columns represented the contributions of variables to the R_max and T_opt, respectively, with the contributions of CO₂, CLI, and LCC shown from top to bottom.)

At the global scale, the contribution of CO₂ to the R_max was approximately $-\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ yr⁻¹, while the contributions of LCC and CLI were approximately $-\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{4}}$ and $-\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{4}}$ yr⁻¹, respectively. The contributions of CO₂, LCC, and CLI to T_opt were 2.3$\times {\mathrm{10}}^{-\mathrm{2}}$, 0.2$\times {\mathrm{10}}^{-\mathrm{2}}$, and $-\mathrm{0.2}\times {\mathrm{10}}^{-\mathrm{2}}$ month yr⁻¹, respectively (Fig. 4). Further analysis across vegetation types and drought gradients revealed that CO₂ generally exerted a negative contribution to R_max, with substantial variation among vegetation types. The strongest negative contribution occurred in Cropland, particularly in Cropland within AR ($-\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{2}}$ yr⁻¹). Conversely, Forest and Shrub in AR and Pasture in HR exhibited relatively high positive contributions of 0.8$\times {\mathrm{10}}^{-\mathrm{3}}$, 1.2$\times {\mathrm{10}}^{-\mathrm{3}}$, and 1.3$\times {\mathrm{10}}^{-\mathrm{3}}$ yr⁻¹. In contrast, the contribution of LCC to R_max was negligible, with an absolute magnitude less than 0.1$\times {\mathrm{10}}^{-\mathrm{3}}$ yr⁻¹. The contribution of CLI to R_max was strongly dependent on drought gradients, with its negative effect gradually intensifying from AR to HR. For T_opt, the contributions of LCC and CO₂ exhibited broadly consistent patterns across vegetation types and drought gradients. Across drought gradients, the largest contributions of CO₂ and LCC occurred in AR, while across vegetation types, their effects were much greater in cropland and pasture, with the maximum observed in pastures within AR, where the contributions of CO₂ and LCC reached 0.0815 and 0.0196 month yr⁻¹, respectively. The contribution of CLI to T_opt was predominantly negative in AR and SAR but positive in more humid regions (SH and HR). Specifically, the strongest negative contribution of CLI to T_opt occurred in pastures within AR (−0.0196 month yr⁻¹), whereas the strongest positive contribution was observed in pastures within SH.

Figure 4 The contributions separated by DGVMs under an 18-year moving window, along with statistical plots for different vegetation types and drought gradients.

Based on the independent contributions of CO₂, LCC and CLI, their spatially dominant regions were further identified (Fig. 5). The results indicated that, for both R_max and T_opt, CO₂ dominated the largest area, followed by CLI, while LCC accounted for the smallest proportion. For R_max, the regions dominated by negative CO₂ contributions accounted for 27.3 % of the total area and were mainly distributed across the eastern Europe, western Siberia, Central Asia, southeastern South America, central and eastern Africa, and Australia. Regions dominated by positive CO₂ contributions covered a slightly smaller proportion (approximately 25.0 %), mainly concentrated in the central and eastern Siberia, south-central Africa, and eastern South America. Regions dominated by CLI accounted for 34.5 % of the total area and were mainly distributed along the margins of the CO₂-dominated zones, whereas those dominated by LCC accounted for only 13.2 %. For T_opt, regions dominated by positive CO₂ contributions accounted for 38.2 % of the total area and were concentrated between 30 and 60° N. In contrast, regions dominated by negative CO₂ contributions accounted for 20.2 % and were mainly distributed in high-elevation regions such as the eastern Siberia, the Tibetan Plateau, and central-southern North America. The regions dominated by CLI and LCC were relatively scattered.

Figure 5 Spatial distribution of the dominant factors of CO₂, LCC, and CLI on (a) the R_max and (b) T_opt under an 18-year moving window. The bar chart shows the proportions of the dominant areas of CO₂, LCC, and CLI.

4.3 The effects of hydrometeorological factors on the WUE-drought coupled relationship

Using XGBoost and SHAP, we quantified the contributions of hydrometeorological factors to the coupling relationship between WUE and drought (Fig. 6). From a global perspective, SMsurf was identified as the most important driving factor for both R_max and T_opt, yet it exhibited nearly opposite influence patterns. Specifically, the effect of SMsurf on R_max displayed a unimodal pattern: it reduced R_max when SMsurf was below 0.22 m³ m⁻³, promoted R_max between 0.22 and 0.44, and peaked around 0.35 m³ m⁻³. SMsurf showed a significant nonlinear negative relationship with T_opt: as SMsurf increased, its effect on T_opt shifted from positive to negative at approximately 0.25 m³ m⁻³. The effects of SMroot on R_max and T_opt were opposite to those of SMsurf, although the relative importance of SMroot was lower. Another key driving factor for both R_max and T_opt was Temp, ranking third and second in relative importance, respectively. Specifically, Temp promoted R_max at temperatures below 5.01 °C, inhibited it between 5.01 and 18.38 °C with a local minimum near 12 °C, and again exerted a positive effect above 18.38 °C. For T_opt, Temp exerted predominantly negative effects below 1.66 and above 22.80 °C, while showing positive effects between these thresholds. Other hydrometeorological factors also exhibited nonlinear patterns and threshold effects in their influences on R_max and T_opt.

Figure 6 Global and local SHAP analyses of the contributions of hydrothermal factors to (a) R_max and (b) T_opt. In the local SHAP analysis, the variables from top to bottom are Pre, Rad, SMroot, SMsurf, Temp, VPD and WS on R_max and T_opt.

The relative importance rankings and influence mechanisms of hydrothermal factors on R_max and T_opt varied slightly across vegetation types and drought gradients (Figs. S20–S27). In Cropland, VPD was the most important driver of R_max, exhibiting a single-threshold effect: VPD exerted a negative influence when exceeding 1.13 hPa and a positive influence when below this threshold. Rad was the dominant driver of T_opt, showing two distinct thresholds: Rad promoted T_opt below 145.69 W m⁻², inhibited it between 145.69 and 204.51 W m⁻², and promoted it again above 204.51 W m⁻², with a local SHAP maximum at 204.51 W m⁻². For Forest, Temp was the most important factor for both R_max and T_opt. Temp reduced R_max between 2.23 and 18.65 °C, with the strongest inhibition near 10 °C, whereas it increased T_opt between −0.99 and 22.00 °C, showing the strongest promotion around 10 °C. For Pasture, Pre and SMsurf were the most critical factors influencing R_max and T_opt, respectively. Pre decreased R_max when below 438.46 mm but increased it above this value, while SMsurf promoted T_opt below 0.24 m³ m⁻³ and inhibited it above. For Shrub, SMsurf was the most important driver for both R_max and T_opt with thresholds of 0.23 and 0.20 m³ m⁻³, respectively, showing opposite effects across the thresholds.

Across drought gradients, Pre and SMsurf were the most important drivers of R_max and T_opt, respectively, in AR. The effect of Pre on R_max shifted from negative to positive at 245.03 mm, while that of SMsurf on T_opt changed from positive to negative at 0.12 m³ m⁻³. In SAR, Temp and Rad were the most critical drivers of R_max and T_opt, respectively. The influence of Temp on R_max exhibited pronounced nonlinearity: Temp decreased R_max below 16.91 °C (with a local minimum near 13 °C) but promoted T_opt above 16.91 °C, showing a local maximum near 23 °C. Rad inhibited T_opt below 206.65 W m⁻² but promoted it above this threshold. In SH, VPD and Rad were the dominant drivers of R_max and T_opt, respectively; VPD increased R_max when exceeding 0.58 hPa, while lower VPD values reduced it, and the positive effects of Rad on T_opt were concentrated between 598.44 and 1472.70 W m⁻². In HR, Temp was the most critical driver for both R_max and T_opt, with their partial dependence plots exhibiting nearly opposite patterns.

The results derived from PCMCI+ revealed a complex causal network between hydrothermal factors, R_max, and T_opt (Fig. 7). At the global scale, WS exhibited a negative causal relationship with R_max, while R_max, in turn, negatively influenced Pre, Rad, and Temp. Although many factors did not exert a direct causal effect on R_max, they influenced it indirectly through various pathways. For instance, VPD affected R_max indirectly via WS, while Rad influenced Temp, which subsequently impacted WS and ultimately affected R_max. For T_opt, Rad and WS exhibited direct negative and positive causal relationships, respectively. Conversely, T_opt exerted a direct positive causal effect on Pre and a direct negative effect on Temp. Additionally, SMroot and SMsurf were causally linked to T_opt, though the direction of their influence remained unclear.

Figure 7 Causal network diagram of hydrothermal factors, R_max, and T_opt across different vegetation types, drought gradients, and the global scale.

To further explore the causal network under different vegetation types and drought gradients, we conducted an in-depth statistical analysis, examining interactions within overlapping regions of vegetation type and drought gradient (Figs. S28–S31). In terms of drought gradients, the causal network structures in AR and SAR were relatively simple. In AR, R_max was directly influenced by Rad, while T_opt showed little causal linkage with other hydrothermal factors. In SAR, R_max was not directly affected by hydrothermal factors but exhibited direct causal relationships with VPD and WS. Meanwhile, T_opt was positively influenced by Temp and Pre, negatively influenced by WS, and had an undefined negative causal relationship with SMsurf. Conversely, the causal relationships in SH and HR were considerably more intricate. Beyond the complex interactions among hydrothermal factors and their indirect effects on R_max and T_opt, hydrothermal factors also directly influenced R_max and T_opt through multiple pathways. For example, in SH, both Pre and Rad exerted direct negative causal effects on R_max, while Temp had a direct positive influence. From the perspective of vegetation types, the causal networks in Cropland, Pasture, and Forest were relatively complex, whereas Shrub exhibited a simpler structure. In Cropland, R_max was directly and negatively influenced by Rad, while T_opt was directly and positively affected by Rad but negatively impacted by SMroot and WS. In Pasture, R_max was negatively influenced by Pre, while T_opt was negatively affected by SMroot. In Forest, R_max was negatively influenced by Rad, whereas T_opt was negatively influenced by VPD and WS but positively affected by Pre and Temp. In Shrub, SMroot exerted a direct negative influence on R_max, while Temp had a direct positive effect. T_opt, on the other hand, was mainly influenced by Pre through a direct positive causal relationship and was negatively influenced by SMsurf and Temp.

5 Discussion

5.1 Robustness of the results

To ensure the robustness of the remotely sensed WUE observations, we averaged three sets of GPP and ET products and validated the results against data from 85 FLUXNET 2015 sites (Fig. S1, Table S1). The results indicated a high consistency between the multi-product average WUE from remote sensing and the site data (R= 0.62, Fig. S2). To confirm the alignment of the selected DGVMs with the remote sensing observations, we performed calculations for all models within TRENDY v12 and identified seven models whose trends under the S3 scenarios were consistent with the remote sensing observations. As shown in the Taylor diagram, the multi-model ensemble results exhibited high correlation with the remote sensing observations and lower standard deviation, supporting the credibility of the separated contributions of CO₂, LCC, and CLI (Fig. S6). To mitigate the impact of interannual variations on trend significance, all calculations were performed within an 18-year moving window. To test the robustness of the R_max and T_opt under different window choices, we repeated all calculations using 16 and 20-year windows, which yielded similar conclusions (Figs. S7–S17). Therefore, we considered our results to be highly reliable.

5.2 Spatiotemporal variation characteristics of the coupling relationship

This study pioneers the identification of optimal time scale for WUE response to drought coupling, demonstrating significant spatiotemporal heterogeneity. Temporally, the optimal WUE-drought response scale exhibited an upward trend (Fig. 1), indicating a delayed WUE response to drought and, consequently, an enhanced buffering capacity of vegetation against drought has been enhanced (Peters et al., 2018; Frank et al., 2015). This phenomenon likely stems from climate change-driven drought intensification, compelling vegetation to utilize slow-response water sources like deep soil moisture or groundwater (Miguez-Macho and Fan, 2012). However, Tang et al. (2024) employed GPP as a proxy for vegetation and found that vegetation sensitivity to drought increased and its response time shortened, which appears to contradict our findings. However, that study focused only on single indicators such as GPP, which reflect vegetation photosynthesis, whereas WUE likely plays a more complex role in regulating water and carbon dynamics. In addition, differences in the definition of sensitivity, data sources, methods, research periods, and research areas might lead to different conclusions. Different types of droughts might also result in different outcomes. Therefore, our findings did not necessarily refute previous views. However, we firmly believe that considering the time scale in the response of WUE to drought is a crucial step in revealing the dynamic relationship between WUE and SPEI.

Across drought gradients, the optimal drought scale in arid regions was significantly shorter than in humid regions. This phenomenon could be explained by the resource limitation hypothesis (Huxman et al., 2004) , which posited that in humid regions, plant growth was generally limited by light or nutrient availability, whereas in arid regions, it was primarily constrained by water scarcity (Maurer et al., 2020; Knapp et al., 2024). Larger T_opt were mainly concentrated in humid, energy-limited regions, where drought conditions caused by insufficient precipitation were often accompanied by high temperatures and intense solar radiation, thereby enhancing the resistance of WUE to drought (Gentine et al., 2019; Walther et al., 2019). For instance, Miller et al. (2023) reported a sustained increase in plant photosynthesis during spring droughts across most vegetated areas of the Northern Hemisphere. Similarly, during recent drought events such as those in Europe in 2018 and 2022, several energy-limited ecosystems also exhibited enhanced plant activity (Bastos et al., 2020). Differences in vegetation adaptation strategies along drought gradients might also have led to distinct response times. Vegetation in arid regions optimized water use through stomatal dynamics, whereas humid ecosystems buffered drought effects via deep rooting systems and canopy shading (Klein et al., 2011). A shorter response time might indirectly indicate higher sensitivity (Vicente-Serrano et al., 2013); thus, our findings also supported previous conclusions that vegetation in arid regions exhibited stronger drought sensitivity (De Keersmaecker et al., 2015; Seddon et al., 2016). However, we found that this pattern was shifting between arid and humid regions. Specifically, the optimal drought scale in arid regions showed the most rapid increase, whereas the rate of change in humid regions was much smaller or even declined. We attributed this discrepancy in the consideration of time scales.

From the perspective of vegetation types, the optimal drought scale increased rapidly in Cropland and Pasture, increased slightly in Shrub, but decreased in Forest. This pattern might result from the stronger influence of human activities on Cropland and Pasture, where irrigation and fertilization altered the water and nutrient conditions of these ecosystems (Jaramillo et al., 2018). On the one hand, irrigation maintained adequate soil moisture at all times. On the other hand, sufficient nitrogen supply ensured high photosynthetic nitrogen-use efficiency in crop leaves, indicating that even under partial stomatal closure and limited intercellular CO₂ concentration, leaves were still able to efficiently utilize the available CO₂ for photosynthesis. In contrast, Shrub and Forest were mostly subject to natural conditions. The slight increase in T_opt observed in Shrub might have reflected their conservative, isohydric water-use strategy (Yang et al., 2016). The decrease in T_opt in Forest might have been attributed to continued transpiration under water stress, which persisted due to their extensive canopies and residual hydraulic function, thereby shortening the response time to drought (Novick et al., 2016). These findings collectively reveal vegetation-drought coupling as a nonlinear outcome of climate forcing, ecological adaptation, and anthropogenic intervention, necessitating multi-scale models and sustained monitoring to unravel critical thresholds (Reichstein et al., 2013).

It is also important to emphasize that drought is a continuously evolving process, and the transition from one type of drought to another is referred to as drought propagation (Zhou et al., 2021). Drought generally begins as meteorological drought caused by reduced precipitation or increased potential evapotranspiration and may subsequently result in soil drought (Ma and Yuan, 2024). However, this propagation does not occur instantaneously because the soil water-storage capacity can buffer and delay the response of WUE to meteorological drought. Therefore, the decoupling of WUE from meteorological drought is closely linked to the propagation of meteorological drought into soil drought. Nevertheless, further research is required to investigate the coupling relationship between soil drought and WUE and to determine how soil drought regulates the coupling between meteorological drought and WUE.

5.3 Influencing factors of coupling relationship

TRENDY multi-model simulations identified CO₂ fertilization as the dominant driver weakening WUE-SPEI coupling, as evidenced by a decrease in the maximum correlation coefficient and an increase in the optimal drought timescale (Fig. 4). On the one hand, elevated atmospheric CO₂ enhances the carboxylation efficiency of Rubisco, thereby directly promoting GPP and enabling vegetation to maintain relatively high carbon assimilation without increased water supply; consequently, the synchrony between decreases in GPP and SPEI is weakened. On the other hand, elevated atmospheric CO₂ induces partial stomatal closure, thereby reducing transpiration water loss (Keenan et al., 2013). These mechanisms play a critical role in preventing rapid declines in WUE under drought conditions, thereby directly reducing its sensitivity to meteorological drought. Moreover, the CO₂ fertilization effect promotes deeper root distribution and adaptive changes in hydraulic traits, leading plants to increasingly rely on water stored from past precipitation in deeper unsaturated soil layers or groundwater to maintain carbon-water balance under drought stress (Zhang et al., 2026). This enhances the ecohydrological buffering capacity of vegetation against short-term meteorological drought, as reflected by an increase in the optimal drought timescale. Theoretically, these processes occur at the leaf scale, enhancing intrinsic water use efficiency (Li et al., 2025), however, cross-scale dependencies allow these effects to propagate to the ecosystem scale. Land-use changes exert secondary but critical impacts, including tropical deforestation, forest degradation, and afforestation programs, all of which gradually altered global vegetation patterns. In addition, agricultural management practices such as irrigation, fertilization, and crop modification did not change vegetation cover types but still influenced vegetation adaptability to environmental conditions. Chen et al. (2019) found that human land-use practices might have contributed to more than one-third of the global vegetation increase, particularly in extensive cropland regions – most notably in China (25 %) and India (6.8 %). Therefore, we reasonably inferred that LCC driven by human interventions enhanced ecosystem resilience to drought and other extreme events, particularly in croplands. This was manifested in our results as LCC contributed a negative effect in R_max and a positive effect in T_opt.

Although TRENDY multi-model simulations indicated that the influence of climate change was minimal, this study employed XGBoost and SHAP to reveal how meteorological factors affected the WUE–SPEI coupling relationship. SMsurf was identified as the most important driving factor, exhibiting nearly opposite effects on R_max and T_opt and a significant threshold response (Fig. 6). In fact, this remained closely associated with increased atmospheric CO₂. Previous studies suggested that although elevated atmospheric CO₂ could reduce plant water loss by decreasing stomatal conductance – thereby theoretically slowing soil moisture depletion (Gray et al., 2016) – the elevated CO₂ was accompanied by increased leaf temperature and accelerated evaporation from shallow soils, which offset this water-saving effect (Wilson et al., 1999), ultimately promoting soil moisture depletion (Kellner et al., 2019). As the direct “water reservoir” for vegetation, surface soil moisture limited plant water supply and reduced photosynthesis, thereby exerting a profound influence on the WUE–SPEI coupling relationship (Liu et al., 2020; Martínez-Vilalta et al., 2014). Along the drought gradient, precipitation was the most important factor influencing the WUE–SPEI coupling in arid regions, with SMsurf ranking third in importance. In contrast, in humid regions, temperature and radiation were the dominant factors, while SMsurf ranked fifth and sixth in importance for R_max and T_opt, respectively (Figs. S23, S26), consistent with previous findings. Vegetation in arid regions was typically water-limited, and precipitation, as its direct water source, played a crucial regulatory role in drought response (Li et al., 2022a). In contrast, in humid regions with sufficient water supply, vegetation responses to drought were primarily controlled by energy-related factors such as temperature and radiation (Liu et al., 2025). Across vegetation types, the dominant factors influencing the WUE–SPEI coupling were more clearly defined in Forest and Shrub, being temperature and surface soil moisture, respectively, which were closely related to their rooting patterns and drought-resistance strategies (Yang et al., 2025). In contrast, the dominant factors in Pasture and Cropland were more complex, likely resulting from vegetation structural diversity and intensive human interventions (Liu et al., 2018).

5.4 Limitations and future prospects

This study provides valuable insights into the dynamic coupling between WUE and meteorological drought; however, several limitations remain. Among various drought types, meteorological drought is the most prevalent, which is why SPEI was selected in this study (Zhou et al., 2024). Previous studies have demonstrated that other drought types are closely linked to meteorological drought, often triggered by precipitation anomalies and exhibiting certain lag effects; in contrast, heatwaves typically precede meteorological drought, indicating that heatwave-drought events evolve dynamically through a sequential triggering process (Yan et al., 2025a). Different drought types operate through distinct mechanisms and exert varying impacts on WUE. For example, soil drought is more directly associated with root-zone water availability for vegetation, whereas atmospheric drought affects WUE by regulating stomatal aperture. Furthermore, this study focuses solely on WUE at the ecosystem scale, while the coupling relationships between drought and WUE at the leaf and canopy scales remain unclear (Wang et al., 2026b). Therefore, future studies should further investigate the dynamic coupling between WUE across multiple scales and different types of heatwaves and droughts, as WUE at different scales responds differently across various stages of heatwave-drought propagation. The relatively short study period represents an additional limitation of this study. Although extending the analysis to more recent years would improve the timeliness and relevance of the results, such an extension was limited by the temporal coverage and consistency of available long-term GPP and ET datasets. Nevertheless, future studies should extend the analysis to more recent years as longer, more reliable, and temporally consistent datasets become available. Future research should also incorporate Coupled Model Intercomparison Project (CMIP) datasets to investigate projected changes in WUE-drought coupling and their driving mechanisms under different development scenarios. Such efforts will facilitate a more comprehensive understanding of vegetation responses to drought and provide a stronger theoretical foundation for evaluating carbon-water cycles under climate change.

6 Conclusions

Understanding the coupling relationship between WUE and drought is crucial for assessing the impact of drought on terrestrial ecosystem carbon-water cycles. To this end, we investigated the spatiotemporal patterns of the WUE-SPEI coupling relationship from 1982 to 2018 globally. Using the DGVMs, we separated the contributions of CO₂ fertilization, CLI, and LCC to the changes in WUE-drought coupling relationship. Furthermore, the XGBoost combined with SHAP and PCMCI+ were employed to elucidate the mechanisms and influencing pathways of hydrothermal factors on the coupling relationship. Our key findings included the following:

1. R_max decreased at a rate of −0.0003 yr⁻¹ (p< 0.01), while T_opt increased at a rate of 0.0155 month yr⁻¹ (p< 0.01), indicating the buffer capacity of WUE against drought was enhanced. The rates of change in R_max and T_opt varied substantially across drought gradients and vegetation types. R_max increased most rapidly in Shrub within AR at 0.0021 yr⁻¹ (p< 0.01), whereas the largest decreases occurred in Shrub within SAR and Forest within SH, both at −0.0020 yr⁻¹ (p< 0.01). T_opt increased most rapidly in Pasture within AR at 0.1181 month yr⁻¹ (p< 0.01), while it decreased most markedly in Forest within SAR at −0.0394 month yr⁻¹ (p< 0.01).

2. CO₂ fertilization was identified as the primary cause for the weakening of the coupling relationship. SMsurf was identified as the most critical driver, exhibiting nearly opposite effects and significant threshold effects on R_max and T_opt.

3. At the global scale, WS exhibited a direct negative causal relationship with R_max, while Rad and WS had direct negative and positive causal effects on T_opt, respectively. The complex causal networks under different vegetation types and drought gradients were also revealed.

Overall, our findings highlighted an enhanced resistance of WUE to drought, manifested as a reduced correlation and a delayed response time. Moreover, vegetation under different drought gradients gradually adjusted its water-use strategies to maintain WUE stability. These findings encourage a reassessment of the WUE–drought relationship, offering new insights to support the sustainable development of ecosystems under climate change.