# Deforestation Increases Vegetation Vulnerability to Drought Across Biomes

Chenwei Xiao, Sönke Zaehle, Stephen Sitch, Gregory Duveiller, Daniel E Pabon-moreno, Anthony P Walker, Jürgen Knauer, Fabienne Maignan, Christiane Schmullius, Ana Bastos

<figure><img src="image_1.png" /></figure>

► To cite this version:

Chenwei Xiao, Sönke Zaehle, Stephen Sitch, Gregory Duveiller, Daniel E Pabon-moreno, et al.. Deforestation Increases Vegetation Vulnerability to Drought Across Biomes. Global Biogeochemical Cycles, 2025, 39 (5), 10.1029/2024gb008378. hal-05081425

HAL Id: hal-05081425

https://hal.science/hal-05081425v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Global Biogeochemical Cycles®

## RESEARCH ARTICLE

10.1029/2024GB008378

**Key Points:**

* Satellite observations show that drought impacts on biomass are more severe when forests are converted to non-forests

* Dynamic global vegetation models generally agree on a more negative drought sensitivity for the forest-grassland transition

* Direct effects from land use and land cover changes increase negative drought impacts on deforestation hotspots, but are partially compensated by indirect effects

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

C. Xiao,
cxiao@bgc-jena.mpg.de

**Citation:**

Xiao, C., Zaehle, S., Sitch, S., Duveiller, G., Pabon-Moreno, D. E., Walker, A. P., et al. (2025). Deforestation increases vegetation vulnerability to drought across biomes. *Global Biogeochemical Cycles*, 39, e2024GB008378. https://doi.org/10.1029/2024GB008378

Received 30 SEP 2024
Accepted 8 MAY 2025

**Author Contributions:**

**Conceptualization:** Chenwei Xiao, Sönke Zaehle, Ana Bastos

**Data curation:** Stephen Sitch

**Formal analysis:** Chenwei Xiao

**Funding acquisition:** Sönke Zaehle, Ana Bastos

**Investigation:** Chenwei Xiao

**Methodology:** Chenwei Xiao, Sönke Zaehle, Gregory Duveiller, Daniel E. Pabon-Moreno, Ana Bastos

**Resources:** Sönke Zaehle, Anthony

P. Walker, Jürgen Knauer, Fabienne Maignan

**Software:** Daniel E. Pabon-Moreno

**Supervision:** Sönke Zaehle, Ana Bastos

**Visualization:** Chenwei Xiao

© 2025. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_4.png" /></figure>

## Deforestation Increases Vegetation Vulnerability to Drought Across Biomes

Chenwei Xiao<fnref n="1" />,²,<fnref n="3" /> Ⓓ, Sönke Zaehle<fnref n="1" /> Ⓓ, Stephen Sitch⁴ Ⓓ, Gregory Duveiller<fnref n="1" /> Ⓓ, Daniel E. Pabon-Moreno<fnref n="1" />, Anthony P. Walker⁵ Ⓓ, Jürgen Knauer⁶ Ⓓ, Fabienne Maignan⁷ Ⓓ, Christiane Schmullius⁸ Ⓓ, and Ana Bastos<fnref n="1" />,<fnref n="3" /> Ⓓ

**Abstract** Land use and land cover changes have altered terrestrial ecosystem carbon storage, but their impacts on ecosystem sensitivity to drought and temperature fluctuations have not been evaluated spatially over the globe. We estimate drought and temperature sensitivities of ecosystems using vegetation greenness from satellite observations and vegetation biomass from dynamic global vegetation model (DGVM) simulations. Using a space-for-time substitution with satellite data, we first illustrate the effects of vegetation cover changes on drought and temperature sensitivity and compare them with the effects estimated from DGVMs. We also compare simulations forced by scenarios with and without land cover changes to estimate the historical land cover change effects. Satellite data and vegetation models both show that converting forests to grasslands results in a more negative or decreased positive sensitivity of vegetation greenness or biomass to drought. Significant variability exists among models for other types of land cover transitions. We identify substantial effects of historical land cover changes on drought sensitivity from model simulations with a generally positive direction globally. Deforestation can lead to either an increased negative sensitivity, as drought-tolerant forests are replaced by grasslands based on model ensemble mean, or a decreased negative sensitivity, since forests under current land cover are predicted to exhibit greater drought resistance compared to those under pre-industrial land cover. Overall, our findings emphasize the critical role of forests in maintaining ecosystem stability and resistance to drought and temperature fluctuations, thereby implying their importance in stabilizing the carbon stock under increasingly extreme climate conditions.

**Plain Language Summary** Human activities, such as agricultural expansion and deforestation, have dramatically changed plant species and their distribution, affecting their growth. Given global warming and more severe and frequent drought events, it is crucial to understand whether land use and land cover changes make the terrestrial ecosystem more resilient or vulnerable to drought events and higher temperatures. We use three satellite observations as proxies for vegetation growth states and process-based models designed for simulating vegetation growth to evaluate whether different vegetation types respond differently to drought events and temperature fluctuations. We find that forests are more resistant to drought and temperature than grasslands and croplands as observed through satellites. Models also predict decreased ecosystem sensitivity to drought when forests become grasslands, but simulate sensitivity differently when vegetation changes to croplands. Land cover changes since the pre-industrial era have substantially altered ecosystem sensitivity to drought and temperature. Overall, our findings highlight the vital role of forests in maintaining ecosystem stability or making the ecosystem benefit from climate. Future land management strategies should consider these insights. We need to improve model performance in simulating croplands' response to drought and temperature changes to better guide future land use.

## 1. Introduction

Terrestrial ecosystems play a key role in the global carbon cycle and have been absorbing around 30% of anthropogenic carbon dioxide (CO₂) emissions (Friedlingstein et al., 2023). Expanding and intensifying land use and land cover changes (LULCC) have dramatically altered the land surface, directly perturbing terrestrial carbon

<footnote n="1">¹Max Planck Institute for Biogeochemistry, Jena, Germany, ²IMPRS Global Biogeochemical Cycles, Jena, Germany,</footnote>

<footnote n="3">³Department of Earth System Science and Remote Sensing, Leipzig University, Leipzig, Germany, ⁴College of Life and Environmental Sciences, University of Exeter, Exeter, UK, ⁵Environmental Sciences Division and Climate Change Science Institute, Oak Ridge National Laboratory, Oak Ridge, TN, USA, ⁶Hawkesbury Institute for the Environment, Western Sydney University, Penrith, NSW, Australia, ⁷Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France, ⁸Department for Earth Observation, Friedrich Schiller University Jena, Jena, Germany</footnote>
**Writing – original draft:** Chenwei Xiao
**Writing – review & editing:**

Sönke Zaehle, Stephen Sitch,
Gregory Duveiller, Daniel E. Pabon-
Moreno, Anthony P. Walker,
Jürgen Knauer, Fabienne Maignan,
Christiane Schmullius, Ana Bastos

storage and resulting in a small net land carbon sink (Pongratz et al., 2021). LULCC has contributed about 10%–15% to total anthropogenic CO₂ emissions (Friedlingstein et al., 2023) and has had substantial regional effects in recent decades, such as a decline in Eastern Europe's terrestrial carbon sink with a ~52% decrease of net carbon sink from LULCC and forestry over 2010–2019 (Winkler et al., 2023) and a carbon loss of 1.18 PgC from deforestation in the Brazilian Amazon over 2010–2019 (Qin et al., 2021). LULCC also indirectly impacts carbon storage by changing the vegetation response to extreme events, which is here referred to as ecosystem sensitivity, which has not been quantified. Under climate change with increasing extreme events such as heatwaves and droughts (IPCC, 2021), it is essential to understand how and to what extent LULCC modulates vegetation carbon sensitivity and stability to environmental stressors at the global scale.

LULCC affects ecosystem response to climate extremes by modulating ecosystem species compositions and forest structure. Various studies have proposed that ecosystem responses during drought vary with plant species (Anderegg et al., 2018; Condit et al., 1995; DeSoto et al., 2020; N. McDowell et al., 2008), biodiversity (Grossiord et al., 2014; Isbell et al., 2015; D. Liu et al., 2022), and vegetation traits such as tree height (Nepstad et al., 2007; L. Liu et al., 2021), tree density (Bottero et al., 2017; Zamora-Pereira et al., 2021) and tree age (F. Liu et al., 2021; N. G. McDowell & Allen, 2015; Phillips et al., 2009). These studies suggest that LULCC, which directly changes the above ecosystem properties, might strongly affect the ecosystem sensitivity to extreme events. For example, forest management practices, such as wood harvest, change tree age and height distributions (Picchio et al., 2020). Forest fragmentation and subsequent edge effects from deforestation increase wildfire susceptibility (Alencar et al., 2004; Cochrane & Laurance, 2002) and tree mortality (Laurance et al., 2000; Qie et al., 2017), change plant species composition and diversity (Haddad et al., 2015), and contribute to forest degradation (Bourgoin et al., 2024).

LULCC effects on ecosystem sensitivity can be partitioned into two classes: (a) direct effects caused by changes in vegetation composition through LULCC, which in turn directly alters ecosystem sensitivity because different vegetation types exhibit distinct responses to drought and temperature; and (b) indirect effects arising from changes in sensitivity of vegetation types between the simulation with and without LULCC. For example, after deforestation, due to the decreased tree population, the competition among trees for water and nutrients is relieved, and trees may suffer less from drought stress, which causes an increase in the drought resistance of forests. Dynamic global vegetation models (DGVMs) can be used to understand these two classes of LULCC effects. DGVMs have been commonly used to assess the contribution of LULCC to the global carbon budget or LULCC biophysical effects through multiple factorial simulations with different forcing (Devaraju et al., 2018; Friedlingstein et al., 2023). These factorial simulations also allow us to extract the effect of LULCC on ecosystem sensitivity from other confounding environmental factors as they are based on the same varying CO₂ and climate conditions but different LULCC conditions in DGVMs.

The availability of new remote sensing data has allowed for a global evaluation of the effect of LULCC on ecosystem sensitivity. These results can be used to compare with the LULCC effects predicted by DGVMs. For example, Xiao et al. (2023) quantified ecosystem resistance as the sensitivities of vegetation indices to drought duration and temperature and found substantial effects of land cover and land management on ecosystem resistance to droughts and temperature, including deforestation, forest harvest, and irrigation. However, it remains difficult to disentangle the land cover and land management effects from confounding environmental factors, such as background climate, which introduce uncertainties in the results. A common approach to isolate the local effect of land cover changes and minimize the potential effects of confounding environmental factors relies on a space-for-time substitution over a spatial moving window when comparing different land cover types, which has been applied to investigate the influence of LULCC on the surface energy budget and redistribution, surface temperature, and temperature diurnal cycle (Duveiller et al., 2018; Y. Li et al., 2015; Z. Liu et al., 2019), and cloud cover through land-atmosphere coupling (e.g., Duveiller et al., 2021). Under the assumption of a similar climate background in each local moving window, the variance of the target variables in the moving window is mostly explained by different fractions of land cover types. To our knowledge, this kind of approach has not yet been applied to disentangle the effects of land cover changes on ecosystem sensitivity from other confounding environmental factors.

In this study, we aim to disentangle the role of LULCC in shaping ecosystem sensitivity to droughts and maximum temperature anomalies at the global scale. We do so by combining a space-for-time substitution approach based on remote sensing products and factorial simulations comparison based on DGVMs. Within this
scope, we compare the spatial pattern of ecosystem sensitivity to droughts and maximum temperature anomalies derived from satellite observations and DGVMs. We then explore the effects of vegetation changes on ecosystem sensitivity and compare the results from satellite data and DGVMs. Finally, we estimate the historical LULCC effects and further attribute the effects to direct and indirect effects.

## 2. Materials and Methods

### 2.1. Satellite-Based Vegetation and Land-Cover Data

We used 21 years of NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data (MOD13C1 v061) for 2001–2020 at a 0.05° spatial resolution (Didan, 2021). NDVI is an empirical index exploiting the differences that green and healthy vegetation portrays in terms of reflectances in the red and near infra-red domains. However, this index typically saturates. Our analysis is thus based on kNDVI, a non-linear generalization of the NDVI that is calculated following the empirical equation proposed by Camps-Valls et al. (2021):

$$kNDVI = \tanh(NDVI^2) \quad (1)$$

While kNDVI can be calculated using more sophisticated formulations that better correlate with variables such as gross primary productivity (Camps-Valls et al., 2021, Q. Wang et al., 2023), the simpler formulation in Equation 1 is here preferred for convenience in terms of deriving it directly from the MODIS NDVI product and still mitigating the effect of saturation. kNDVI is also chosen here because it can be acquired at a relatively high spatial resolution compared to global products related to aboveground vegetation biomass with a period longer than 10 years, such as vegetation optical depth, which is directly proportional to the vegetation water content of the aboveground canopy biomass (X. Li et al., 2021). kNDVI enables us to use smaller local moving window sizes, which will better satisfy the assumption of a homogeneous climate background for space-for-time substitution. We calculated the annual mean kNDVI to investigate the inter-annual vegetation dynamics.

