Emerging signals of declining forest resilience under climate change

Abstract

Forest ecosystems depend on their capacity to withstand and recover from natural and anthropogenic perturbations (that is, their resilience)¹. Experimental evidence of sudden increases in tree mortality is raising concerns about variation in forest resilience², yet little is known about how it is evolving in response to climate change. Here we integrate satellite-based vegetation indices with machine learning to show how forest resilience, quantified in terms of critical slowing down indicators^3-5, has changed during the period 2000-2020. We show that tropical, arid and temperate forests are experiencing a significant decline in resilience, probably related to increased water limitations and climate variability. By contrast, boreal forests show divergent local patterns with an average increasing trend in resilience, probably benefiting from warming and CO₂ fertilization, which may outweigh the adverse effects of climate change. These patterns emerge consistently in both managed and intact forests, corroborating the existence of common large-scale climate drivers. Reductions in resilience are statistically linked to abrupt declines in forest primary productivity, occurring in response to slow drifting towards a critical resilience threshold. Approximately 23% of intact undisturbed forests, corresponding to 3.32 Pg C of gross primary productivity, have already reached a critical threshold and are experiencing a further degradation in resilience. Together, these signals reveal a widespread decline in the capacity of forests to withstand perturbation that should be accounted for in the design of land-based mitigation and adaptation plans.

Fig. 1. Temporal variations of forest resilience and its key drivers.

a, Spatial map of the temporal trend of TAC (δTAC). Positive δTAC values (for example, tropical forests) suggest a reduction in recovery rates and thus a decline in resilience, and vice versa for negative δTAC values (for example, boreal forests). The values are averaged over a 1° × 1° moving window for visual purposes. b, δTAC as in a binned as a function of climatological temperature and precipitation. The black dots indicate bins with average values that are statistically different from zero (two-sided Student’s t-test; P value ≤ 0.05). c, Frequency distribution of the differences in TAC computed for two independent temporal windows (2011–2020 minus 2000–2010) and shown separately for different climate regions. The numbers refer to the percentage of the observations lower and greater than zero; the asterisks indicate distributions with averages that are statistically different from zero (two-sided Student’s t-test; P value ≤ 0.05). The thin vertical line in each plot shows the distribution average. d, The cover fraction corresponding to each climate region and colour code reported in c and shown over the latitudinal gradient. e, The zonal average of the trend in TAC (δTAC) as determined by the three drivers (X) at 5° latitudinal resolution and the corresponding 95% confidence interval shown as a coloured line and shaded band, respectively. The colours reflect the three different driver categories: forest density, background climate and climate variability.

Fig. 2. Effect of forest management on forest resilience and interplay with GPP.

a, Frequency distributions of long-term TAC(2000–2020) for managed forests (MF) and intact forests (IF) located in a similar background climate. The coloured numbers report the respective averages, the top labels refer to the mean of the differences (diff.) in long-term TAC between managed and intact forests, and the asterisk indicates distributions that are statistically different (two-sided Student’s t-test; P value ≤ 0.05). b, The same as for a but for δTAC; the coloured numbers refer to the percentage of the observations lower and greater than zero (on the left and right of 0 on the x-axis, respectively). c, The same as for b but for the temporal correlation between annual GPP and TAC, denoted as ρ(GPP,TAC). d, A spatial map of ρ(GPP,TAC). e, ρ(GPP,TAC) binned as a function of climatological precipitation and temperature. The black dots indicate bins with average values that are statistically different from zero (two-sided Student’s t-test; P value ≤ 0.05). f, A spatial map of the areas, with different colours for the four combinations of positive/negative δGPP and δTAC. The cover fractions of each of the four classes for managed and intact forests are reported in stacked bars.

Fig. 3. Early-warning signals of ADs in intact forests.

a, Probability of occurrence of AD conditional on the values of δTAC for different severities of AD (expressed as anomaly n-times local standard deviation below the local mean, σ) shown separately for three different climate regions. The asterisks indicate probabilities statistically different from 0.5 (two-sided Student’s t-test; P value ≤ 0.05). b, TAC retrieved in the year preceding the occurrence of an AD (TAC_AD) binned as a function of climatological precipitation and temperature. c, Tolerance to TAC_AD (the absolute increase in TAC that an ecosystem in equilibrium can tolerate before reaching critical conditions) across a gradient of aridity index. The circle and whiskers refer to the average value and its 95% confidence interval; colours refer to the corresponding TAC_AD. Each binned aridity index ranging from 0 to 500 mm °C⁻¹ counts 10,868, 16,799, 728, 59 and 13 sampled pixels. d, Proximity to TAC_AD (proximity of present intact forests to their critical condition threshold) binned as a function of climatological precipitation and temperature. The black dots indicate bins with average values that are statistically different from zero (two-sided Student’s t-test; P value ≤ 0.05). Negative values of proximity to TAC_AD represent areas where the threshold resilience for AD (TAC_AD) has been already overpassed, and vice versa for positive values. e, Frequency distributions of proximity to TAC_AD shown separately for different climate regions and computed over the whole domain (blue) and over those areas experiencing a concomitant positive δTAC (red). The coloured numbers refer to the percentage of the frequency distribution lower and greater than zero (on the left and right of 0 on the x -axis, respectively) with respect to the whole domain.

