Emerging signals of declining forest resilience under climate change

Climate drivers

To explore the impact of climate on forest resilience (see the following sections), we used monthly averaged total precipitation, 2-m air temperature, evapotranspiration deficit and surface solar radiation downwards acquired from the ERA5-Land reanalysis product at 0.1° spatial resolution for the 2000–2020 period (https://cds.climate.copernicus.eu/cdsapp#!/home). Evapotranspiration deficit was quantified as the total precipitation minus evapotranspiration. In this study, we referred to climate regions as defined by the Köppen–Geiger world map of climate classification⁵¹ (http://koeppen-geiger.vu-wien.ac.at/present.htm). The original 31 climatic zones were merged into major zones and only those characterized by vegetation cover were included in our study (tropical, arid, temperate and boreal; Extended Data Fig. 8).

Vegetation dynamics

NDVI data acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument aboard the Terra satellite was used to derive changes in global vegetation for the period 2000–2020. We used cloud-free spatial composites provided at 16-day temporal resolution and 0.05° spatial resolution (MOD13C1 Version 6; https://lpdaac.usgs.gov/products/mod13c1v006/) and retained only pixels with good and marginal overall quality. The MODIS-derived NDVI dataset represents a state-of-the-art product of vegetation state whose retrieval algorithm is constantly improved⁵², and being derived from a unique platform and sensor, it is temporally and spatially consistent. Vegetation dynamics were analysed in terms of kNDVI, a nonlinear generalization of the NDVI based on ref. ²² and derived as follows:

kNDVI has recently been proposed as a strong proxy for ecosystem productivity that shows high correlations with both plot level measurements of primary productivity and satellite retrievals of sun-induced fluorescence²². In addition, kNDVI has been documented to be more closely related to primary productivity, to be resistant to saturation, bias and complex phenological cycles, and to show enhanced robustness to noise and stability across spatial and temporal scales compared to alternative products (for example, NDVI and near-infrared reflectance of vegetation). For these reasons, it has been retained in this study as the preferred metric to describe the state of the forest ecosystem.

To obtain an accurate estimate of resilience indicators, vegetation time series need to be stationary without seasonal periodic patterns or long-term trends⁵³. To this aim, vegetation anomalies were obtained from kNDVI data by first subtracting the multi-year 16-day sample mean and then removing linear trends from the resulting time series. Missing data, due for instance to snow cover affecting the retrieval of reflectance properties, have been gap-filled by climatological kNDVI values. The time series of kNDVI-based vegetation anomalies was used to derive resilience indicators and assess their spatial and temporal variations (see next sections).

Interannual changes in vegetation were assessed in terms of growing-season-averaged kNDVI. To this end, a climatological growing season that spanned months with at least 75% of days in the greenness phase was derived from the Vegetation Index and Phenology satellite-based product⁵⁴ (https://vip.arizona.edu/) and acquired for the 2000–2016 period at 0.05° spatial resolution. In addition, forest cover (FC) fraction was derived from the annual land-cover maps of the European Space Agency’s Climate Change Initiative (https://www.esa-landcover-cci.org/)⁵⁵ over the 2000–2018 period at 300-m spatial resolution. FC was retrieved by summing the fraction of broadleaved deciduous, broadleaved evergreen, needle leaf deciduous and needle leaf evergreen forest. FC was resampled to 0.05° to match the kNDVI spatial resolution.

Spatial patterns of slowness and its dependence on environmental factors

In this study, we quantified the resilience of forest ecosystems—their ability to recover from external perturbations—by the use of the 1-lag TAC (refs. ^3,4,5). Such an indicator was initially computed on the whole time series of vegetation anomalies (2000–2020) for forest pixels with less than 50% missing data in the original NDVI and FC greater than 0.05 and referred to in the text as long-term TAC. This analysis was used to assess the spatial patterns of the forest slowness mediated by environmental factors that affect plant growth rates and capacity to recover from perturbations. The long-term TAC was explored both in the geographic and climate space (Extended Data Fig. 1). In the climate space, long-term TAC was binned in a 50 × 50 grid as a function of average annual precipitation and temperature, both computed over the 2000–2020 period, using the average as an aggregation metric weighted by the areal extents of each record. We retained only bins with at least 50 records.

To explore the potential drivers of long-term TAC, we developed an RF regression model²³ and predicted the observed long-term TAC (response variable) based on a set of environmental features (predictors). The use of machine learning in general and of RF in particular, being nonparametric and nonlinear data-driven methods, does not require a priori assumptions about the functional form relating the key drivers and the response functions. The environmental variables include vegetation properties (FC and growing-season-averaged kNDVI) and climate variables (total precipitation, 2-m air temperature, evapotranspiration deficit and surface solar radiation downwards). Each of the climate variables was expressed in terms of average, coefficient of variation and 1-lag autocorrelation and resampled to 0.05° spatial resolution to match the spatial resolution of kNDVI. All environmental variables were computed annually and then averaged over time, except the autocorrelation that was computed directly for the whole period, analogously to the long-term TAC. This resulted in a set of 14 predictors representing the forest density, the background climate, the climate variability and its TAC in the observational period (Extended Data Table 1). The RF model was developed by splitting the observed long-term TAC into two separate samples: 60% of records were used for model calibration, and the remaining 40% were used to validate model performances in terms of coefficient of determination (R²), mean squared error and percentage bias (PBIAS). Each record refers to a 0.05° pixel. The RF implemented here uses 100 regression trees, whose depth and number of predictors to sample at each node were identified using Bayesian optimization. The general model formulation is as follows:

in which f is the RF regression model, X are the environmental predictors and ε_f are the residuals. We found that the model explains 87% of the spatial variance (R²) of the observed long-term TAC with a mean squared error of 0.007 and an average overestimation of 0.058 (PBIAS; Extended Data Fig. 2a). By definition, machine learning methods are not based on the mechanistic representation of the phenomena and therefore cannot provide direct information on the underlying processes influencing the system response to drivers. However, some model-agnostic methods can be applied to gain insights into the outputs of RF models. Here we used variable importance metrics to quantify and rank how individual environmental factors influence TAC (Extended Data Fig. 2b). Furthermore, using partial dependence plots derived from the machine learning algorithm RF, we explored the ecosystem response function (TAC) across gradients of vegetation and climate features (Supplementary Discussion 1 and Extended Data Fig. 2c–f).

CSD indicators

To explore the temporal variation in forest resilience, we used CSD indicators, here quantified in terms of temporal changes in TAC retrieved for two consecutive and independent periods ranging from 2000 to 2010 and from 2011 to 2020, and assessed the significance of the change in the sampled mean aggregated for different climate regions through a two-sided t-test (Fig. 1c). This analysis was complemented by the computation of TAC on the annual scale over a 2-year lagged temporal window (3-year window size) to track the temporal changes in CSD. This resulted in a time series of TAC with an annual time step.

We point out that temporal dynamics of annual TAC are driven by two processes: the changes in the resilience of the system that affect the velocity of the recovery from external perturbations and the confounding effects of the changes in autocorrelation of the climate drivers (X_ac) that directly affect the autocorrelation of NDVI. Given the specific goals of this study, we factored out the second process from the total TAC signal to avoid that an increasing autocorrelation in the drivers would affect our analysis and conclusions about the resilience and the potential increase in instability⁵⁶. For this purpose, we disentangled the temporal changes in TAC due to variations in autocorrelation in the climate drivers (({rm{TAC}}| {X}_{{rm{ac}}})) by adopting the space-for-time analogy and applied the RF model (f) at an annual time step (t) in a set of factorial simulations as follows:

The first term on the right side of equation (3) is the RF model simulation obtained by accounting for the dynamics of all predictors, and the second term is the RF model simulation generated by considering all predictors dynamic except the factors of autocorrelation in climate that are kept constant to their first-year value (year 2000). For such runs, we used predictors computed on an annual scale over a 2-year lagged temporal window, consistently to the TAC time series. We found that the direct effects of autocorrelation in climate have led to a positive trend of TAC in dry zones (due to the increasing autocorrelation of the drivers in these regions) and to an opposite effect in temperate humid forests (Supplementary Fig. 3). To remove these confounding effects, the estimated term ({{rm{TAC}}}^{t}| {X}_{{rm{ac}}}) is factored out from the TAC^t by subtraction to derive an enhanced estimate of annual resilience that is independent of autocorrelation in climate (Extended Data Fig. 3).

Long-term linear trends computed on the resulting enhanced TAC time series (δTAC) represent our reference CSD indicator used in this study to explore the changes in forest resilience. δTAC was quantified for each grid cell (Fig. 1a) and represented in the climate space following the methodology previously described (Fig. 1b). We then assessed the significance of the trends at bin level by applying a two-sided t-test for the sampled trend distributions within each bin. This significance test is independent from the structural temporal dependencies originating from the use of a 2-year lagged temporal window to compute the TAC time series.

Following an analogous approach described in equation (3), we disentangled the effect of the variation in forest density, background climate and climate variability on temporal changes in TAC (Fig. 1d,e). We recognize that other environmental factors not explicitly accounted for in our RF model could play a role in modulating the temporal variations in TAC. However, given the comprehensiveness of the suite of predictors used in equation (2) (Extended Data Table 1), it seems plausible that residuals mostly reflect the intrinsic forest resilience, the component intimately connected to the short-term responses of forests to perturbations, which is not directly related to climate variability. Forest ecosystem evolutionary processes could also play a role, but longer time series would be required to reliably capture these dynamics. Furthermore, abrupt declines (ADs) in the vegetation state and following recoveries, similarly to those potentially originating from forest disturbances (for example, wildfires and insect outbreaks), could influence the TAC changes. However, such occurrences, being distributed across the globe throughout the whole period, are expected to only marginally affect the resulting trend in TAC time series.

Sensitivity analysis

To assess the robustness of our results with respect to the modelling choices described above, we performed a series of sensitivity analyses for the difference in TAC retrieved for the two independent periods (2000–2010 and 2011–2020). To this aim, we tested their dependence on: the quality flag of the NDVI data used for the analyses (good, good and marginal); the gap-filling procedure tested on different periods (year and growing season); the inclusion or exclusion of forest areas affected by ADs; the threshold on the maximum percentage of missing NDVI data allowed at the pixel level (20%, 50% and 80%); the threshold on the minimum percentage of FC allowed at the pixel level (5%, 50% and 90%); and the pixel spatial resolution used for the analyses (0.05°, 0.25° and 1°). In addition, we tested the sensitivity of the trend in total TAC signal on the moving temporal window length used to calculate autocorrelation at lag 1. Results obtained for the different configurations were compared in terms of frequency distributions, separately for climate regions (Extended Data Fig. 4), and further explored in the climate space (Extended Data Figs. 5 and 6). Outcomes of the sensitivity analysis are discussed in Supplementary Discussion 2.

Interplay between GPP and CSD

Resilience and GPP interact with each other through mutual causal links. On one hand, a reduction in forest resilience makes the system more sensitive to perturbations with potential consequent losses in GPP (ref. ²⁶). On the other hand, a reduction in GPP may lead to a decline in resilience according to the carbon starvation hypothesis, and may be associated with increasing hydraulic failure⁴⁶. To explore the link between forest resilience and primary productivity, we quantified the correlation between TAC and GPP. Estimates of GPP were derived from the FluxCom Model Tree Ensemble for the 2001–2019 period at 8-daily temporal resolution and 0.0833° spatial resolution and generated using ecosystem GPP fluxes from the FLUXNET network and MODIS remote sensing data as predictor variables³⁶ (http://www.fluxcom.org/). Annual maps of GPP were quantified and resampled to 0.05° to match the temporal and spatial resolution of TAC time series. The Spearman rank correlation (ρ) was then computed between annual GPP and TAC over a 1° spatial moving window to better sample the empirical distribution of the two variables (Fig. 2d). The significance of ρ(GPP,TAC) was assessed over the climate space separately for each bin (Fig. 2e), similarly to the approach used to test the significance of δTAC. Furthermore, we explored the relationships between the trend in GPP (δGPP) and the trend in TAC (δTAC) by clustering the globe according to the directions of the long-term trajectories of the above-mentioned variables (Fig. 2f).

Disentangling the impact of forest management

To characterize TAC on different forest types and disentangle the potential effects originating from forest management, results were separately analysed for intact forests and managed forests. Intact forests were considered those forest pixels constituting the Intact Forest Landscapes⁵⁷ dataset (https://intactforests.org/). Intact Forest Landscapes identifies the forest extents with no sign of significant human activity over the period 2000–2016 based on Landsat time series. The remaining forests pixels—not labelled as intact—were considered as managed forests (Extended Data Fig. 8). The resulting forest type map is consistent with those used for United Nations Framework Convention on Climate Change reporting⁵⁸, although with more conservative estimates of intact forests in the boreal zone due to the masking based on FC and percentage of missing data applied in this study.

We analysed the differences in long-term TAC (computed for the whole 2000–2020 period) between managed and intact forests by masking out the potential effect of climate background. To this aim, we compared the climate spaces generated separately for managed and intact forests by extracting only those bins that are covered by both forest classes. The resulting distributions—one for each forest class—have the same sample size, and each pair of elements shares the same climate background. Potential confounding environmental effects on average recovery rates are, therefore, minimized. We then applied a two-sided t-test for analysing the significance of the difference in the sampled means (Fig. 2a). An analogous approach was used to test the differences in δTAC and ρ(GPP,TAC) between managed and intact forests (Fig. 2b,c).

Early-warning signals of abrupt forest declines

When forest ecosystems are subject to an extended and progressive degradation, the loss of resilience can lead to an AD (refs. ^3,4,5). Such abrupt changes can trigger a regime shift (tipping point) depending on the capacity of the system to recover from the perturbations (Supplementary Methods 1 and 2). We investigated the potential of changes in TAC as early-warning signals of ADs in intact forests over the 2010–2020 period. To this aim, we quantified at the pixel level ADs as the events occurring on a certain year when the corresponding growing-season average kNDVI was more than n-times local standard deviation below the local mean. Local mean and standard deviation (σ) were computed over the 10-year antecedent temporal window (undisturbed) period and n varies between 1 and 6 with higher values reflecting more severe changes in the state of the system. For each pixel and for each fixed n value, we recorded only the first AD occurrence, thus imposing a univocal record for each abrupt change in the state of the system.

We then explored whether the retrieved ADs were statistically associated with antecedent high values of δTAC. To avoid confusion with the attribution of causality, for each AD that occurred at time t (over the 2010–2020 period), we derived the δTAC over the temporal window 2000 − (t − 1). The resulting trend in δTAC is therefore antecedent and independent of the changes in vegetation associated with the AD. Then, for each pixel with an AD at time t, we also extracted randomly one of the undisturbed (with no AD) adjacent pixels and retrieved δTAC over the same temporal window. This analysis produced two distributions of δTAC associated with pixels with and without ADs (AD and no AD, respectively). The two distributions have the same size and each pair of elements shares similar background climate. We calculated the probability of occurrence of AD conditional on the trend in δTAC (({rm{AD}}| delta {rm{TAC}})) as the frequency of ADs for which (delta {rm{TAC}}left(mathrm{AD}right)| > delta {rm{TAC}}left(mathrm{no; AD}right)), and the significance of the difference in the two sampled means (AD and no AD) was evaluated through a two-sided t-test. Probability and significance were assessed for different climate regions and severity of ADs (Fig. 3a). High statistically significant probabilities suggest that the AD is following the drifting towards a critical resilience threshold plausibly associated with changes in environmental drivers.

We complemented the aforementioned analyses by retrieving the tolerance and proximity to AD, hereafter determined for a 3σ severity. We first quantified the TAC that proceeded the occurrence of an AD and followed a progressive loss of resilience as captured by positive δTAC. This value, hereafter referred to as abrupt decline temporal autocorrelation (TAC_AD), reflects the TAC threshold over which we observed an abrupt change in the forest state (Fig. 3b). The tolerance to AD was quantified as the difference between the local TAC_AD and the TAC value averaged over the 2000–2009 period to characterize the pre-disturbance conditions. The tolerance metric was explored across a gradient of aridity index⁵⁹ (Fig. 3c).

TAC_AD can be directly retrieved only on those forest pixels that have already experienced an AD. As a considerable fraction of undisturbed forests could potentially be close to their critical TAC threshold, or even have already passed it, it is important to determine their TAC_AD. To this aim, we developed an RF regression model that expresses the TAC_AD as a function of the set X of environmental variables used in model f (equation (2)) but excluding the autocorrelation in climate drivers (X_reduced) already disentangled in the TAC signal. The general formulation is as follows:

in which g is the RF regression model, X_reduced are the environmental predictors and ε_g are the residuals. Implementation, calibration and validation of g follow the same rationale described before for the f model. We found that the RF model explains 50% of the variance (R2) of the observed TAC_AD, with a mean squared error of 0.019 and an average underestimation of 0.86 (PBIAS).

The RF model was then used to predict the TAC_AD over the whole domain of intact forests and served as input to quantify the proximity to AD of undisturbed forest pixels at the end of the observational period (year 2020). Here we defined the proximity metric as the difference between the value of TAC in 2020 and TAC_AD. Proximity takes negative or zero values when TAC_AD has already been reached (({{{rm{TAC}}}^{2020}ge {rm{TAC}}}_{{rm{AD}}})) and positive values when there are still margins before reaching the critical threshold (({{{rm{TAC}}}^{2020} < {rm{TAC}}}_{{rm{AD}}})). Together (delta {rm{TAC}} > 0) and ({{{rm{TAC}}}^{2020}ge {rm{TAC}}}_{{rm{AD}}}) therefore represent the most critical conditions, as they indicate that the critical resilience threshold for AD has already been reached and the ecosystem is continuing to lose its capacity to respond to external perturbations. We finally quantified the amount of GPP potentially exposed to such critical conditions by linearly extrapolating the GPP for the year 2020 (available GPP data stop in 2019) and overlaying it on the map of critical conditions (proximity to ({rm{AD}} < 0) and (delta {rm{TAC}} > 0)).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

Source: Ecology - nature.com