We used the MODIS MCD12Q1 v061 land cover data set based on the International Geosphere-Biosphere Program classification. The data set provides land cover fractions of 17 land-cover classes at 0.05° spatial resolution for 2001–2020 (Friedl & Sulla-Menashe, 2022). We then aggregated vegetation cover classes into four main vegetation cover classes (Forest, Shrub, Grass, and Crop) as described in Table S1 in Supporting Information S1. We also tested the robustness of our results to the choice of land cover maps by using the same method but with the plant functional type (PFT) maps derived from the European Space Agency (ESA)'s Climate Change Initiative (CCI) annual land cover maps (Harper et al., 2023) from 2001 to 2020 to account for uncertainties in land cover classification (Hartley et al., 2017; W. Li et al., 2018). The original spatial resolution is 300 m and we perform a first-order conservative remapping (area-weighted mean) to regrid them to 0.05° to fit the resolution of kNDVI. For our analysis, we exclude pixels with less than 5% vegetation cover or more than 50% water body cover over the selected period.

### 2.2. Climatic Variables

To estimate ecosystem sensitivity to drought and high temperatures, we use monthly soil moisture (SM) and monthly 2-m maximum temperature ($T_{max}$) variables from TerraClimate, which is available at ~4 km spatial resolution for the period 1958–2022 (Abatzoglou et al., 2018). TerraClimate employs climatically aided interpolation by merging high-spatial resolution climatological normals from the WorldClim data set with coarser resolution time-varying (i.e., monthly) data from Climatic Research Unit (CRU) Ts4.0, JRA-55, and root zone storage capacity data set (Abatzoglou et al., 2018). We perform a first-order conservative remapping to regrid the data to 0.05° spatial resolution and select the period of 2001–2020. TerraClimate SM is derived from monthly surface water balance data sets using a water balance model that incorporates reference evapotranspiration, precipitation, temperature, and interpolated plant extractable soil water capacity at a higher resolution.

### 2.3. Models and Simulations

We use simulations from DGVMs performed within the project “Trends and drivers of the regional-scale sources and sinks of carbon dioxide” (TRENDY-v11; Friedlingstein et al., 2022; Le Quéré et al., 2013; Sitch et al., 2024) that provide outputs of monthly total SM content, yearly mean carbon mass in vegetation (cVeg) at the pixel level as well as disaggregated per PFT (cVeg$_{PFT}$) between 2001 and 2020 at 0.5° or 1° spatial resolution. The models
fulfilling these requirements were CABLE-POP, CLASSIC, OCN, ORCHIDEE, and SDGVM (see Table S3 in Supporting Information S1 for a comparison of the model output and relevant processes included).

The models are forced by either the 6-hourly historical climate forcing over 1901–2021 from the CRU and Japanese Reanalysis forcing data or the monthly CRU data set on a 0.5° × 0.5° grid (Harris et al., 2014, 2020) and annual time-series of CO₂ over 1700–2021, derived from ice core CO₂ data merged with NOAA annual resolution from 1958 onwards (Lan et al., 2024). Land cover and land-use maps are based on the harmonized land-use change data LUH2-GCB2022 for the period 1700–2022 (an update to the LUH2 v2h data set; Hurtt et al., 2020), which includes land-use change data, including subgrid-scale transitions between primary forest, secondary forest, primary non-forest, secondary non-forest, cropland, pasture, rangeland, and urban land (Chini et al., 2021; Hurtt et al., 2020).

To extract the LULCC effects on ecosystem sensitivity to drought and temperature, we use two simulations: (a) simulation S2, which is forced by the time-invariant land use map corresponding to the year 1700, and reconstructed and observed historical environmental conditions of climate, N deposition, and atmospheric CO₂ concentrations; and (b) simulation S3, which is forced by observed historical land use and cover distribution in addition to the aforementioned time-variant environmental conditions. The simulated SM was calculated from the hydrological model in each DGVM.

Global vegetation is represented using different PFTs in different DGVMs. For a consistent comparison, we group the PFTs of each DGVM into four main vegetation classes (Forest, Shrub, Grass, and Crop) as described in Table S1 in Supporting Information S1 and calculate the cVeg for each main vegetation class (cVegclass) from cVegPFT.

## 2.4. Estimating Ecosystem Sensitivity

We use annual means of observed kNDVI and modeled cVeg at the pixel level as a measure for vegetation dynamics to calculate ecosystem sensitivity to droughts and Tmax. We also used simulated cVegclass to calculate ecosystem sensitivity to droughts and Tmax for different vegetation types for S2 and S3 simulations. We calculate the annual mean kNDVI to compare to results calculated from cVeg and cVegPFT reported on an annual basis in the TRENDY-v11 data set and the annual mean values reflect the vegetation conditions and disturbance from drought and heat stress during the whole year.

The annual mean kNDVI, cVeg, cVegclass, and Tmax, were detrended by subtracting the long-term trend of the linear fit. These time series anomalies are then z-score standardized by dividing by their standard deviation over the period 2001–2020. The detrended annual time series is used in the analysis because cVeg and cVegPFT are reported on an annual basis from the TRENDY v11 data set.

To identify drought periods, we first remove the long-term mean seasonal cycle (monthly mean SM over the period 1979–2020) and the linear trends of monthly SM data from TerraClimate and DGVM output. We define extremes as the occurrence of a value of SM below a threshold value (10th percentile) near the tails of the probability distribution of SM anomalies, that is, P(X ≤ SM(t)) < 0.1 (Seneviratne et al., 2012). For each pixel, we determine the distribution of monthly SM anomalies from 1979 to 2020 using kernel density estimation (KDE; Parzen, 1962; Harmeling et al., 2006). We use drought months from the S3 simulation to identify droughts in both the S3 and S2 simulations for each selected DGVM for a consistent comparison of drought events between the two simulations.

We then apply a multivariate linear regression model for each pixel and grid cell over 2001–2020 following:

$$Y_{anom}(t) = \alpha N(t) + \beta T_{anom}(t) + \phi Y_{anom}(t-1) + C + \epsilon(t) \quad (2)$$

where $Y_{anom}(t)$ and $Y_{anom}(t-1)$ denote the annual anomaly of kNDVI, cVeg, or cVegclass at a given year (t) and the previous year (t-1). In our study, year (t) is defined as a calendar year since cVeg is reported as the annual mean for each calendar year, as provided in the TRENDY model output. $N(t)$ denotes the number of drought months in year t, where the drought months are defined as the occurrence of a value of SM below the 10th percentile threshold value, given the probability distribution of SM anomalies, that is, P(X ≤ SM(t)) < 0.1. $T_{anom}(t)$ represents the deseasonalized and detrended $T_{max}$. $\alpha$ corresponds to the vegetation sensitivity to drought duration and $\beta$ indicates the vegetation sensitivity to $T_{max}$. The term $\phi$ corresponds to the dependence of
<table><caption>**Table 1**
Symbols of the Calculated Sensitivity to Drought and Temperature and Their Changes in This Study</caption><thead><tr><th>Category</th><th>Data source</th><th>Ecosystem sensitivity to drought</th><th>Ecosystem sensitivity to temperature</th></tr></thead><tbody><tr><td rowspan="3">Pixel level sensitivity</td><td>kNDVI from MODIS</td><td>α<sub>kNDVI</sub></td><td>β<sub>kNDVI</sub></td></tr><tr><td>DGVMs ensemble mean</td><td>α<sub>DGVMs</sub></td><td>β<sub>DGVMs</sub></td></tr><tr><td>Individual DGVM (OCN as an example)</td><td>α<sub>OCN</sub></td><td>β<sub>OCN</sub></td></tr><tr><td rowspan="3">Vegetation cover type sensitivity (F: forests; G: grasslands; C: croplands)</td><td>kNDVI from MODIS</td><td>X<sub>F</sub>, X<sub>G</sub>, X<sub>C</sub> (X could be any of the terms we defined for pixel-level sensitivity)</td><td></td></tr><tr><td>DGVMs ensemble mean</td><td></td><td></td></tr><tr><td>Individual DGVM (OCN as an example)</td><td></td><td></td></tr><tr><td rowspan="3">Sensitivity change (F → G: forests to grasslands; F → C: forests to croplands; G → C: grasslands to croplands)</td><td>kNDVI from MODIS</td><td>ΔX<sub>F→G</sub>, ΔX<sub>F→C</sub>, ΔX<sub>G→C</sub> (X could be any of the terms we defined for pixel level sensitivity)</td><td></td></tr><tr><td>DGVMs ensemble mean</td><td></td><td></td></tr><tr><td>Individual DGVM (OCN as an example)</td><td></td><td></td></tr></tbody></table>

vegetation on its state in the previous year $t-1$, and aims to control for memory effects of the previous year (Cranko Page et al., 2022; De Keersmaecker et al., 2016; Y. Liu et al., 2019) when estimating $\alpha$ and $\beta$. C is the intercept and represents the expected anomaly of vegetation greenness or biomass over 2001–2020 when all predictors are zero. $\epsilon$ is the residual term. We include the terms $\phi$ and C considering their physical relationship with vegetation growth and their contribution to the predictability of our Equation 2. More discussion can be found in Text S2 in Supporting Information S1.

A positive value of $\alpha$ implies an increase in vegetation greenness or biomass when there is a longer drought. Conversely, a negative value of $\alpha$ represents the detrimental impacts during drought periods, which we use as a measure of ecosystem resistance to droughts. A more negative $\alpha$ (i.e., with a higher absolute magnitude) corresponds to a lower ecosystem resistance, while an $\alpha$ close to zero implies high resistance. Similarly, we use negative $\beta$ as a measure of ecosystem resistance to temperature and positive $\beta$ implies beneficial responses of the ecosystem to higher temperatures. Based on Equation 2, drought sensitivity and temperature sensitivity derived from kNDVI with MODIS land cover maps are referred to as $\alpha_{kNDVI}$ and $\beta_{kNDVI}$. They are also calculated from cVeg at the pixel level and then averaged for the five selected DGVMs (i.e., the model ensemble mean), because all model results are based on cVeg, we omit the symbols of cVeg and subsequently refer to these sensitivities as $\alpha_{DGVMs}$ and $\beta_{DGVMs}$ (see Table 1).

## 2.5. Estimating Vegetation Transition Effects From Remote Sensing Data

We quantify the effects of pairwise land cover transitions on ecosystem sensitivity to droughts and temperatures ($\alpha$ and $\beta$). We do not investigate the pixels with abrupt land cover changes because the areas with strong and abrupt vegetation transitions from 2001 to 2020 can be limited and dividing the 20-year time series into multiple periods shortens the record and makes the estimate of ecosystem sensitivity to drought duration and temperature unreliable. Large changes in land cover compositions might directly affect kNDVI and result in overestimated ecosystem sensitivities. Instead, to estimate these effects from satellite data, we use the space-for-time substitution method established by Duveiller et al. (2018) and implemented in the Julia YAXArrays Toolbox package (Pabon-Moreno et al., 2023). This method allows us to infer the effects of land cover changes without focusing on the pixels with abrupt land cover changes or shifts in vegetation types in time. We use the ecosystem sensitivity estimated for pixels with stable land cover by excluding pixels with over 25% changes in any of the four vegetation types over 2001–2020 before applying the space-for-time substitution. We used a 5 by 5 moving window on 0.05° pixels. For every moving window, we use the linear regression model (Equation 2) to establish a relationship between the ecosystem sensitivity ($\alpha$ and $\beta$) estimated over 2001–2020 and mean land cover fractions X from MODIS land cover maps over 2001–2020, as described in Equation 3.
$$y = X\omega \Rightarrow \begin{cases} y_1 = \omega_1 \overline{x_{11}} + \omega_2 \overline{x_{12}} + \dots + \omega_m \overline{x_{1m}} \\ y_2 = \omega_1 \overline{x_{21}} + \omega_2 \overline{x_{22}} + \dots + \omega_m \overline{x_{2m}} \\ \vdots \\ y_n = \omega_1 \overline{x_{n1}} + \omega_2 \overline{x_{n2}} + \dots + \omega_m \overline{x_{nm}} \end{cases} \quad (3)$$

We use the four main vegetation types (Forest, Shrub and Savannas, Grass, and Crop) and five other non-vegetation land cover types (Permanent Wetlands, Urban and Built-up Lands, Permanent Snow and Ice, Barren, and Water Bodies) because the method requires all land cover fraction predictors to sum to 1. For each window of n (n = 25) pixels, the matrix X(n, m) comprises the fractions of each of the m (m = 9) land cover types as explanatory variables $\overline{x}_{nm}$ and the response variable y_n corresponds to ecosystem sensitivity to drought (α_n) or temperature (β_n), respectively. The vector of ω coefficients (ω_n) corresponds to the predicted sensitivity of each land cover type. However, X cannot be directly used in this analysis because the fractions of the 9 land cover types necessarily sum to one on land, which can lead to spurious correlations between different land cover fractions, and/or between land cover fractions and the response variable. Therefore, in this procedure, an additional transformation is applied to reduce the dimensionality of X using singular value decomposition. More technical details can be found in Duveiller et al. (2018).