Extended Data Fig. 1. Spatial variation of forest slowness.

(a) Spatial map of long-term $\mathit{TAC}$ computed for the whole 2000-2020 period. (b) Long-term $\mathit{TAC}$ binned as a function of climatological temperature and precipitation.

Extended Data Fig. 2. Performance and response functions of the resilience model.

(a) Observed versus modelled long-term $\mathit{TAC}$. Number of binned records (N), coefficient of determination (R²), mean squared error (MSE) and percent bias (PBIAS) are shown in labels, while the frequency distribution in color. (b) Predictors of $\mathit{TAC}$ and corresponding variable importance based on the random forest regression model of forest resilience. The four categories of environmental predictors are identified with hatched fill patterns; whereas the colors distinguish the different variables. (c) Dependence of $\mathit{TAC}$ on predictors of forest density. (d), (e) and (f) as (c) but for predictors of background climate, climate variability and autocorrelation, respectively.

Extended Data Fig. 3. Temporal trajectories of forest resilience.

Temporal changes in $\mathit{TAC}$ computed over a 3-year moving window and displayed with respect to the reference year 2002 separately for the global (a), tropical (b), arid (c), temperate (d) and boreal (e) regions. Continuous lines refer to the regional averages, whereas shaded areas show their 95% confidence interval magnified by a factor of 10 for visual purposes.

Extended Data Fig. 4. Sensitivity analysis of temporal changes in forest resilience (frequency distributions).

(a–e) Frequency distribution of the differences in $\mathit{TAC}$ computed for two independent temporal windows (2011-2020 and 2000-2010) shown separately for different climate regions and for the use of different quality flags of NDVI data (QF). Numbers refer to the percentage of the frequency distribution lower and greater than zero (on the left and right y-axis, respectively). (f–j), (k–o), (p–t), (u–y) and (z–ad) as (a–e) but computed for different gap filling analyses (GF), inclusion/exclusion of areas affected by abrupt declines (AD), percentages of missing data (PMD), percentages of forest cover (PFC) and spatial resolution (PSR), respectively.

Extended Data Fig. 5. Sensitivity analysis of temporal changes in forest resilience (climate spaces).

(a–b) Differences in $\mathit{TAC}$ computed for two independent temporal windows (2011-2020 and 2000-2010), separately shown for different quality flags (QF), binned as a function of climatological temperature and precipitation. Black dots indicate bins with average values that are statistically different from zero (two-sided Student’s t-test; P-value ≤ 0.05). (c–d), (e–f), (g–i), (j–l) and (m–o) as (a–b) but computed for different gap filling analyses (GF), inclusion/exclusion of areas affected by abrupt declines (AD), percentages of missing data (PMD), percentages of forest cover (PFC) and spatial resolution (PSR), respectively.

Extended Data Fig. 6. Effects of varying lagged temporal windows lengths.

(a–c) Trend in total $\mathit{TAC}$ binned as a function of climatological temperature and precipitation, separately shown for different temporal window lengths (TWL). Black dots indicate bins with average values that are statistically different from zero (two-sided Student’s t-test; P-value ≤ 0.05).

Extended Data Fig. 7. Temporal variations in environmental predictors.

(a-b) Temporal trends in environmental predictors of the category ‘forest density’ computed over a 3-year moving window and binned as a function of climatological precipitation and temperature. Black dots indicate bins with average values that are statistically different from zero (two-sided Student’s t-test; P-value ≤ 0.05). (c–f) and (g–j) as (a,b) but for environmental predictors of the categories ‘background climate’ and ‘climate variability’, respectively. Predictor acronyms are reported in Extended Data Table 1.

Extended Data Fig. 8. Climate and forest domains.

(a) Spatial map of climate regions. (b) Cover fraction of managed and intact forests for each climate region. (c) Spatial map of managed and intact forests. (d) Cover fraction of climate regions for each forest domain.