After establishing linear regression models for all local moving windows over the globe, we estimate the effects of land cover transitions from vegetation cover type A to another vegetation cover type B with:

$$\Delta \hat{y}_{A \to B} = \hat{y}_{B} - \hat{y}_{A} \quad (4)$$

Here, y_A and y_B are calculated from the relationship established by Equation 3 by assuming that the pixel is covered by only land cover type A or B (x_A = 1 or x_B = 1). We refer to the calculated ecosystem sensitivity to drought for forests, grasslands, and croplands as α_kNDVI_F, α_kNDVI_G and α_kNDVI_C, and ecosystem sensitivity to temperature as β_kNDVI_F, β_kNDVI_G and β_kNDVI_C. The space-for-time substitution method does not assume the direction of the transition and treats the effect of forest-to-grassland and the effect of grassland-to-forest transitions as opposite. For conciseness, we only show the results for forest-to-grassland, forest-to-cropland and grassland-to-cropland transitions. We refer to the effects of the transitions between them as Δα_kNDVI_{F→C}, Δα_kNDVI_{F→G}, and Δα_kNDVI_{G→C} and Δβ_kNDVI_{F→C}, Δβ_kNDVI_{F→G}, and Δβ_kNDVI_{G→C} (see Table 1).

The method assumes that the average land cover compositions ($\overline{x}_{nm}$) is the major driver of the spatial variation of the ecosystem sensitivity to drought and temperature (y_n) over a local moving window. This assumption might be violated in some areas, leading to a low R² of the fitted linear model from Equation 3. To make sure that the vegetation transition effects are appropriately captured by Equations 3 and 4, we filter these effects using a threshold of R² > 0.2 of the fitted linear model from Equation 3. To calculate the pairwise change in sensitivities from land cover transitions, the method also requires enough co-occurrences of vegetation classes within the local moving window and tends to provide more reliable results if there are large and evenly balanced presences of both vegetation classes of interest. We use an index quantifying mutual presence (Ic) for each pair of vegetation classes and mask out those pixels whose Ic < 0.5 (local windows do not provide enough co-occurrences of two classes of vegetation of interest for the transition). We further exclude local moving windows with significant topographic variability, as topographical relief often corresponds to climatic gradients, undermining the space-for-time approach. Pixels are filtered based on three criteria: (a) the standard deviation of elevation within the local window must not exceed 50 m; (b) the difference between the mean elevation of the central pixel and the mean elevation across the local window should be under 50 m; and (c) the difference in elevation standard deviation between the central pixel and the entire local window must be less than 50 m. Finally, to address potential uncertainties of our space-for-time substitution approach, we only use the values of the effects of land cover transitions on ecosystem sensitivities (α and β) estimated from Equations 3 and 4 that fall inside the 20th and 80th percentiles for further analysis. We do not use a standard deviation-based outlier filtering because the standard deviation is sensitive to extreme values and is greater than 1 in these cases. Applying an SD-based outlier filtering criterion still retains physically unrealistic values. By filtering with the 20th–80th percentile, we ensure a more
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Illustration for partitioning historical land cover effects on ecosystem sensitivity into direct (i) and indirect effects (ii). Land cover types are represented using patches within a grid cell. Land cover changes occur in areas indicated by the red box, where the forest is converted to cropland. (a) A grid cell in the S2 simulation with fractions of forest and cropland identical to the pre-industrial land cover. (b) A grid cell where fractions of forest and cropland are identical to the pre-industrial land cover with ecosystem sensitivity to drought and temperature derived from the S3 simulation. (c) A grid cell in the S3 simulation with forest and cropland fractions averaged for 2001–2020. The thin red arrow indicates direct effects (i) caused by different ecosystem sensitivity between forest and cropland. The thin black arrow represents indirect effects (ii) caused by varying ecosystem sensitivity in the forest between S3 and S2, and the dashed black arrow represents indirect effects (ii) caused by varying ecosystem sensitivity in the cropland between S3 and S2.</figcaption></figure>

reliable representation of the vegetation transition effects on ecosystem sensitivity. After the above filtering processes, we keep 3% of the pixels with valid values.

## 2.6. Estimating Vegetation Transition Effects From the Model Output

We use cVeg from S3 simulation as Y, Tmax from model forcing data, and the number of drought months N in each year calculated from SM from each model output to calculate the global pattern of ecosystem sensitivity to drought and Tmax following Equation 2. We use S3 simulation because it is forced by historical LULCC, so that it is more comparable to the kNDVI than S2 simulation without LULCC. DGVMs also simulate cVeg per PFT and allow us to directly calculate the ecosystem sensitivity to droughts and Tmax for each main vegetation class by using cVeg\_class as Y, instead of using the space-for-time substitution with satellite data. We then calculate the effect of vegetation cover changes on the ecosystem sensitivity to drought in each pixel using:

$$ \Delta\alpha_{i \to j} = \alpha_{S3_j} - \alpha_{S3_i} \quad (5) $$

where $\alpha$ represents drought sensitivity, with $i$ and $j$ denoting two distinct vegetation cover classes, where vegetation cover class $i$ is converted to $j$. Similarly, we estimate the effect on ecosystem sensitivity to Tmax $\beta$. We refer to the model ensemble mean ecosystem sensitivity to drought for forests, grasslands, and croplands as $\alpha_{DGVMs_F}$, $\alpha_{DGVMs_G}$ and $\alpha_{DGVMs_C}$, and ecosystem sensitivity to temperature as $\beta_{DGVMs_F}$, $\beta_{DGVMs_G}$ and $\beta_{DGVMs_C}$ and the model ensemble mean effect of vegetation cover changes as $\Delta\alpha_{DGVMs_{F \to C}}$, $\Delta\alpha_{DGVMs_{F \to G}}$, and $\Delta\alpha_{DGVMs_{G \to C}}$ and $\Delta\beta_{DGVMs_{F \to C}}$, $\Delta\beta_{DGVMs_{F \to G}}$, and $\Delta\beta_{DGVMs_{G \to C}}$ (see Table 1).

We also estimate historical LULCC effects on ecosystem sensitivity to drought and temperature, that is, those that occurred, as opposed to the potential effects that can be expected from space-for-time substitution or from looking only at differences within S3 simulation. We further partition them into direct and indirect contributions by contrasting S3 and S2 simulations. Figure 1 illustrates the division of these two effects. Under S3 and S2 simulations, DGVMs are forced with the same transient climate and CO₂ concentrations; the only difference is that S3 has land cover changes but S2 has a fixed pre-industrial land use and land cover map, so we can compare them to extract the historical LULCC effects. We calculate the reconstructed ecosystem sensitivity $\hat{\alpha}$ and $\hat{\beta}$ for each pixel based on $\alpha_j$, $\beta_j$ from Equation 2, and the corresponding land cover fraction $f_j$ of each main vegetation class j
from S3 and S2 simulations following Equations 6 and 7. We then estimate the historical LULCC effects by comparing the $\hat{\alpha}$ and $\hat{\beta}$ between S3 and S2 simulations and partition the total effects into (a) direct effects, where the changes arise from changes in f; and (b) indirect effects, where the changes arise from changes in α following Equation 8:

$$ \hat{\alpha}_{S3} = \sum_{j=1}^{n} \alpha_{S3_j, \overline{fS3_j}} \qquad (6) $$

$$ \hat{\alpha}_{S2} = \sum_{j=1}^{n} \alpha_{S2_j, \overline{fS2_j}} \qquad (7) $$

$$ (\hat{\alpha}_{S3} - \hat{\alpha}_{S2}) = \underbrace{\left(\sum_{j=1}^{n} \alpha_{S3_j, \overline{fS3_j}} - \sum_{j=1}^{n} \alpha_{S3_j, \overline{fS2_j}}\right)}_{(i)} + \underbrace{\left(\sum_{j=1}^{n} \alpha_{S3_j, \overline{fS2_j}} - \sum_{j=1}^{n} \alpha_{S2_j, \overline{fS2_j}}\right)}_{(ii)} \qquad (8) $$

In these equations, $\hat{\alpha}$ denotes the reconstructed drought sensitivity. S2 and S3 refer to the two simulations and f denotes the fraction of each main vegetation class relative to the total vegetation cover fraction and j represents each of the n main vegetation classes (n = 3). Similarly, we also evaluate the historical LULCC effects for $\hat{\beta}$.

In Equation 8, we group the terms with the same α but different f into term (i), and terms with varying α but the same f into term (ii). Term (i) refers to the direct effects attributable to the difference in land cover fractions between S3 and S2. The difference arises from different α in different vegetation types. Term (ii) represents the indirect effects arising from differences in α per vegetation type j between S3 and S2. Different α estimated in S2 and S3 for the same vegetation types might be attributed to the differing water resource needs of the tree PFT following land cover changes. In S3, the reduction in tree cover can result in a decreased water requirement for the tree PFT to endure drought conditions, leading to a less negative sensitivity to drought at the grid cell level. Similarly, we also evaluated the impact of vegetation cover changes on temperature sensitivity.

We estimate α, β, Δα, Δβ, $\hat{\alpha}$ and $\hat{\beta}$ first for each model and then calculate their model ensemble mean. For pixel-wise comparison between DGVMs and kNDVI and among models, all the above metrics are regridded to 1° spatial resolution using first-order conservative remapping. We employ the two-sided Wilcoxon signed-rank test (Wilcoxon, 1945) at a 5% confidence level to evaluate the statistical significance of differences in the median of α, β and the effects of vegetation cover changes on α and β between DGVMs and kNDVI. We use the one-sided Wilcoxon signed-rank test (Wilcoxon, 1945) at a 5% confidence level to evaluate whether medians of α, β in forests are significantly higher than those in grasslands or croplands. We only show results aggregated for three main vegetation classes (Forest, Grass, and Crop) because only CABLE-POP has the Shrub PFT.

# 3. Results

## 3.1. Observation-Based Sensitivity and Comparison to Models

Globally, ecosystem sensitivity to drought estimated from kNDVI predominantly exhibits negative values, with 65% of the vegetated area showing negative $\alpha_{kNDVI}$. Exceptions were observed in some tropical forests and boreal regions (Figures 2a and S7b in Supporting Information S1). Compared to $\alpha_{kNDVI}$, ecosystem sensitivity to temperature $\beta_{kNDVI}$ shows a smoother spatial pattern that strongly depends on latitudes (Figures 2d and 2f). Vegetation north of 45°N and vegetation in tropical forests, such as Amazon rainforests, show positive sensitivity to yearly mean monthly maximum temperature ($\beta_{kNDVI}$), which indicates that vegetation kNDVI increases under higher temperatures. Negative $\beta_{kNDVI}$ is detected in other tropical regions such as Central America, East Brazil, and South Asia (Figure S7c in Supporting Information S1). The spatial pattern of $\beta_{kNDVI}$ in the tropics closely aligns with the distribution of dominant vegetation cover (Figure S4 in Supporting Information S1). Positive $\beta_{kNDVI}$ values are predominantly found in tropical forests, while negative $\beta_{kNDVI}$ values are mainly associated with shrublands, savannas, grasslands, and croplands. This pattern reflects the varying temperature ranges optimal for vegetation growth across different vegetation types. Areas with positive $\alpha_{kNDVI}$ and $\beta_{kNDVI}$ coincide with energy-limited areas (Denissen et al., 2022; Jiao et al., 2021). Whereas negative $\alpha_{kNDVI}$ and $\beta_{kNDVI}$, indicating negative impacts on vegetation greenness during droughts or hotter periods, are mostly in water-limited
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Estimated drought sensitivity α and temperature sensitivity β from kNDVI (a, d) and the modeled cVeg (b, e) over 2001–2020. The sensitivities are then averaged across Dynamic global vegetation models. Areas with a coefficient of variation (CV) > 3 between models are marked as gray stipplings. The averages for different latitudes and their standard deviations are shown on the right (c, f).</figcaption></figure>

regions. This is consistent with the previous findings that ecosystem response to water availability depends on whether the area is energy-limited or water-limited (Ciais et al., 2005; Flach et al., 2018).

To compare models to kNDVI, the global results are used, including the stippled areas with low agreement among models in Figures 2b and 2e. cVeg sensitivity to drought from DGVMs ensemble mean shows a similar spatial pattern to the estimation of drought sensitivity using kNDVI ($α_{kNDVI}$): DGVMs estimate negative values of $α_{DGVMs}$ in 83% of the vegetated areas and positive values in limited areas in tropical forests and boreal regions (Figure 2b). Drought sensitivity from DGVMs shows a smoother spatial pattern than drought sensitivity from kNDVIs, partly due to kNDVI's higher spatial resolution. Even after regridding drought sensitivity from kNDVI to a 1° spatial resolution, $α_{kNDVI}$ still shows higher spatial variability than $α_{DGVMs}$ (Figure S6 in Supporting Information S1). This may be attributed to kNDVI capturing fine scale processes not necessarily represented in models, as well as noise due to, for example, cloud cover and shadows. Additionally, the regression model may attribute non-drought effects, which are not accounted for in the predictors, to the coefficients including $α_{kNDVI}$. The spatial correlation between $α_{DGVMs}$ and regridded $α_{kNDVI}$ is quite low ($r = 0.19$), with $α_{OCN}$ showing the closest agreement ($r = 0.20$) of the five DGVMs (Figure S5 in Supporting Information S1). Models and kNDVI show a better agreement for the temperature sensitivity than drought sensitivity, with a spatial correlation of 0.67, and CLASSIC shows the closest agreement ($r = 0.64$) (Figure S5 in Supporting Information S1). Similar to $β_{kNDVI}$, $β_{DGVMs}$ tends to be positive in areas north of 45°N, but the simulated effects are weaker than $β_{kNDVI}$ (Figures 2e and 2f). In the tropics, $β_{kNDVI}$ shows divergent signals among areas dominated by different vegetation types. In the tropical forests, $β_{kNDVI}$ shows positive values, but $β_{DGVMs}$ is negative or near zero, showing no beneficial effects to higher temperatures. This discrepancy may suggest that models have difficulty simulating the positive responses of tropical forests to climate anomalies, including increased radiation and temperature, during periods of higher temperatures.

Our linear autoregressive model Equation 2 fits the modeled cVeg better than kNDVI observations over the globe (Figure S3 in Supporting Information S1). $R^2$ of Equation 2 fitted for kNDVI is higher than 0.5 in 12% of the areas, with a global mean of 0.27. Our model performs better in some areas in eastern South America, southern Africa, eastern Australia, and some boreal regions, where the $R^2$ is higher than 0.5. The model mean $R^2$ of Equation 2 fitted for modeled cVeg is higher than 0.5 in 51% of the areas with a global mean of 0.50.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Estimated drought sensitivity distributions (a) and temperature sensitivity distributions (b) from S3 cVeg\_class for different vegetation types from different models, the model ensemble mean, and kNDVI, with medians as black lines, interquartile ranges as colored shading, and the 10th and 90th percentile range as whiskers. The depicted distributions are limited to pixels where forest, grass, and crop coexist. Colored stars indicate vegetation types with sensitivity significantly lower than that calculated in forests, as determined by the one-sided Wilcoxon signed-rank test at a 5% confidence level.</figcaption></figure>

## 3.2. Contrasting Sensitivity for Different Vegetation Cover

We compare ecosystem sensitivity for different vegetation types, forests, grasslands, and croplands with results from kNDVI. In this study, we define higher sensitivity as values that are either less negative or more positive, rather than just larger in magnitude. Ecosystem kNDVI sensitivity to drought and temperature is less negative or more positive in forests than in grasslands and croplands. The predicted $\alpha_{kNDVI_f}$ shows a global median of -0.06 ($p$-value <0.05) with the lower and upper quartile range [Q1, Q3] = [-0.22, 0.13] per drought month, significantly higher than the median of $\alpha_{kNDVI_G}$ (-0.18 per drought month, [Q1, Q3] = [-0.36, 0.05], $p$-value <0.05) and $\alpha_{kNDVI_C}$ (-0.15 per drought month, [Q1, Q3] = [-0.35, 0.02], $p$-value <0.05) over the areas where forests, croplands, and grasslands coexist (one-sided Wilcoxon signed-rank test; $p$-value <0.05). Approximately 40% of the regions exhibit positive $\alpha_{kNDVI_f}$, particularly along the edges of the Amazon, Eastern Europe, and Russia, where $\alpha_{kNDVI}$ is also predominantly positive. In contrast, positive values of $\alpha_{kNDVI_G}$ and $\alpha_{kNDVI_C}$ are less common, covering only 28% and 26% of the areas, mainly at the edges of tropical forests in the Amazon and Central Africa, as well as in Eastern Europe and Russia. $\beta_{kNDVI_f}$ shows a median of 0.30 over the globe ([Q1, Q3] = [-0.03, 0.52], $p$-value <0.05), significantly higher than the median of $\beta_{kNDVI_G}$ (0.04, [Q1, Q3] = [-0.26, 0.28], $p$-value <0.05) and $\beta_{kNDVI_C}$ (0.06, [Q1, Q3] = [-0.25, 0.26], $p$-value <0.05) (one-sided Wilcoxon signed-rank test; $p$-value <0.05). We also find less negative sensitivity (i.e., higher resistance) to drought in forests than grasslands from most DGVMs and their ensemble mean. When we compare forests to croplands, DGVMs have divergent predictions of whether forests exhibit reduced sensitivity to drought than croplands. The median of $\alpha_{DGVMs_f}$ over the globe (-0.09 per drought month, [Q1, Q3] = [-0.14, -0.02] per drought month, $p$-value <0.05) is significantly higher than the median of $\alpha_{DGVMs_G}$ (-0.13 per drought month, [Q1, Q3] = [-0.21, -0.07] per drought month, $p$-value <0.05), indicating that forests tend to be more drought-resistant. The distribution of $\alpha_{DGVMs_C}$ (median of -0.07 per drought month, [Q1, Q3] = [-0.17, -0.02] per drought month, $p$-value <0.05) is similar to the distribution of $\alpha_{DGVMs_F}$ (Figure 3a). The median of $\beta_{DGVMs_f}$ is significantly higher than the median of $\beta_{DGVMs_G}$ (one-sided Wilcoxon signed-rank test; $p$-value <0.05) but their absolute difference is small. The difference between $\beta_{DGVMs_f}$ and $\beta_{DGVMs_C}$ is also small.

Similar to the spatial pattern of $\alpha_{DGVMs}$, the model ensemble mean estimates $\alpha_{DGVMs_f}$, $\alpha_{DGVMs_G}$, and $\alpha_{DGVMs_C}$ are predominantly negative with the exceptions of limited positive values in boreal and tropical regions (Figure S8b in Supporting Information S1). $\alpha_{DGVMs_f}$ is less negative than $\alpha_{DGVMs_G}$ and $\alpha_{DGVMs_C}$ in most areas. $\beta_{DGVMs_f}$, $\beta_{DGVMs_G}$ and $\beta_{DGVMs_C}$ show negative values mostly in tropical regions and positive values in areas at around 40°N to 60°N latitudes, Eastern US, and East Asia (Figure S8d in Supporting Information S1). $\beta_{DGVMs_G}$ is more negative compared to $\beta_{DGVMs_f}$ in the tropics but more positive in boreal regions, as reflected by a wider interquartile and 10th–90th percentile range in Figure 3b. Over the globe, there is only a small difference in the spatial medians of temperature sensitivity between forests and grasslands.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Land cover transition effects on drought sensitivity (a) and temperature sensitivity (b) retrieved from local moving windows based on kNDVI. The distributions of $\Delta\alpha$ and $\Delta\beta$ are shown inside each map. The median is shown as the blue number, indicated by the blue line. If the median significantly differs from zero (Wilcoxon signed-rank test; *p*-value <0.05), the corresponding number is displayed in bold. The zero line is shown as the black line.</figcaption></figure>

Among the five DGVMs, significantly higher medians of $\alpha_{DGVMS_F}$ than $\alpha_{DGVMS_G}$ are simulated in CABLE-POP, CLASSIC, and OCN (one-sided Wilcoxon signed-rank test; *p*-value <0.05). Only OCN simulates significantly higher $\alpha_{DGVMS_F}$ than $\alpha_{DGVMS_C}$. CABLE-POP and ORCHIDEE simulate significantly higher medians of $\beta_{DGVMS_F}$ than $\beta_{DGVMS_G}$. Only ORCHIDEE exhibit significantly higher medians of $\beta_{DGVMS_F}$ than $\beta_{DGVMS_C}$.

### 3.3. Vegetation Transition Effects on Ecosystem Sensitivity

By applying the local moving window method on $\alpha_{kNDVI}$ and $\beta_{kNDVI}$, we further estimate the potential vegetation transition effects on ecosystem sensitivity to droughts and temperatures (Figure 4). In this study, we define increased sensitivity as values that are either less negative or more positive, rather than just an increase in magnitude. Similarly, decreased sensitivity is defined as values that are either more negative or less positive. Our analysis reveals that when forests are converted to grasslands or croplands (also referred to as non-forests), the median of $\Delta\alpha_{kNDVI_{F-G}}$ is -0.08 per drought month and its lower and upper quartile ([Q1, Q3]) is [-0.24, 0.07] per drought month. The median of $\Delta\alpha_{kNDVI_{F-C}}$ is around -0.08 per drought month ([Q1, Q3] = [-0.19, 0.02] per drought month). $\alpha_{kNDVI}$ decreases when forests are converted to non-forests, predominantly in North America, Europe, East Asia, and the tropics (Figure 4a). The transition from grasslands to croplands results in small changes in $\alpha_{kNDVI}$, with a global median of 0.02 per drought month ([Q1, Q3] = [-0.08, 0.13] per drought month). Specifically for areas with negative $\alpha_{kNDVI}$, ecosystem sensitivity to drought increases in 63% of the areas when forests are converted to grasslands and in 70% of the areas when forests are converted to croplands. For temperature sensitivity changes, the median $\Delta\beta_{kNDVI_{F-G}}$ is -0.15 ([Q1, Q3] = [-0.34, 0.03]), and $\Delta\beta_{kNDVI_{F-C}}$ exhibits a median change of -0.19 ([Q1, Q3] = [-0.33, -0.06]). The transitions from forests to non-forests lead to a decrease in $\beta_{kNDVI}$ in over 71% of the areas, predominantly in North America, Europe, East Asia, and the tropics, as illustrated in Figure 4b. There is no significant change in $\beta_{kNDVI}$ for the transition between grasslands and croplands. The contrasts of $\alpha_{kNDVI}$ and $\beta_{kNDVI}$ between forests and non-forests underscore that deforestation or afforestation activities, in addition to their direct impacts on vegetation, also play a significant role in changing the ecosystem sensitivity.

Based on the $\alpha$ and $\beta$ calculated for different vegetation types from the five DGVMs, we analyze the effects of the transition between them on ecosystem sensitivity to droughts and temperatures. We calculated the model ensemble mean of these transition effects and compared them to the effects inferred from kNDVI (Figure 5). When the forests are converted to grasslands, $\alpha_{DGVMs}$ decreases globally with a median of -0.04 per drought month ([Q1, Q3] = [-0.10, 0.02] per drought month), consistent with the kNDVI estimate (Figure 5c). More than 66% of the areas exhibit a decrease in $\alpha_{DGVMs}$ when forests are converted to grasslands. This pattern is consistent for all models, except ORCHIDEE, which estimates a small positive median $\Delta\alpha_{DGVMS_{F-G}}$. For other transitions, the agreement between models and kNDVI is weaker. The transition from forests to croplands shows a median effect of 0.01 per drought month ([Q1, Q3] = [-0.06, 0.08] per drought month) on $\alpha_{DGVMs}$, which decreases in only 45% of the areas (Figure 5a), whereas the median of effects for the same transition estimated from kNDVI is around -0.08 per drought month and 70% of the areas show a decrease. OCN also shows a similar negative
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Model-based transition effects and their comparison with kNDVI. Global pattern of land cover transition effects on drought sensitivity (a) and temperature sensitivity (b) estimated from model mean. The distributions of $\Delta\alpha$ and $\Delta\beta$ are shown inside each map. The median is denoted in blue, indicated by the blue line. If the median significantly differs from zero (two-sided Wilcoxon signed-rank test; *p*-value <0.05), the corresponding number is displayed in bold. The zero line is shown as the black line. The median of the corresponding transition effects estimated from kNDVI is shown as the orange line. The distributions of transition effects estimated from kNDVI and different models are also shown in panels (c, d) with medians as black lines, interquartile ranges as colored shading, and ranges between the 10th and 90th percentiles as whiskers. The horizontal blue lines in panels (c, d) indicate the median of transition effects calculated from kNDVI. Colored stars mark models whose median is not significantly different from kNDVI (two-sided Wilcoxon signed-rank test; *p*-value <0.05).</figcaption></figure>

median of the forest-crop transition effects on $\alpha$ (-0.08 per drought month, [Q1, Q3] = [-0.22, 0.05] per drought month) to kNDVI but other models do not capture the negative impacts. The transition between grasslands and croplands has negligible effects on $\alpha_{kNDVI}$, but all models except ORCHIDEE predict increased $\alpha$.

Transition effects on $\beta$ inferred from the model ensemble mean and kNDVI are less consistent than the transition effects on $\alpha$ (Figures 4b and 5b). kNDVI indicates a decrease in $\beta_{kNDVI}$ but the model ensemble mean shows a close to 0 change in $\beta_{DGVMs}$ ([Q1, Q3] = [-0.11, 0.12]) when forests are converted to grasslands and a positive median change of 0.04 ([Q1, Q3] = [-0.09, 0.18]) when forests are converted to croplands (Figure 5d). Less than 50% of the areas exhibit a decreased $\beta_{DGVMs}$ (Figure 5b). Among the five DGVMs, CABLE-POP and ORCHIDEE show negative medians of $\Delta\beta_{F \to G}$, with values close to kNDVI. When grasslands are converted to croplands, the model ensemble mean simulates a positive median of $\Delta\beta_{DGVMs_{G \to C}}$ of 0.03 ([Q1, Q3] = [-0.07, 0.12]), but the median of $\Delta\beta_{kNDVI_{G \to C}}$ is close to 0 ([Q1, Q3] = [-0.11, 0.12]). ORCHIDEE and SDGVM also simulate a transition effect with a median close to 0 ([Q1, Q3] = [-0.12, 0.13] in ORCHIDEE and [Q1, Q3] = [-0.17, 0.15] in SDGVM), but CLASSIC simulates a strong positive effect on $\beta_{CLASSIC}$ for grass-to-crop transition with a median of 0.18 ([Q1, Q3] = [-0.02, 0.41]).

To conclude, ecosystem sensitivity to drought and temperature is estimated to be less positive or more negative (i.e., weaker resistance) when forests are converted to grasslands or croplands, while grass-to-crop transition does not result in a systematic change in sensitivity. To test the robustness of our results to land cover maps, we also used the ESA CCI PFT maps to calculate the kNDVI-based transition effects. The main conclusion remains the same (Figures 5c and 5d). Therefore, we consistently use the results from the MODIS MCD12Q land cover map for discussion in this study.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Global pattern of total, direct, and indirect effects of land cover change since pre-industry on ecosystem sensitivity to droughts (a) and Tmax (b), estimated from the dynamic global vegetation model mean. The averages of the effects on ecosystem sensitivity to droughts (c) and ecosystem sensitivity to Tmax (d) over the deforested areas are shown for all DGVMs, with error bars indicating the standard deviation.</figcaption></figure>

3.4. Direct and Indirect Effects From Historical LULCC

The effects of vegetation changes shown in Figure 5 point to potential substantial effects of historical LULCC on ecosystem sensitivity to drought and temperature. However, the space-for-time substitution method applied to satellite products is limited to identifying only the potential direct effects of LULCC. DGVMs, on the other hand, also allow for the inference of indirect LULCC effects arising from changes in α or β for the same vegetation types within a grid cell in the simulation with LULCC compared to the simulation with static pre-industrial land cover maps. By contrasting the reconstructed ā and ß between S3 and S2 simulations, we evaluate the total historical LULCC effects on ecosystem sensitivity to droughts and temperature. The total global LULCC effects on the model ensemble mean ā and ß are shown in Figures 6a and 6b. We estimate the uncertainty range of the total, direct and indirect historical LULCC effects by plus/minus one standard deviation. The global mean effect on ā is 0.008 ± 0.04 per drought month and the global average effect on ß is 0.013 ± 0.06. However, decreased ecosystem sensitivity to drought and temperature is simulated in deforestation hotspots where forests are converted to grasslands, including Southern Brazil, Western US, and Australia.

We further separate the total effects into direct and indirect effects as outlined in Equation 8 to understand their
relative importance. Direct effects on ā estimated by the model ensemble mean exhibit a global average of around
0 ± 0.02 per drought month. 52% of the areas show negative direct LULCC effects. They are negative in some
deforested areas where forests are dominantly converted to grasslands at the edge of the Amazon, in Southeast
Asia, Central America, and Australia but positive in areas with higher forest cover in 2001–2020 compared to
1700, such as Western Europe (Figures 6a and S9a in Supporting Information S1). Indirect effects, on the
contrary, contribute to increasing ā with a global average of 0.008 ± 0.04 per drought month. 57% of the areas
show positive indirect effects. Specifically, in contrast to negative direct effects estimated at the edge of the
Amazon, indirect effects are largely positive in the Amazon. The mean direct effect on ā in Brazil is
−0.013 ± 0.026 per drought month, but the mean indirect effect is 0.018 ± 0.053 per drought month. Among
different DGVMs, mean direct effects are negative in CABLE-POP and OCN but indirect effects contribute
positively to the total LULCC effects for all models in deforested areas (Figure 6c).
The direct effects on $\hat{\beta}$ exhibit a mean of around 0 $\pm$ 0.04 and they are negative in 45% of the areas, particularly in deforested areas at the edge of the Amazon, in South Asia, Western US, and part of Eastern Europe (Figures 6b and S9b in Supporting Information S1). The indirect effects contribute to an increased $\hat{\beta}$ (Figure S9b in Supporting Information S1). The mean indirect effect on $\hat{\beta}$ is 0.013 $\pm$ 0.06 and they are positive in 55% of the areas. Mean direct effects on $\hat{\beta}$ are negative in three DGVMs (CABLE-POP, OCN, and ORCHIDEE), and mean indirect effects are positive in all models except ORCHIDEE in deforested areas. The two effects compensate for each other in CABLE-POP and OCN.

In summary, DGVMs suggest that while direct effects result in less positive or more negative ecosystem sensitivity (weaker resistance) to drought and temperature over many deforested areas, due to the impacts of the predominant transition from forests to grasslands in 2001–2020 compared to 1700, indirect effects remain important and compensate for direct effects in many regions.

# 4. Discussions

## 4.1. Contrasting Ecosystem Sensitivity and Resistance to Drought and Temperature Between Forests and Non-Forest Vegetations

In this study, we assess kNDVI sensitivity to drought and temperature. Drought and temperature sensitivity can be positive or negative. Positive responses during drought periods are also observed in some ecosystems, reflected in increased GPP (Bastos et al., 2020) or vegetation greenness (Xiao et al., 2023). Negative values of sensitivities of the vegetation state to drought and temperature reflect ecosystem resistance, which measures the ability of ecosystems to persist and maintain their functioning during a disturbance (Gessler et al., 2020). More negative sensitivities indicate weaker resistance. Forests exhibit higher positive or less negative sensitivity (i.e., stronger resistance) to droughts than non-forests in more than 63% of the areas in all climate zones. This finding is supported by previous experimental studies (Anderegg et al., 2018), and research using GPP data sets across Europe (Bastos et al., 2020; Zhang et al., 2016), SIF (Walther et al., 2019), and L-VOD and EVI globally (Xiao et al., 2023). Higher resistance to droughts in forests compared to non-forest vegetations can be explained by the intrinsic structural and physiological differences between trees and crops, such as deeper rooting depth of trees (Canadell et al., 1996), greater water storage capacity in forest stems (Matheny et al., 2015), a weaker decline in light use efficiency with SM for forests compared to non-forest vegetation (Walther et al., 2019) and distinct water use strategies in response to drought between forest and grassland or cropland (Buras et al., 2020; Fu et al., 2020; Teuling et al., 2010; Wolf et al., 2013).

Similarly, we find stronger positive sensitivity and less negative sensitivity to higher temperatures in forests than in non-forests from kNDVI, particularly in areas where forests benefit from higher temperatures ($\beta$ > 0). In temperate regions, forests might benefit from higher temperatures through warming-induced phenological changes, while grasses and crops may exhibit a negative response in photosynthesis due to a nonlinear response of photosynthesis to temperature, resulting in a less positive or more negative sensitivity (weaker resistance) to hot extreme days. A notable trend toward earlier spring growing season onset and delayed autumn senescence in the temperate forests in the eastern USA has been observed (Keenan et al., 2014). In Russia during the heatwave in 2010, strong reductions in photosynthesis and greenness in grasslands and croplands were found due to heat effects on canopy structure, while regions of predominantly mixed forests showed small to insignificant reductions with even enhanced absorbed radiation (Yoshida et al., 2015).

Our estimates of ecosystem sensitivity to drought and temperature derived from kNDVI may be affected by uncertainties in the original NDVI data set, MOD13C1 v061. This data set employs a QA filtering scheme to address cloud contamination when processing input 1 km pixels into a 0.05° geographic product (Didan & Munoz, 2019). This approach reduces bias toward sunny days in annual mean kNDVI values (Didan & Munoz, 2019).

Severe tree mortality events following a strong drought year can impact ecosystem carbon fluxes in subsequent years (Yu et al., 2022), potentially influencing kNDVI and its estimated sensitivity to drought and temperature. However, the impact of tree mortality on vegetation greenness may be partially or entirely obscured at coarser spatial scales, such as 250 m or 8 km (Yan et al., 2024). In addition, the lag-1 term in Equation 2 accounts for legacy effects from disturbances in the previous year, making our estimation of kNDVI sensitivity to drought and temperature less affected by long-term tree mortality events.
Apart from vegetation types, land management practices, such as irrigation and wood harvest, can significantly influence ecosystem sensitivity to drought and temperature. However, the absence of explicit representation of irrigation in the TRENDY v11 DGVMs used in this study limits our ability to directly assess the effects of irrigation in these models. Future incorporation of irrigation processes in those DGVMs without explicit representation of irrigation, could provide more insights into the role of land management in shaping ecosystem resistance to climate extremes (Pongratz et al., 2021).

## 4.2. Comparison Between Observed and Modeled Transition Effects

Less positive or more negative sensitivity to drought $\alpha$, when forests are converted to grasslands, is evident from kNDVI. This effect is also simulated by the DGVM ensemble mean. However, when forests are converted to croplands, in contrast to kNDVI, CLASSIC and SDGVM predict a less negative sensitivity (higher resistance) in croplands than in forests. This discrepancy may stem from diverse representations of complex crop management and varying performances of DGVMs in agricultural lands. SDGVM tends to underpredict $\alpha$ across Earth's major agricultural regions, probably attributable to the absence of agricultural management practices such as improved seeds, fertilization, and irrigation (Walker et al., 2017). In CLASSIC, the incorporation of an N cycle in CLASSIC is in development and there is no representation of fertilization and irrigation in croplands. Furthermore, the benchmark for cropland is limited to one FLUXNET site (Melton et al., 2020). All five DGVMs also do not simulate irrigation explicitly (Table S3 in Supporting Information S1) and might assume no or less water stress for crops, which might explain less negative drought sensitivity of crops than trees.

Land cover transition effects on $\beta$ estimated from kNDVI are generally negative, indicating a decreased sensitivity to temperature when forests are converted to grasslands. Models do not simulate these negative effects on $\beta$, particularly at around 50°N and in the central Amazon (Figures 4b and 5b). The main reason for this discrepancy between modeled and observed effects lies in the underestimation of $\beta_{DGVMS_f}$ in these regions (Figures S8c and S8d in Supporting Information S1). Although both models and kNDVI exhibit positive sensitivity to temperature in forests in mid-latitudes, suggesting a beneficial effect of higher temperatures, the model mean underestimates this effect. In the tropics, $\beta_{kNDVI_f}$ estimated from kNDVI is positive in 61% of the area, whereas $\beta_{DGVMS_f}$ estimated from the model mean is only positive in 29% of the area, also contributing to the underestimated transition effects between forests and grasslands.

Another potential explanation for the discrepancies, aside from the suboptimal performance of DGVMs in simulating crops and vegetation sensitivity to temperature, is the inherent differences between vegetation carbon and kNDVI. The latter measures only the variations in greenness of the vegetation canopy. While there is a good linear correlation between L-VOD, a proxy for vegetation biomass, and kNDVI in ecosystems with relatively low biomass density, such as grasslands and croplands, divergent signals appear in dense tropical forests where the vertical canopy structure becomes complex (Xiao et al., 2023). Therefore, the comparison between ecosystem sensitivity calculated from kNDVI and model-simulated vegetation carbon might yield less accurate results in ecosystems with dense biomass, such as tropical forests. The choice of different vegetation indices, such as EVI or NIRv, might also influence the results because they reflect slightly different vegetation properties.

A few uncertainties exist in the estimated vegetation transition effects using the space-for-time substitution. First, we do not consider spatial gradients in general soil properties within the local moving window of 0.25° when applying the space-for-time substitution and we assume that vegetation cover is the only driver of the changes in ecosystem sensitivity to drought and temperature in the local moving window. Second, we cautiously filter our vegetation transition effects retrieved from kNDVI and the space-for-time substitution and only keep around 3% of the total pixels with valid kNDVI sensitivity to drought and temperature. The filtered results are sparse in certain regions, particularly those with large-scale homogeneous vegetation types and low land cover co-occurrence, such as boreal regions dominated by shrubs and savannas, and the Central Amazon dominated by forests. Despite this, we still retained over 110,000, 190,000, and 170,000 pixels for the three transition types, respectively, as shown in Figures 4a and 4b. The high-resolution kNDVI and climate data set ensure homogeneous climate conditions within a moving window of 0.25°. They also provide enough pixels for analysis even after our strict filtering. We also recognize that the space-for-time substitution approach may not capture the potential differences between forest-to-grassland/cropland and grassland/cropland-to-forest transitions, whose ecological impacts can be asymmetric. Deforestation typically results in a significant loss of biodiversity, while afforestation practices often favor a few commercially successful tree species, leading to forests with less
structural complexity (Ma et al., 2023; Naudts et al., 2016). Deforestation may more severely impact vegetation's sensitivity to drought and high temperatures compared to the mitigating effects of afforestation.

### 4.3. Direct and Indirect Effects of Historical LULCC on Ecosystem Sensitivity to Drought

By contrasting model simulations S3 and S2, we compute the historical LULCC effects on drought sensitivity. Subsequently, we decompose the total effects into direct and indirect components. Direct effects contribute to a decreased drought sensitivity in deforested areas in CABLE-POP and OCN (Figure 6c), reflecting negative transition effects on drought sensitivity when forests are transformed into grasslands or croplands. In contrast, CLASSIC and SDGVM exhibit positive mean direct effects because they simulate a strong increase in $\alpha$ when forests are converted to croplands.

Nevertheless, all models show positive mean indirect effects in deforested areas, resulting in a more positive or less negative drought sensitivity under LULCC. These positive indirect effects are primarily contributed by the increased drought sensitivity of forests in these areas (Figures S10 and S11 in Supporting Information S1). Over 56% of the forests under identical CO2 and climate conditions exhibit less negative sensitivity to drought (higher drought resistance) in the S3 simulation with large-scale deforestation compared to the S2 simulation (Figure S11 in Supporting Information S1). Additionally, the difference in drought sensitivity between S3 and S2 shows a positive global mean (Figure S10 in Supporting Information S1). Such an increased resistance has also been reported in previous observation-based studies and can be partly attributed to reduced competition for resources within forests. Higher resource acquisition capacity per tree with increasing growing space has been found to reduce drought-induced mortality (Allen et al., 2010; N. McDowell et al., 2008). Forest thinning has been assessed as an effective management strategy to mitigate drought impacts (D'Amato et al., 2013; Gavinet et al., 2020; Sohn et al., 2016), especially on sites where water is the main growth-limiting factor (Sohn et al., 2016). Forest stands with less dense canopies following thinning are often characterized by higher soil water availability for the residual trees (Brooks & Mitchell, 2011; McDowell et al., 2003).

DGVMs, however, are not capable of detecting the edge effects and might overestimate the beneficial indirect effects on ecosystem sensitivity to drought. In reality, deforestation also leads to higher forest fragmentation and an increased proportion of forests being located in close proximity to the forest edge. Forest edges alter microclimate conditions, such as temperature, humidity, radiation, and wind speed, by exposing parts of the forest environment to external climatic conditions, reducing the ability of a forest to buffer its internal microclimate from those more extreme macroclimate conditions (Ewers & Banks-Leite, 2013). Edge effects have also been reported to increase wildfire susceptibility (Alencar et al., 2004; Cochrane & Laurance, 2002) and tree mortality (Laurance et al., 2000; Qie et al., 2017).

### 5. Conclusions

In conclusion, our study combines MODIS kNDVI data and TRENDY-v11 DGVMs to conduct a comprehensive assessment of the effects of vegetation cover transitions and historical LULCC on ecosystem sensitivity to both droughts and maximum temperatures at the global scale. Our findings underscore less positive or more negative sensitivity (a decrease in resistance) to drought and temperature when forests are converted to grasslands and croplands from kNDVI. Despite large variability, most of the DGVMs agree on higher positive or less negative sensitivity (higher drought resistance) of forests than grasslands. We find a less negative or more positive drought sensitivity due to LULCC from DGVMs. By distinguishing between the direct and indirect effects of LULCC on ecosystem sensitivity, our study highlights the significant role of indirect effects, which arise from changes in the sensitivity of the same vegetation type, especially forests. These results are based on DGVMs, which are not coupled to the atmosphere, the atmospheric feedback induced by LULCC should also be considered in the future to comprehensively assess the vulnerability of ecosystems to climate extremes under climate change. Overall, our findings emphasize the critical role of forests in maintaining ecosystem stability and resistance. Understanding the complex effects of vegetation transition aids in developing land management and conservation strategies that enhance ecosystem stability under ongoing climate change. Our study emphasizes the need to improve vegetation sensitivity representation in DGVMs to better represent emergent changes in vegetation resistance due to interactions between human activities and climate change.
## Data Availability Statement

The MODIS NDVI data set is available from https://lpdaac.usgs.gov/products/mod13c1v061/ (Didan, 2021). The MODIS MCD12Q1 v061 land cover data set is available from https://lpdaac.usgs.gov/products/mcd12q1v061/ (Friedl & Sulla-Menashe, 2022). The TerraClimate data set (Abatzoglou et al., 2018) can be accessed from https://www.climatologylab.org/terraclimate.html. The TRENDY-v11 data set (Sitch et al., 2024) can be accessed from https://globalcarbonbudgetdata.org/closed-access-requests.html upon request. The Julia YAXArrays Toolbox package used for space-for-time substitution in this study is available at https://zenodo.org/records/7989936 (Pabon-Moreno et al., 2023). The code of this study is available at https://doi.org/10.5281/zenodo.15366417 (Xiao, 2025).

## Acknowledgments

Chenwei Xiao acknowledges support from the International Max Planck Research School for Global Biogeochemical Cycles. The authors would like to thank Ulrich Weber for the preparation of TerraClimate data, MODIS land cover, MODIS NDVI, and ESA CCI PFT maps. The authors thank Dr. Peter Anthoni for the insightful discussions related to LPJ-GUESS simulations. ORNL is managed by UT-Battelle LLC for the DOE under contract DE-AC05-1008 00R22725. Sonke Zaehle acknowledges funding by the European Commission H2020 research project Carbon-Climate Interaction in the Coming Century (4C) under Grant 821003. Ana Bastos acknowledges funding by the European Union (ERC StG, ForExD, Grant 101039567). Gregory Duveiller acknowledges funding by the European Research Council (ERC) Synergy Grant 'Understanding and modelling the Earth System with Machine Learning (USMLE)' under the Horizon 2020 research and innovation programme (Grant Agreement No. 855187). Gregory Duveiller and Daniel E. Pabon-Moreno acknowledge funding from the European Union's Horizon Europe Research and Innovation programme under grant agreement no. 101059548 (OEMC project). Open Access funding enabled and organized by Projekt DEAL.

## References

Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A., & Hegewisch, K. C. (2018). TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. *Scientific Data*, 5(1), 170191. Article 1. https://doi.org/10.1038/sdata.2017.191

Alencar, A. A. C., Solórzano, L. A., & Nepstad, D. C. (2004). Modeling forest understory fires in an eastern Amazonian landscape. *Ecological Applications*, 14(sp4), 139–149. https://doi.org/10.1890/01-0629

Allen, C. D., Macalady, A. K., Chenchouni, H., Bachelet, D., McDowell, N., Vennetier, M., et al. (2010). A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. *Forest Ecology and Management*, 259(4), 660–684. https://doi.org/10.1016/j.foreco.2009.09.001

Anderegg, W. R. L., Konings, A. G., Trugman, A. T., Yu, K., Bowling, D. R., Gabbitas, R., et al. (2018). Hydraulic diversity of forests regulates ecosystem resilience during drought. *Nature*, 561(7724), 538–541. Article 7724. https://doi.org/10.1038/s41586-018-0539-7

Bastos, A., Fu, Z., Ciais, P., Friedlingstein, P., Sitch, S., Pongratz, J., et al. (2020). Impacts of extreme summers on European ecosystems: A comparative analysis of 2003, 2010 and 2018. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 375(1810), 20190507. https://doi.org/10.1098/rstb.2019.0507

Bottero, A., D'Amato, A. W., Palik, B. J., Bradford, J. B., Fraver, S., Battaglia, M. A., & Asherin, L. A. (2017). Density-dependent vulnerability of forest ecosystems to drought. *Journal of Applied Ecology*, 54(6), 1605–1614. https://doi.org/10.1111/1365-2664.12847

Bourgoin, C., Ceccherini, G., Girardello, M., Vancutsem, C., Avitabile, V., Beck, P. S. A., et al. (2024). Human degradation of tropical moist forests is greater than previously estimated. *Nature*, 1–7(8021), 570–576. https://doi.org/10.1038/s41586-024-07629-0

Brooks, J. R., & Mitchell, A. K. (2011). Interpreting tree responses to thinning and fertilization using tree-ring stable isotopes. *New Phytologist*, 190(3), 770–782. https://doi.org/10.1111/j.1469-8137.2010.03627.x

Buras, A., Rammig, A., & Zang, C. S. (2020). Quantifying impacts of the 2018 drought on European ecosystems in comparison to 2003. *Bio-geosciences*, 17(6), 1655–1672. https://doi.org/10.5194/bg-17-1655-2020

Camps-Valls, G., Campos-Taberner, M., Moreno-Martínez, A., Walther, S., Duveiller, G., Cescatti, A., et al. (2021). A unified vegetation index for quantifying the terrestrial biosphere. *Science Advances*, 7(9), eabc7447. https://doi.org/10.1126/sciadv.abc7447

Canadell, J., Jackson, R. B., Ehleringer, J. B., Mooney, H. A., Sala, O. E., & Schulze, E.-D. (1996). Maximum rooting depth of vegetation types at the global scale. *Oecologia*, 108(4), 583–595. https://doi.org/10.1007/BF00329030

Chini, L., Hurtt, G., Sahajpal, R., Frolking, S., Klein Goldewijk, K., Sitch, S., et al. (2021). Land-use harmonization datasets for annual global carbon budgets. *Earth System Science Data*, 13(8), 4175–4189. https://doi.org/10.5194/essd-13-4175-2021

Ciais, P., Reichstein, M., Viovy, N., Granier, A., Ogée, J., Allard, V., et al. (2005). Europe-wide reduction in primary productivity caused by the heat and drought in 2003. *Nature*, 437(7058), 529–533. Article 7058. https://doi.org/10.1038/nature03972

Cochrane, M. A., & Laurance, W. F. (2002). Fire as a large-scale edge effect in Amazonian forests. *Journal of Tropical Ecology*, 18(3), 311–325. https://doi.org/10.1017/s0266476402002237

Condit, R., Hubbell, S., & Foster, R. (1995). Mortality-rates of 205 neotropical tree and shrub species and the impact of a severe drought. *Ecological Monographs*, 65(4), 419–439. https://doi.org/10.2307/2963497

Cranko Page, J., De Kauwe, M. G., Abramowitz, G., Cleverly, J., Hinko-Najera, N., Hovenden, M. J., et al. (2022). Examining the role of environmental memory in the predictability of carbon and water fluxes across Australian ecosystems. *Biogeosciences*, 19(7), 1913–1932. https://doi.org/10.5194/bg-19-1913-2022

D'Amato, A. W., Bradford, J. B., Fraver, S., & Palik, B. J. (2013). Effects of thinning on drought vulnerability and climate response in north temperate forest ecosystems. *Ecological Applications*, 23(8), 1735–1742. https://doi.org/10.1890/13-0677.1

De Keersmaecker, W., van Rooijen, N., Lhermitte, S., Tit's, L., Schaminée, J., Coppin, P., et al. (2016). Species-rich semi-natural grasslands have a higher resistance but a lower resilience than intensively managed agricultural grasslands in response to climate anomalies. *Journal of Applied Ecology*, 53(2), 430–439. https://doi.org/10.1111/1365-2664.12595

Denissen, J. M. C., Teuling, A. J., Pitman, A. J., Koirala, S., Migliavacca, M., Li, W., et al. (2022). Widespread shift from ecosystem energy to water limitation with climate change. *Nature Climate Change*, 12(7), 677–684. Article 7. https://doi.org/10.1038/s41558-022-01403-8

DeSoto, L., Cailleret, M., Sterck, F., Jansen, S., Kramer, K., Robert, E. M. R., et al. (2020). Low growth resilience to drought is related to future mortality risk in trees. *Nature Communications*, 11(1), 545. https://doi.org/10.1038/s41467-020-14300-5

Devaraju, N., de Noblet-Ducoudré, N., Quesada, B., & Bala, G. (2018). Quantifying the relative importance of direct and indirect biophysical effects of deforestation on surface temperature and teleconnections. *Journal of Climate*, 31(10), 3811–3829. https://doi.org/10.1175/JCLI-D-17-0563.1

Didan, K. (2021). MODIS/Terra vegetation indices 16-day L3 global 0.05Deg CMG V061 [Dataset]. NASA EOSDIS Land Processes Distributed Active Archive Center. https://doi.org/10.5067/MODIS/MOD13C1.061

Didan, K., & Munoz, A. B. (2019). MODIS vegetation index user's guide (MOD13 series).

Duveiller, G., Filipponi, F., Ceglar, A., Bojanowski, J., Alkama, R., & Cescatti, A. (2021). Revealing the widespread potential of forests to increase low level cloud cover. *Nature Communications*, 12(1), 4337. https://doi.org/10.1038/s41467-021-24551-5

Duveiller, G., Hooker, J., & Cescatti, A. (2018). The mark of vegetation change on Earth's surface energy balance. *Nature Communications*, 9(1), 679. Article 1. https://doi.org/10.1038/s41467-017-02810-8

Ewers, R. M., & Banks-Leite, C. (2013). Fragmentation impairs the microclimate buffering effect of tropical forests. *PLoS One*, 8(3), e58093. https://doi.org/10.1371/journal.pone.0058093
Flach, M., Sippel, S., Gans, F., Bastos, A., Brenning, A., Reichstein, M., & Mahecha, M. D. (2018). Contrasting biosphere responses to hydrometeorological extremes: Revisiting the 2010 western Russian heatwave. *Biogeosciences*, *15*(20), 6067–6085. https://doi.org/10.5194/bg-15-6067-2018

Friedl, M., & Sulla-Menashe, D. (2022). MODIS/Terra+Aqua land cover type yearly L3 global 500m SIN grid V061 [Dataset]. NASA EOSDIS Land Processes Distributed Active Archive Center. https://doi.org/10.5067/MODIS/MCD12Q1.061

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Bakker, D. C. E., Hauck, J., et al. (2023). Global carbon budget 2023. *Earth System Science Data*, *15*(12), 5301–5369. https://doi.org/10.5194/essd-15-5301-2023

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Gregor, L., Hauck, J., et al. (2022). Global carbon budget 2022. *Earth System Science Data*, *14*(11), 4811–4900. https://doi.org/10.5194/essd-14-4811-2022

Fu, Z., Ciais, P., Bastos, A., Stoy, P. C., Yang, H., Green, J. K., et al. (2020). Sensitivity of gross primary productivity to climatic drivers during the summer drought of 2018 in Europe. *Philosophical Transactions of the Royal Society B: Biological Sciences*, *375*(1810), 20190747. https://doi.org/10.1098/rstb.2019.0747

Gavinet, J., Ourcival, J.-M., Gauzere, J., García de Jalón, L., & Limousin, J.-M. (2020). Drought mitigation by thinning: Benefits from the stem to the stand along 15 years of experimental rainfall exclusion in a holm oak coppice. *Forest Ecology and Management*, *473*, 118266. https://doi.org/10.1016/j.foreco.2020.118266

Gessler, A., Bottero, A., Marshall, J., & Arend, M. (2020). The way back: Recovery of trees from drought and its implication for acclimation. *New Phytologist*, *228*(6), 1704–1709. https://doi.org/10.1111/nph.16703

Grossiord, C., Granier, A., Ratcliffe, S., Bouriaud, O., Bruelheide, H., Chečko, E., et al. (2014). Tree diversity does not always improve resistance of forest ecosystems to drought. *Proceedings of the National Academy of Sciences*, *111*(41), 14812–14815. https://doi.org/10.1073/pnas.1411970111

Haddad, N. M., Brudvig, L. A., Clobert, J., Davies, K. F., Gonzalez, A., Holt, R. D., et al. (2015). Habitat fragmentation and its lasting impact on Earth's ecosystems. *Science Advances*, *1*(2), e1500052. https://doi.org/10.1126/sciadv.1500052

Harmeling, S., Dornhege, G., Tax, D., Meinecke, F., & Müller, K.-R. (2006). From outliers to prototypes: Ordering data. *Neurocomputing*, *69*(13), 1608–1618. https://doi.org/10.1016/j.neucom.2005.05.015

Harper, K. L., Lamarche, C., Hartley, A., Peylin, P., Ottlé, C., Bastrikov, V., et al. (2023). A 29-year time series of annual 300 m resolution plant-functional-type maps for climate models. *Earth System Science Data*, *15*(3), 1465–1499. https://doi.org/10.5194/essd-15-1465-2023

Harris, I., Jones, P. d., Osborn, T. j., & Lister, D. h. (2014). Updated high-resolution grids of monthly climatic observations – The CRU TS3.10 dataset. *International Journal of Climatology*, *34*(3), 623–642. https://doi.org/10.1002/joc.3711

Harris, I., Osborn, T. J., Jones, P., & Lister, D. (2020). Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. *Scientific Data*, *7*(1), 109. Article 1. https://doi.org/10.1038/s41597-020-0453-3

Hartley, A. J., MacBean, N., Georgievski, G., & Bontemps, S. (2017). Uncertainty in plant functional type distributions and its impact on land surface models. *Remote Sensing of Environment*, *203*, 71–89. https://doi.org/10.1016/j.rse.2017.07.037

Hurtt, G. C., Chini, L., Sahajpal, R., Frolking, S., Bodirsky, B. L., Calvin, K., et al. (2020). Harmonization of global land use change and management for the period 850–2100 (LUH2) for CMIP6. *Geoscientific Model Development*, *13*(11), 5425–5464. https://doi.org/10.5194/gmd-13-5425-2020

IPCC. (2021). *Climate change 2021 – The physical science basis: Working group I contribution to the sixth assessment report of the inter-governmental panel on climate change*. Cambridge University Press. https://doi.org/10.1017/9781009157896

Isbell, F., Craven, D., Connolly, J., Loreau, M., Schmid, B., Beierkuhnlein, C., et al. (2015). Biodiversity increases the resistance of ecosystem productivity to climate extremes. *Nature*, *526*(7574), 574–577. Article 7574. https://doi.org/10.1038/nature15374

Jiao, W., Wang, L., Smith, W. K., Chang, Q., Wang, H., & D'Odorico, P. (2021). Observed increasing water constraint on vegetation growth over the last three decades. *Nature Communications*, *12*(1), 3777. https://doi.org/10.1038/s41467-021-24016-9

Keenan, T. F., Gray, J., Friedl, M. A., Toomey, M., Bohrer, G., Hollinger, D. Y., et al. (2014). Net carbon uptake has increased through warming-induced changes in temperate forest phenology. *Nature Climate Change*, *4*(7), 598–604. Article 7. https://doi.org/10.1038/nclimate2253

Lan, X., Tans, P., & Thoning, K. W. (2024). Trends in globally-averaged CO2 determined from NOAA global monitoring laboratory measurements [Dataset]. *Version 2024-06*. https://doi.org/10.15138/9N0H-ZH07

Laurance, W. F., Delamôrica, P., Laurance, S. G., Vasconcelos, H. L., & Lovejoy, T. E. (2000). Rainforest fragmentation kills big trees. *Nature*, *404*(6780), 836. https://doi.org/10.1038/35009032

Le Quéré, C., Andres, R. J., Boden, T., Conway, T., Houghton, R. A., House, J. I., et al. (2013). The global carbon budget 1959–2011. *Earth System Science Data*, *5*(1), 165–185. https://doi.org/10.5194/essd-5-165-2013

Li, W., MacBean, N., Ciais, P., Defourny, P., Lamarche, C., Bontemps, S., et al. (2018). Gross and net land cover changes in the main plant functional types derived from the annual ESA CCI land cover maps (1992–2015). *Earth System Science Data*, *10*(1), 219–234. https://doi.org/10.5194/essd-10-219-2018

Li, X., Wigneron, J.-P., Frappart, F., Fan, L., Ciais, P., Fensholt, R., et al. (2021). Global-scale assessment and inter-comparison of recently developed/reprocessed microwave satellite vegetation optical depth products. *Remote Sensing of Environment*, *253*, 112208. https://doi.org/10.1016/j.rse.2020.112208

Li, Y., Zhao, M., Motesharrei, S., Mu, Q., Kalnay, E., & Li, S. (2015). Local cooling and warming effects of forests based on satellite observations. *Nature Communications*, *6*(1), 6603. https://doi.org/10.1038/ncomms7603

Liu, D., Wang, T., Peñuelas, J., & Piao, S. (2022). Drought resistance enhanced by tree species diversity in global forests. *Nature Geoscience*, *15*(10), 800–804. Article 10. https://doi.org/10.1038/s41561-022-01026-w

Liu, F., Liu, H., Xu, C., Shi, L., Zhu, X., Qi, Y., & He, W. (2021). Old-growth forests show low canopy resilience to droughts at the southern edge of the taiga. *Global Change Biology*, *27*(11), 2392–2402. https://doi.org/10.1111/gcb.15605

Liu, L., Chen, X., Ciais, P., Yuan, W., Maignan, F., Wu, J., et al. (2021). Tropical tall forests are more sensitive and vulnerable to drought than short forests. *Global Change Biology*, *28*, 1583–1595. https://doi.org/10.1111/gcb.16017

Liu, Y., Schwalm, C. R., Samuels-Crow, K. E., & Ogle, K. (2019). Ecological memory of daily carbon exchange across the globe and its importance in drylands. *Ecology Letters*, *22*(11), 1806–1816. https://doi.org/10.1111/ele.13363

Liu, Z., Ballantyne, A. P., & Cooper, L. A. (2019). Biophysical feedback of global forest fires on surface temperature. *Nature Communications*, *10*(1), 214. https://doi.org/10.1038/s41467-018-08237-z

Ma, Q., Su, Y., Niu, C., Ma, Q., Hu, T., Luo, X., et al. (2023). Tree mortality during long-term droughts is lower in structurally complex forest stands. *Nature Communications*, *14*(1), 7467. https://doi.org/10.1038/s41467-023-43083-8

Matheny, A. M., Bohrer, G., Garrity, S. R., Morin, T. H., Howard, C. J., & Vogel, C. S. (2015). Observations of stem water storage in trees of opposing hydraulic strategies. *Ecosphere*, *6*(9), art165-13. https://doi.org/10.1890/ES15-00170.1
Mcdowell, N., Brooks, J. R., Fitzgerald, S. A., & Bond, B. J. (2003). Carbon isotope discrimination and growth response of old Pinus ponderosa trees to stand density reductions. *Plant, Cell and Environment*, **26**(4), 631–644. https://doi.org/10.1046/j.1365-3040.2003.00999.x

McDowell, N., Pockman, W. T., Allen, C. D., Breshears, D. D., Cobb, N., Kolb, T., et al. (2008). Mechanisms of plant survival and mortality during drought: Why do some plants survive while others succumb to drought? *New Phytologist*, **178**(4), 719–739. https://doi.org/10.1111/j. 1469-8137.2008.02436.x

McDowell, N. G., & Allen, C. D. (2015). Darcy's law predicts widespread forest mortality under climate warming. *Nature Climate Change*, **5**(7), 669–672. Article 7. https://doi.org/10.1038/nclimate2641

Melton, J. R., Arora, V. K., Wisernig-Cojoc, E., Seiler, C., Fortier, M., Chan, E., & Teckentrup, L. (2020). CLASSIC v1.0: The open-source community successor to the Canadian Land Surface Scheme (CLASS) and the Canadian Terrestrial Ecosystem Model (CTEM) – Part 1: Model framework and site-level performance. *Geoscientific Model Development*, **13**(6), 2825–2850. https://doi.org/10.5194/gmd-13-2825- 2020

Naudts, K., Chen, Y., McGrath, M. J., Ryder, J., Valade, A., Otto, J., & Luyssaert, S. (2016). Europe's forest management did not mitigate climate warming. *Science*, **351**(6273), 597–600. https://doi.org/10.1126/science.aad7270

Nepstad, D. C., Tohver, I. M., Ray, D., Moutinho, P., & Cardinot, G. (2007). Mortality of large trees and lianas following experimental drought in an Amazon forest. *Ecology*, **88**(9), 2259–2269. https://doi.org/10.1890/06-1046.1

Pabon-Moreno, D. E., Duveiller, G., Reichstein, M., Gans, F., Cremer, F., & Winkler, A. (2023). dpabon/YAXArraysToolbox.jl: Zenodo publication (Version v0.3.2). [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.7989936

Parzen, E. (1962). On estimation of a probability density function and mode. *The Annals of Mathematical Statistics*, **33**(3), 1065–1076. https://doi. org/10.1214/aoms/1177704472

Phillips, O. L., Aragão, L. E. O. C., Lewis, S. L., Fisher, J. B., Lloyd, J., López-González, G., et al. (2009). Drought sensitivity of the Amazon rainforest. *Science*, **323**(5919), 1344–1347. https://doi.org/10.1126/science.1164033

Picchio, R., Mederski, P. S., & Tavankar, F. (2020). How and how much, do harvesting activities affect forest soil, regeneration and stands? *Current Forestry Reports*, **6**(2), 115–128. https://doi.org/10.1007/s40725-020-00113-8

Pongratz, J., Schwingshackl, C., Bultan, S., Obermeier, W., Havermann, F., & Guo, S. (2021). Land use effects on climate: Current state, recent progress, and emerging topics. *Current Climate Change Reports*, **7**(4), 99–120. https://doi.org/10.1007/s40641-021-00178-y

Qie, L., Lewis, S. L., Sullivan, M. J. P., Lopez-Gonzalez, G., Pickavance, G. C., Sunderland, T., et al. (2017). Long-term carbon sink in Borneo's forests halted by drought and vulnerable to edge effects. *Nature Communications*, **8**(1), 1966. https://doi.org/10.1038/s41467-017-01997-0

Qin, Y., Xiao, X., Wigneron, J.-P., Ciais, P., Brandt, M., Fan, L., et al. (2021). Carbon loss from forest degradation exceeds that from deforestation in the Brazilian Amazon. *Nature Climate Change*, **11**(5), 442–448. Article 5. https://doi.org/10.1038/s41558-021-01026-5

Seneviratne, S. I., Nicholls, N., Easterling, D., Goodess, C. M., Kanae, S., Kossin, J., et al. (2012). Changes in climate extremes and their impacts on the natural physical environment. In C. B. Field, Q. Dahe, T. F. Stocker, & V. Barros (Eds.), *Managing the risks of extreme events and disasters to advance climate change adaptation: Special report of the intergovernmental panel on climate change* (pp. 109–230). Cambridge University Press. https://doi.org/10.1017/CBO9781139177245.006

Sitch, S., O'Sullivan, M., Robertson, E., Friedlingstein, P., Albergel, C., Anthoni, P., et al. (2024). Trends and drivers of terrestrial sources and sinks of carbon dioxide: An overview of the TRENDY project. *Global Biogeochemical Cycles*, **38**(7), e2024GB008102. https://doi.org/10. 1029/2024GB008102

Sohn, J. A., Saha, S., & Bauhus, J. (2016). Potential of forest thinning to mitigate drought stress: A meta-analysis. *Forest Ecology and Management*, **380**, 261–273. https://doi.org/10.1016/j.foreco.2016.07.046

Teuling, A. J., Seneviratne, S. I., Stöckli, R., Reichstein, M., Moors, E., Ciais, P., et al. (2010). Contrasting response of European forest and grassland energy exchange to heatwaves. *Nature Geoscience*, **3**(10), 722–727. https://doi.org/10.1038/ngeo950

Walker, A. P., Quaife, T., van Bodegom, P. M., De Kauwe, M. G., Keenan, T. F., Joiner, J., et al. (2017). The impact of alternative trait-scaling hypotheses for the maximum photosynthetic carboxylation rate (Vcmax) on global gross primary production. *New Phytologist*, **215**(4), 1370– 1386. https://doi.org/10.1111/nph.14623

Walther, S., Duveiller, G., Jung, M., Guanter, L., Cescatti, A., & Camps-Valls, G. (2019). Satellite observations of the contrasting response of trees and grasses to variations in water availability. *Geophysical Research Letters*, **46**(3), 1429–1440. https://doi.org/10.1029/2018GL080535

Wang, Q., Moreno-Martínez, Á., Muñoz-Mari, J., Campos-Taberner, M., & Camps-Valls, G. (2023). Estimation of vegetation traits with kernel NDVI. *ISPRS Journal of Photogrammetry and Remote Sensing*, **195**, 408–417. https://doi.org/10.1016/j.isprsjprs.2022.12.019

Wilcoxon, F. (1945). Individual comparisons by ranking methods. *Biometric Bulletin*, **1**(6), 80–83. https://doi.org/10.2307/3001968

Winkler, K., Yang, H., Ganzenmüller, R., Fuchs, R., Ceccherini, G., Duveiller, G., et al. (2023). Changes in land use and management led to a decline in Eastern Europe's terrestrial carbon sink. *Communications Earth & Environment*, **4**(1), 237. Article 1. https://doi.org/10.1038/ s43247-023-00893-4

Wolf, S., Eugster, W., Ammann, C., Häni, M., Zielis, S., Hiller, R., et al. (2013). Contrasting response of grassland versus forest carbon and water fluxes to spring drought in Switzerland. *Environmental Research Letters*, **8**(3), 035007. https://doi.org/10.1088/1748-9326/8/3/035007

Xiao, C. (2025). Deforestation increases vegetation vulnerability to drought across biomes [Software]. Zenodo. https://doi.org/10.5281/zenodo. 15366417

Xiao, C., Zaehle, S., Yang, H., Wigneron, J.-P., Schmullius, C., & Bastos, A. (2023). Land cover and management effects on ecosystem resistance to drought stress. *Earth System Dynamics*, **14**(6), 1211–1237. https://doi.org/10.5194/esd-14-1211-2023

Yan, Y., Piao, S., Hammond, W. M., Chen, A., Hong, S., Xu, H., et al. (2024). Climate-induced tree-mortality pulses are obscured by broad-scale and long-term greening. *Nature Ecology & Evolution*, **8**(5), 1–12. https://doi.org/10.1038/s41559-024-02372-1

Yoshida, Y., Joiner, J., Tucker, C., Berry, J., Lee, J.-E., Walker, G., et al. (2015). The 2010 Russian drought impact on satellite measurements of solar-induced chlorophyll fluorescence: Insights from modeling and comparisons with parameters derived from satellite reflectances. *Remote Sensing of Environment*, **166**, 163–177. https://doi.org/10.1016/j.rse.2015.06.008

Yu, X., Orth, R., Reichstein, M., Bahn, M., Klosterhalfen, A., Knohl, A., et al. (2022). Contrasting drought legacy effects on gross primary productivity in a mixed versus pure beech forest. *Biogeosciences*, **19**(17), 4315–4329. https://doi.org/10.5194/bg-19-4315-2022

Zamora-Pereira, J. C., Yousefpour, R., Cailleret, M., Bugmann, H., & Hanewinkel, M. (2021). Magnitude and timing of density reduction are key for the resilience to severe drought in conifer-broadleaf mixed forests in Central Europe. *Annals of Forest Science*, **78**(3), 68. https://doi.org/10. 1007/s13595-021-01085-w

Zhang, Y., Xiao, X., Zhou, S., Ciais, P., McCarthy, H., & Luo, Y. (2016). Canopy and physiological controls of GPP during drought and heat wave. *Geophysical Research Letters*, **43**(7), 3325–3333. https://doi.org/10.1002/2016GL068501
## References From the Supporting Information

Asaadi, A., Arora, V. K., Melton, J. R., & Bartlett, P. (2018). An improved parameterization of leaf area index (LAI) seasonality in the Canadian Land Surface Scheme (CLASS) and Canadian Terrestrial Ecosystem Model (CTEM) modelling framework. *Biogeosciences*, *15*(22), 6885-6907. https://doi.org/10.5194/bg-15-6885-2018

de Rosnay, P., & Polcher, J. (1998). Modelling root water uptake in a complex land surface scheme coupled to a GCM. *Hydrology and Earth System Sciences*, *2*(2/3), 239-255. https://doi.org/10.5194/hess-2-239-1998

Haverd, V., Smith, B., Nieradzik, L., Briggs, P. R., Woodgate, W., Trudinger, C. M., et al. (2018). A new version of the CABLE land surface model (Subversion revision r4601) incorporating land use and land cover change, woody vegetation demography, and a novel optimisation-based approach to plant coordination of photosynthesis. *Geoscientific Model Development*, *11*(7), 2995-3026. https://doi.org/10.5194/gmd-11-2995-2018

Haverd, V., Smith, B., Nieradzik, L. P., & Briggs, P. R. (2014). A stand-alone tree demography and landscape structure module for earth system models: Integration with inventory data from temperate and boreal forests. *Biogeosciences*, *11*(15), 4039-4055. https://doi.org/10.5194/bg-11-4039-2014

Haberstroh, S., Werner, C., Grün, M., Kreuzwieser, J., Seifert, T., Schindler, D., & Christen, A. (2022). Central European 2018 hot drought shifts scots pine forest to its tipping point. *Plant Biology*, *24*(7), 1186-1197. https://doi.org/10.1111/plb.13455

Kattge, J., & Knorr, W. (2007). Temperature acclimation in a biochemical model of photosynthesis: A reanalysis of data from 36 species. *Plant, Cell and Environment*, *30*(9), 1176-1190. https://doi.org/10.1111/j.1365-3040.2007.01690.x

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogee, J., Polcher, J., Friedlingstein, P., et al. (2005). A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Global Biogeochemical Cycles*, *19*(1). https://doi.org/10.1029/2003GB002199

Kumarathunge, D. P., Medlyn, B. E., Drake, J. E., Tjoelker, M. G., Aspinwall, M. J., Battaglia, M., et al. (2019). Acclimation and adaptation components of the temperature dependence of plant photosynthesis at the global scale. *New Phytologist*, *222*(2), 768-784. https://doi.org/10.1111/nph.15668

Medlyn, B. E., Duursma, R. A., Eamus, D., Ellsworth, D. S., Prentice, I. C., Barton, C. V. M., et al. (2011). Reconciling the optimal and empirical approaches to modelling stomatal conductance. *Global Change Biology*, *17*(6), 2134-2144. https://doi.org/10.1111/j.1365-2486.2010.02375.x

Meyerholt, J., Zaehle, S., & Smith, M. J. (2016). Variability of projected terrestrial biosphere responses to elevated levels of atmospheric CO₂ due to uncertainty in biological nitrogen fixation. *Biogeosciences*, *13*(5), 1491-1518. https://doi.org/10.5194/bg-13-1491-2016

Melton, J. R., & Arora, V. K. (2016). Competition between plant functional types in the Canadian Terrestrial Ecosystem Model (CTEM) v. 2.0. *Geoscientific Model Development*, *9*(1), 323-361. https://doi.org/10.5194/gmd-9-323-2016

Picard, G., Quegan, S., Delbart, N., Lomas, M. R., Le Toan, T., & Woodward, F. I. (2005). Bud-burst modelling in Siberia and its impact on quantifying the carbon budget. *Global Change Biology*, *11*(12), 2164-2176. https://doi.org/10.1111/j.1365-2486.2005.01055.x

Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., et al. (2003). Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. *Global Change Biology*, *9*(2), 161-185. https://doi.org/10.1046/j.1365-2486.2003.00569.x

Thonicke, K., Venevsky, S., Sitch, S., & Cramer, W. (2001). The role of fire disturbance for global vegetation dynamics: Coupling fire into a dynamic global vegetation model. *Global Ecology and Biogeography*, *10*(6), 661-677. https://doi.org/10.1046/j.1466-822X.2001.00175.x

verseghy, D. (2017). CLASS – The Canadian land surface scheme. Climate Research Division, Science and Technology Branch, Environment Canada.

Vuichard, N., Messina, P., Luyssaert, S., Guenet, B., Zaehle, S., Ghattas, J., et al. (2019). Accounting for carbon and nitrogen interactions in the global terrestrial ecosystem model ORCHIDEE (trunk version, rev 4999): Multi-scale evaluation of gross primary production. *Geoscientific Model Development*, *12*(11), 4751-4779. https://doi.org/10.5194/gmd-12-4751-2019

Wang, Y. P., Law, R. M., & Pak, B. (2010). A global model of carbon, nitrogen and phosphorus cycles for the terrestrial biosphere. *Bio-geosciences*, *7*(7), 2261-2282. https://doi.org/10.5194/bg-7-2261-2010

Wang, Y. P., & Leuning, R. (1998). A two-leaf model for canopy conductance, photosynthesis and partitioning of available energy I: Model description and comparison with a multi-layered model. *Agricultural and Forest Meteorology*, *91*(1), 89-111. https://doi.org/10.1016/S0168-1923(98)00061-6

Woodward, F. I., & Lomas, M. R. (2004). Vegetation dynamics – Simulating responses to climatic change. *Biological Reviews*, *79*(3), 643-670. https://doi.org/10.1017/S1464793103006419

Woodward, F. I., Smith, T. M., & Emanuel, W. R. (1995). A global land primary productivity and phytogeography model. *Global Biogeochemical Cycles*, *9*(4), 471-490. https://doi.org/10.1029/95GB02432

Zaehle, S., Ciais, P., Friend, A. D., & Prieur, V. (2011). Carbon benefits of anthropogenic reactive nitrogen offset by nitrous oxide emissions. *Nature Geoscience*, *4*(9), 601-605. Article 9. https://doi.org/10.1038/ngeo1207

Zaehle, S., Friedlingstein, P., & Friend, A. D. (2010). Terrestrial nitrogen feedbacks may accelerate future climate change. *Geophysical Research Letters*, *37*(1). https://doi.org/10.1029/2009GL041345

Zaehle, S., & Friend, A. D. (2010). Carbon and nitrogen cycle dynamics in the O-CN land surface model: 1. Model description, site-scale evaluation, and sensitivity to parameter estimates. *Global Biogeochemical Cycles*, *24*(1). https://doi.org/10.1029/2009GB003521