Enhanced structural diversity increases European forest resilience and potentially compensates for climate-driven declines

Abstract

Ecosystem resilience represents the capacity to withstand and recover from perturbations, an important property for ecosystem function in an era of escalating disturbances. Whilst studies relate large-scale forest resilience to abiotic factors, knowledge gaps remain regarding the link to aspects of functional diversity, such as forest structural diversity (FSD), a factor controllable through management. Here we use spaceborne lidar-derived FSD to quantify its role in determining resilience in Southern and Central European forests. A random forest isolates the FSD-resilience interplay by disentangling confounding environmental factors. Results show a positive FSD-resilience relationship in 80% of forests. Canopy complexity is a better predictor of resilience than spatial variability in canopy height. This emergent relationship is explored as an adaptation measure to preserve resilience under warming scenarios, with the decline associated with 1 ^∘C of warming compensated by a 10% increase in canopy complexity. The findings suggest that management promoting canopy complexity potentially compensates warming-driven declines in resilience.

Forest structural diversity is positively associated with forest resilience in about 80% of Southern and Central European forests, implying that increasing canopy complexity may offset warming-driven declines in resilience, based on an analysis of remote sensing data and random forest modelling.

Introduction

Forests are a fundamental element of the global carbon cycle, provide a protective local effect against climate extremes, and are biodiversity hotspots^{1, 2–3}. Moreover, forests are a key element to reach the climate neutrality, Forest Strategy and Nature Restoration goals set by the European Union, whilst the Biodiversity Strategy for 2030 aims to build resilience to future threats^{4, 5, 6–7}. However, recent evidence shows a decline of forest resilience under increasing global change and disturbances, with rising tree mortality in Europe, particularly from drought-driven bark beetle infestations^{8, 9, 10, 11, 12, 13–14}.

Ecosystems constantly undergo perturbations from their equilibria¹⁵. Their capacity to withstand and recover from such perturbations has been characterised by resilience metrics, relating to the speed of recovery or time period between the disruption and restoration of full or partial system function^16,17. Recently, critical transition theory has developed early warning signals to detect the bifurcation points of natural systems, where the equilibrium-restoring forces of the system to natural perturbations weaken^{18, 19, 20–21}. A critical slowing down is observed through rising 1-lag AutoCorrelation (AC1) (the correlation between consecutive points in a time series) and variance in the system response to small perturbations, indicators that function as proxies for the recovery rate from larger perturbations and disturbances^19,21,22.

Many large-scale ecosystem studies use AC1- and variance-based resilience metrics calculated on time series of remotely-sensed vegetation indices that include NDVI (Normalized Difference Vegetation Index) and kNDVI (kernel-NDVI)^{23, 24, 25, 26–27}. In forests, these metrics have demonstrated utility as early warning signals of mortality, providing insights into forest recovery from specific perturbations and revealing declining resilience across many regions^{8, 9–10,28, 29–30}.

Increasing biodiversity has long been suggested as a forest management technique to build forest stability and resilience^{31, 32–33}. Forest Structural Diversity (FSD) is a measure of variability within the canopy structure relating to aspects such as species richness, functional diversity and age distribution, and, importantly, may be shaped by management^34,35. FSD may be described horizontally, reflecting spatial variability, vertically, reflecting diversity in density among canopy layers, or as a combination of the two. At local scales, there are studies relating some specific aspects of forest functional diversity (such as diversity of hydraulic traits) to various measurements of vegetation resistance and recovery time from specific disturbances (such as droughts). However, there remain knowledge gaps in the relationship and the role of diversity in recovery, particularly at large scale, and the effectiveness of management techniques^36,37.

New spaceborne LIDAR (LIght Detection And Ranging) sensors, such as NASA’s Global Ecosystem Dynamics Investigation (GEDI) mission, enable the monitoring of high-resolution quasi-global FSD, offering unprecedented opportunities to assess the resilience-diversity relationship at large-scale^38,39. Such Earth observation technologies offer unprecedented opportunities to assess the resilience diversity relationship at a large scale. This is particularly relevant to help forest managers, stakeholders and policymakers, in identifying suitable adaptation actions to maintain or even improve diversity and resilience in the face of climate change.

In this paper, we assess the relationship between resilience and FSD in Central and Southern European forests. AC1- and variance-derived resilience metrics are developed from a 0.05^∘ and 8-day kNDVI time series, whilst the work utilises three separate GEDI-derived FSD metrics relating to horizontal (FSDH), vertical (FSDV) and combined (FSDH+V) canopy vegetation distribution. A random forest model isolates these FSD-resilience relationships from confounding, spatially distinct, environmental factors, enabling a deeper understanding of the local and Europe-wide relationships and identifying which forest management techniques most effectively promote resilience. Finally, the FSD-resilience relationship is examined as a function of temperature to investigate whether forest management may compensate for warming-induced declines in resilience.

This study uses large-scale remote sensing data and reanalysis climate products to directly link recent developments in forest diversity mapping and resilience. The findings may support optimal land management strategies at a time of increasing forest vulnerability.

Results

The forest resilience and structural diversity metrics

Temporal autocorrelation (AC1) derived and variance-derived restoration rate magnitudes are calculated from a 19-year kNDVI time series of Central and Southern European forests (Fig. 1a and Supplementary Information A, respectively). These metrics, established from critical transition theory, refer to the absolute magnitude of the forest recovery rate to perturbations, whereby a larger value signifies a faster recovery (Methods)^{16,18,19,21,25}. The AC1-derived metric relates to the ecosystem memory via the 1-lag autocorrelation of perturbations in the kNDVI time series, whilst the variance-derived metric relates to the ecosystem stability^21,40.

Fig. 1 Maps of forest restoration rate and structural diversity. [Images not available. See PDF.]

a The spatial distribution of the AC1-derived absolute restoration rate magnitude of forests. b The spatial distribution of the combined horizontal and vertical forest structural diversity (FSDH+V) metric. The other metrics of forest restoration rate and forest structural diversity are available in Supplementary Information A.

Whilst the figure shows variability in forest recovery speeds across Europe, the 0.05^∘ spatial resolution pixels are not directly comparable to each other in terms of ecosystem resilience, mainly due to spatial differences in the background climate. For example, arid zones of Mediterranean Europe are more prone to a higher temporal autocorrelation and variability in the time series of the vegetation indices compared to Central Europe, due to the importance and variability of precipitation on vegetation growth (Supplementary Information G)^8,23,29,30. Only once the confounding effects of different local conditions are accounted for can the recovery rates be interpreted as resilience.

GEDI data is used to derive three metrics of structural diversity in the forest canopy³⁹. These metrics, described in the Methods and summarised in Table 1 are: FSDH, relating to the horizontal diversity and measured as the spatial variability in canopy height; FSDV, relating to the vertical diversity and measured as the inverse excess kurtosis in the vertical profile, with higher values corresponding to a more even distribution of vegetation in the profile; and FSDH+V, a combined horizontal and vertical diversity, measured via the Shannon entropy of relative heights of the vertical profile and canopy cover, and therefore a measure of variability and complexity in canopy substructure. The GEDI observations are constrained to the latitudinal limit  ≲50^∘, and therefore Northern Europe is not included in the analysis.

Table 1. Summary of the three forest structural diversity metrics considered in the analysis

Diversity metric	Calculation	Description
Horizontal diversity (FSDH)	Standard deviation in the top (98th percentile) of the GEDI vertical profile across the 0.05∘ pixel	The spatial variability in the forest canopy height within each pixel
Vertical diversity (FSDV)	Inverse excess kurtosis in the GEDI vertical profile, relative to a Gaussian distribution, averaged across the 0.05∘ pixel	A measure of the excess vegetation distributed in the upper and lower canopy and therefore related to the evenness of the vertical vegetation distribution
Combined horizontal and vertical diversity (FSDH+V)	The Shannon Entropy of the relative forest heights and canopy cover over the 0.05∘ pixel	A quantification of the distribution of vegetation both spatially and in terms of canopy layers and therefore a proxy for variability and complexity in canopy substructure

The metrics are defined in such a way that higher values correspond to higher levels of structural diversity.

Higher levels of FSD are found in central European forests, in particular in Germany, whilst lower levels are found in southern, Mediterranean regions (FSDH+V, Fig. 1b; FSDH & FSDV, Supplementary Information A). However, there is variability between the FSD metrics and correlation between them is low (Supplementary Information I). For example, we observe that areas with high horizontal structural diversity do not necessarily have high vertical diversity, and vice-versa, potentially driven by the effect of precipitation on canopy complexity⁴¹. This indicates that the different metrics relate to different aspects of structural diversity.

The random forest models for resilience estimation

A non-linear, multi-feature Random Forest (RF) model (described in Methods) is used to separate and isolate the effects of FSD on the restoration rate, accounting for confounding environmental factors and differing local perturbations, and enabling us to interpret the relationship in terms of resilience. These factors include the mean, variability (C.V.) and temporal autocorrelation (AC1) in the temperature (2m temp), precipitation (Precip.), Surface Solar Radiation (SSR) and Vapour Pressure Deficit (VPD), as well as soil quality (Soil Carbon 30cm), mean kNDVI, topography (Elevation S.D.), and forest cover (Supplementary Information H).

A series of 6 RF models are trained to predict the two (AC1- and variance-driven) restoration rates from these factors, with each model including one of the three FSD metrics. The models have a strong performance in predicting the test data as observed through the high coefficient of determination (R² ≃ 0.7 for AC1-derived restoration rates), low mean squared error (MSE  ≃ 0.095) and low percentage bias (PBIAS  ≃ 7.2%) between observed and modelled values (Fig. 2a, Supplementary Information B).

Fig. 2 The Random forest model of forest resilience. [Images not available. See PDF.]

a Performance metrics of the Random Forest model measured through comparison of the testing dataset modelled and observed AC1-derived restoration rate for the combined structural diversity metric (FSDH+V). b Predictor importance of that same model measured through the increase in model mean square error (MSE) expected from the removal of each model predictor, with 95% confidence intervals provided via a bootstrapping process. Models for other metrics of forest restoration rate and forest structural diversity are available in Supplementary Information B.

The importance of model predictors can be measured through the expected average increase in the MSE resulting from removing each feature from the model (Methods). Due to sampling variability and the stochastic nature of RF models, some variability in the order of predictors between models may be expected, and the variability between models is calculated via 95% confidence intervals. The FSDH+V metric, representing diversity and complexity within the canopy substructure, is observed as the second most important predictor of the resilience, after forest cover percentage (Fig. 2b). FSDV, relating to the evenness of vegetation distribution in the vertical profile, is slightly less important than FSDH+V, however FSDH, relating to variability in canopy height, is one of the least important features in predicting the resilience (Supplementary Information B). These results suggest that different FSD metrics show large differences in their capacity to explain the variability in forest resilience. As such, understanding the differences in the underlying forest structure represented by the different metrics is required to properly characterise the relationship between forest resilience and forest structural diversity, and ultimately identify and design the most effective management strategy to enhance forest resilience.

The emergent relationship between European forest resilience and structural diversity

To elucidate the emerging relationship between forest resilience and structural diversity in Central and Southern Europe, the partial dependence of the restoration rate on the different FSD metrics is explored. Partial dependence figures disentangle the dependence of the predicted variable, the restoration rate magnitude, from the FSD metrics, whilst accounting for confounding variables (Methods). Taking into account the different environmental conditions at each pixel enables the derivation of the functional relationship between the forest recovery rate and the FSD metrics and therefore enables the relationship to be interpreted in terms of resilience.

As each of the three FSD metrics increases, the partial dependence of the AC1-derived and variance-derived resilience increases (Fig. 3 and Supplementary Information C, respectively). To enable direct comparison, the FSD metrics are rescaled to display 95% of their variability across the studied European forests. The greatest increases in resilience relative to changes in FSD occur around the 50^th to 75^th percentiles in the variability across Europe, with diminishing returns above that, and less dependency in the relationship below. For each FSD metric, the FSD-resilience relationship is significant compared to a distribution of randomly permuted FSD metric values (p < 0.01).

Fig. 3 The relationship between European forest resilience and structural diversity. [Images not available. See PDF.]

The partial dependence of the AC1-derived restoration rate of selected European forests, as a function of three metrics of forest structural diversity: horizontal (FSDH), vertical (FSDV) and combined horizontal + vertical (FSDH+V). The diversity metrics are scaled to display 95% of the variability observed across Europe. A 95% confidence-level band of model and FSD metric uncertainty is displayed via a bootstrapping process. The equivalent figure for the variance-derived resilience metric can be found in Supplementary Information C.

Across the variability observed in the diversity metrics, the predicted AC1-derived restoration rate range is 5 times larger for the combined horizontal and vertical (FSDH+V) metric than for the horizontal metric (FSDH), and almost twice as large compared to the vertical metric (FSDV). The steepest slopes in the partial relationships are also found in the order of FSDH+V, followed by FSDV, and finally FSDH. The greatest resilience-FSD dependency is therefore observed in FSDH+V, suggesting that resilience has the strongest dependency on complexity within the canopy substructure. FSDV shows the next largest dependency, suggesting the greater distribution of vegetation at the upper and lower sections of the canopy vertical profile is also relatively important. Meanwhile, FSDH, which is connected to spatial variability in the canopy height, shows a much lower impact on forest resilience, suggesting that the structural diversity within individual stands is much more important for forest resilience than the structural variability between stands.

Unmanaged forests as well as stable (unharvested during the data-collection period) managed forests are included in the data, and the observed relationship is driven by both the natural variability in forest structure as well as the impact of forest management on structural features. Horizontal variability in canopy height (but low structural complexity) across a 0.05^∘ spatial resolution region is associated with intensive management practices such as clear-cutting⁴². The results therefore suggest that clear-cutting methods of forest harvesting and monospecific, even-aged stands are less effective in enhancing resilience than management practices enhancing canopy complexity or increasing vegetation distribution in the vertical profile, such as selective logging, or other silvicultural systems leading to uneven-aged and uneven-size stands.

The local relationship between forest resilience and structural diversity

The relationship between resilience and FSD can be isolated from the confounding variables to present the expected local relationship given the specific environmental conditions at each pixel. The local partial dependence is produced using an Individual Conditional Expectation (ICE) method (Methods), and the pixelwise derivative is extracted as the magnitude of the local FSD-resilience dependence. A positive value signifies that as the FSD increases, the restoration rate (and therefore resilience) increases.

There is a positive relationship between AC1-derived resilience and the combined FSDH+V metric in pixels representing over 80% of the forest area considered in the analysis (Fig. 4). FSDV and FSDH show a slight reduction in the proportion of forest area exhibiting a positive relationship, with 73% and 70% respectively, and a similar spatial distribution (Supplementary Information D). The strength of the relationship, as measured by the size of the derivative, is, however, significantly smaller in magnitude for FSDH.

Fig. 4 The local relationship between forest resilience and structural diversity. [Images not available. See PDF.]

The local-level partial dependence between the AC1-derived restoration rate and the combined horizontal and vertical diversity (FSDH+V) metric. The relationship direction and magnitude are determined at each pixel by the derivative of the individual conditional expectation curves between the resilience and diversity metric. The percentage of forest area showing a positive and negative relationship is also displayed. The figures for the alternative resilience and diversity metrics can be found in Supplementary Information D.

The forest areas that show the strongest positive relationship tend to be in central Europe, where the highest FSD is found, however there are large areas exhibiting a positive relationship throughout the map. The majority of pixels exhibit a relatively small positive relationship, reflecting the fact that whilst the relationship is net positive, there is some uncertainty in the application of a globally derived model to specific local circumstances. Indeed, some limited areas of forest show a slight inverse relationship, for example, in Northern France. However, the relationship is relatively weak and may simply emerge as a statistical artefact, as the pattern is not replicated in the variance-derived resilience metric or in the FSDV models (Supplementary Information D).

The interplay between forest resilience and structural diversity under climate warming

As Europe is a climate change hotspot with a warming rate about twice the global mean, it is essential to explore the link between forest resilience and diversity within a context of rising temperatures⁴³. Therefore the partial dependence of the AC1-derived restoration rate magnitude on the FSDH+V metric is expressed as a function of the mean daily 2m temperature (Methods).

Southern and Central European forests with faster recovery rates are more structurally diverse and tend to be located in areas with lower background temperatures (Fig. 5a and Supplementary Information F). This link between resilience and temperature is consistent with previous studies, and likely results from the increased risk of tree mortality associated with extreme events such as droughts and heatwaves^{8,23,27,29,30}. Of particular interest are the central regions of the figures, where the most data lies and therefore where the model behaviour is averaged from a larger sample (and less influenced by small patches of pixels that may contain outlier values). White isolines show the balance between temperature and FSD that maintains a constant predicted resilience, and suggests that—within the model—declines in forest resilience associated with increasing temperatures may be compensated by a corresponding increase in structural diversity.

Fig. 5 The interplay between forest resilience, structural diversity and temperature. [Images not available. See PDF.]

a Two-dimensional partial dependence plot of the AC1-derived restoration rate of Central and Southern European forests as a function of the mean background temperature and the combined horizontal and vertical structural diversity (FSDH+V) metric. b Distribution of the separate biogeographical regions considered in the analysis. c Partial dependence of the AC1-derived restoration rate on the FSDH+V metric for the three different biogeographical regions of European forests. Two different mean background temperature scenarios are shown, the 2003-2021 mean (line) and a 1 ^∘C increase on this (dashed). A bootstrapping process is used to derive 95% confidence-level bands displaying the model and FSD metric uncertainty. The figures for the alternative resilience and diversity metrics can be found in Supplementary Information F.

The partial dependence of the AC1-derived restoration rate on forest structural diversity, shown in Fig. 3, can be reproduced with lines separately representing three biogeographical groupings of European forests: alpine, temperate (including continental areas) and Mediterranean (Fig. 5b, c). An incremental change of 1 ^∘C is also added to each pixel’s mean daily temperature to simulate the short-term temperature rises associated with global warming, and the resulting resilience-FSD dependence shown. In order to maintain simplicity, non-temperature changes associated with global warming, such as increased variability of precipitation, are not considered.

Figure 5c can be interpreted in terms of the change (required by the model) in the combined FSDH+V metric in order to compensate increases in the regional temperature and maintain the current observed levels of resilience (assuming other environmental factors stay constant). In each region, as average temperatures rise, resilience is maintained when FSDH+V increases beyond the 95% confidence interval. In the near-term, in order to offset a warming of 1 ^∘C the metric would need to shift from the current FSDH+V median regional value to values 0.054 (0.018–0.107), 0.082 (0.041–0.123), and 0.120 (0.040–0.200) higher for the respective alpine, temperate and Mediterranean regions (95% confidence interval in brackets). Expressed in terms of the changes in the diversity metric relative to the distribution in each biogeographical region, offsetting 1 ^∘C of warming requires respective 4.8% (1.5–9.4%), 7.1% (3.6–10.6%) and 10.2% (3.2–17.4%) increases in the median combined FSDH+V diversity metric, relative to the current observed alpine, temperate and Mediterranean regional median diversity. For the variance-derived restoration rate (Supplementary Information F) the corresponding relative increases are, respectively, 4.8% (1.5–9.4%), 5.3% (0.1–8.8%) and 17.4% (8.3–23.4%) It should be noted that the scaling in all models is non-linear, and the exercise pushes the models towards areas of sparse training data where it is less stable. However, for the purposes of projecting climate warming in the short term, the data coverage is reliable.

Discussion

This paper assesses the relationship between resilience and forest structural diversity in selected European forests. Three different remotely-sensed metrics representing different aspects of forest structural diversity are investigated: the horizontal variability in canopy height, FSDH, the evenness of vegetation distribution in the vertical profile, FSDV, and the complexity and variability in the canopy substructure, FSDH+V. By alternately analysing three distinct FSD metrics we show the Europe-wide relationship between forest resilience and structural diversity, the relative importance of each predictor for resilience, and understand the implications for forest management.

For all three FSD metrics, in over 70% of the selected forests, there is a positive relationship between resilience and FSD, demonstrating that forests with higher structural diversity are associated with a faster recovery from perturbations. However, the findings show that forest structural complexity and diversity in the vertical canopy profile, are much stronger predictors of resilience than spatial variability in the canopy height, as seen from the model feature importance (Fig. 2), the partial dependence relationship (Fig. 3), and the local partial dependence maps (Fig. 4). Canopy complexity is particularly important, associated with a positive resilience relationship in 81% of Central and Southern European forests, and predicts a variability in the recovery rate 5x larger than spatial variability in canopy height. The results appear to be robust, given the use of multiple FSD indicators and the fact that both AC1-derived and variance-derived resilience indicators show agreement, in line with critical transition theory ^20,21. The analysis is also conducted at two different spatial resolutions, with the overall results remaining consistent at both 0.05^∘ and 0. 1^∘.

Significantly during an era of global warming, the work shows that the strength of the resilience-FSD relationship varies as a function of temperature. Increasing temperatures and more frequent droughts are expected to impact forest ecosystems, which already exhibit evidence of declining resilience^8,10. This work indicates that it may be possible to compensate for some of these resilience declines by increases in FSD (Fig. 5). A 10% relative increase in the combined FSDH+V diversity corresponds to an offset of at least 1 ^∘C of warming in each selected major European biogeographical region. Therefore, reversing the decline in biodiversity and promoting forest management that ensures canopy complexity may help mitigate the impacts of climate change on forest ecosystems.

It is important to acknowledge that whilst RF models predict forest recovery rates from the input variables, they do not necessarily demonstrate causation. Therefore, whilst it is clear there is a strong link between forest resilience and FSD, and a theoretical basis for diversity promoting ecosystem resilience, further analysis is needed to demonstrate a causational effect. Additionally, the RF model displays averaged conditions, which may not always hold at local levels, for example, across uncommon or extreme climate domains or specialised forest types, and where there is limited training data. It is also important to note that whilst temporal AC1- and variance-based indicators have demonstrated utility as early warning signals in bifurcation-induced tipping points—whereby gradual changes push the system to a threshold point between two regimes—there is some debate over their utility for other tipping scenarios^18,22,40. Noise-induced tipping scenarios—driven by stochastic forcing and without accompanying change in drivers—are, by definition, less predictable and therefore less likely to exhibit critical slowing down, though rising variance may still be an useful predictor of transition⁴⁴. In rate-induced tipping—whereby the pace of change is too rapid for the system to track its stable state—some studies support the utility of traditional early-warning signals, whilst others suggest that alternative approaches, e.g. via deep learning, may be more appropriate^45,46.

There are several reasons for the mechanism through which greater forest biodiversity reduces the size and duration of disturbances. It is established that forest fires are less likely in heterogeneous landscapes, featuring mixed age and species forests, with discontinuous flammability^47,48. Similarly, uneven-aged, mixed species stands with complex canopy structure are more resistant to windthrows^49,50. In these disturbance regimes, diversity generally acts to disrupt the propagation of the disturbance. However, diversity can also function in a way that reduces the overall density and quantity of trees vulnerable to specific disturbances. For example, forests containing a mixture of tree species, tree ages and sizes are less vulnerable to pests due to a greater dispersion of vulnerable hosts, with the difficulty in locating potential targets acting to both reduce and slow the attack⁵¹. Indeed, a recent study finds that enhancing ecosystem heterogeneity in European forests could reduce biomass loss due to climate-driven disturbances—such as wildfires, strong winds and insect outbreaks—by approximately 18%¹⁴.

The functional redundancy of diverse pre-disturbance ecological traits as well as the complementary use of resources and partitioning of the niche between species, may enable a portfolio effect of varying responses to disturbances^{52, 53–54}. Variability in hydraulic traits buffers variation in ecosystem transpiration during dry periods, stabilising forest productivity in the face of drought^{55, 56, 57, 58–59}. Meanwhile, simulated plot level studies of post-disturbance tropical forest growth show that increased functional diversity reduces recovery time⁶⁰. After disturbances, productive, pioneer species rapidly occupy canopy gaps and create a shaded environment for successional species with more conservative light-response, and this ‘response diversity’ drives the overall forest recovery^61,62. Finally, for unusual or novel perturbations, diversity increases the chance that one species or configuration will survive, the so-called ‘insurance hypothesis’¹⁵.

This study builds on these more specific plot-level studies, indicating that resilience, in general, directly relates to the structural diversity that often acts as a proxy for various other forms of vegetation diversity. Whilst the majority of existing studies measure recovery times and resistance to specific, localised disturbances or plot re-growth (e.g. via dendrochronology), by using remotely-sensed data, significantly more forest area can be included in the analysis, allowing for a generalisation of the relationship for large-scale emergent behaviours and the examination of spatial patterns. Such large-scale studies are particularly useful from a common forest policy perspective, such as the European Green Deal and the biodiversity strategy for 2030^5,6. The work also demonstrates the utility of resilience, as exhibited via novel AC1 and variance timeseries metrics, in detecting forest vulnerability to disturbances at a global scale, and the analysis provides a framework for future studies of forest tipping points and recovery.

Our study finds that increased canopy complexity is particularly important in building resilience, representing a portfolio of functional traits with differing disturbance responses and recovery times. Structural complexity is an indicator that reflects different aspects, including variability in tree species, functional traits, and tree height, and potentially differences in rooting depths that increase the resilience of the forest by allowing access to resources such as water during droughts⁶³. It is therefore important to understand the different aspects of diversity that the GEDI FSD metrics represent^64,65. Excess vegetation at the extremes of the vertical profile acts mostly as a proxy for alpha-diversity. Alpha diversity represents average species richness in a small-scale area, though in this case the diversity metric is more closely aligned with functional diversity, with the more even vegetation vertical profile reflecting diversity in ecological roles and traits. Spatial variability in the canopy height acts as a proxy for beta-diversity, indicating differences in species composition and traits between habitats. The combined (FSDH+V) canopy complexity metric, meanwhile, acts as a combination of horizontal and vertical components of forest structure, measuring the spatial variability and complexity in the vertical profile. At the scale of  ∽5km, this analysis suggests that beta-diversity and the combined horizontal and vertical diversity metric are particularly important to the resilience of landscapes. The broader microclimate exposure of multiple leaf layers and more diverse rooting of the vertical diversity reduces cavitation and drought risk, and increases survival via the insurance hypothesis, whilst the horizontal diversity likely acts to reduce the propagation of disturbances.

Overall, the study underlines the importance of preserving and enhancing structural diversity within forest management practices. It is suggested that changes in European forest extent, structure and composition through certain management methods have created conditions more vulnerable to disturbances, with the magnitude of damage similar to that resulting from climate change⁴⁹. Such disturbances can push forests towards tipping points, and there is evidence to suggest that less intensively managed forests tend to be more resistant to disturbances⁶⁶. Our study shows that, at large scales, practices that promote structural complexity, such as selective logging, mixed-species forestry and uneven-aged management, are likely to be more effective in maintaining or even increasing forest resilience than practices that homogenise vertical forest structure, such as clear-cutting and monospecific stands. Whilst resistance to fire disturbances—not considered in our analysis—may require alternative strategies, such as reducing stand density and fuel continuity, the work aligns with many local studies showing that certain proactive management techniques such as assisted migration and increasing species diversity and genetic exchange, can build resistance to specific local disturbances relating to climate, drought and pests^67,68. Indeed, a recent review of natural disturbance forest management strategies suggests that whilst situation-specific evaluation of the effectiveness of adaptation is important, a new research avenue of large-scale experiments is required to address ‘tipping point management’ to avoid forest resilience thresholds⁶⁸. Therefore, in the context of increasing forest disturbances and anthropogenic pressures, our study provides hope for improving forest resilience and avoiding forest tipping points.

The results therefore have implications for forest management and climate change mitigation strategies. In particular, the European Union’s biodiversity strategy for 2030 aims to build resilience to future threats such as climate change. Selective logging and uneven-aged management, which preserve or enhance canopy complexity, are likely to be more effective in maintaining and increasing forest resilience compared to homogenising practices like clear-cutting. More generally, monitoring, promoting and maintaining forest diversity is a policy to protect against the worst impacts of temperature rises in order to maintain resilient forests.

Methods

Vegation data

The analysis uses the kernel Normalized Difference Vegetation Index (kNDVI) as a proxy for the status and productivity of forests. The daily time series of Surface Reflectance at 500m spatial resolution from MODIS Terra (MOD09GA) and Aqua (MYD09GA) datasets, ranging from 2003 to 2021^69,70 are processed in Google Earth Engine (GEE)⁷¹ to retrieve a long-term time series of kNDVI. The kNDVI is a nonlinear generalisation of the Normalized Difference Vegetation Index (NDVI) showing stronger correlation with key forest parameters compared to NDVI and Near-Infrared Reflectance of Vegetation (NIRv).

The kNDVI is defined as²⁶

1

$$\text{kNDVI}=tanh\left({\left(\frac{\text{NIR}-\text{red}}{\sigma }\right)}^2\right)$$ In this work we define the length-scale parameter, σ, as

2

$$\sigma =\frac{1}{2}\left(\text{NIR}+\text{red}\right)$$ NIR is the Near-InfraRed band (MODIS band 2 at 841-876 nm), and red is the Red band (MODIS band 1 at 620-670 nm).

Only highest quality pixels are selected from daily MODIS data, where observations are masked for clouds, cloud shadows, water, snow and ice, by referring to the Quality Assessment (QA) bit flags. After quality filtering and after retaining only positive values of NDVI, daily kNDVI are computed and averaged into a time series with an 8-day time step and reprojected to 0.005^∘ spatial resolution.

Each pixel in the time series is then masked using a binary forest mask to ensure the time series is associated with a high percentage of forest cover. The forest mask is derived processing the University of Maryland Global Forest Change product version 1.9, at 1 arc-second spatial resolution (approximately 30 metres at the equator), as follows⁷². The tree cover layer for year 2000 is used to select clusters of at least 6 connected pixels containing ≥30% tree cover. These clusters are defined as ‘forest covered’ in the analysis. Any pixel undergoing a loss in the time period is removed. A forest loss in the dataset is defined as a stand-replacement disturbance, meaning a change from a forest to non-forest state. Forest loss pixels are removed in order to exclude from the analysis managed or disturbed forest patches. The identified clusters of forest covered pixels are used to compute the percentage of forest cover in the final grid at 0.005^∘ spatial resolution. This forest cover is masked to retain only 0.005^∘ grid cells with at least 50% of forest cover.

The masked kNDVI is then aggregated to 0.05^∘ resolution via the mean kNDVI value. For each pixel, the observations within the temporal distribution with a z-score below -2.5 are removed from the time series to account for potential residual cloud contaminations. This value is chosen due to highly asymmetric tails below -2.5 and above 2.5 in the z-score distribution, determined to result from cloud contamination. As a temporally-consistent kNDVI time series is required for the computation of resilience metrics, only pixels with a kNDVI time series where at least 200 values have another consecutive kNDVI value are considered (though on average pixels contain around 500 kNDVI values with a consecutive value). This finally results in a 0.05^∘ resolution time series of 8-days averaged kNDVI, derived from forested pixels.

The vegetation dynamics are assessed only within the growing season, focusing the analysis on the period in which forests are most responsive to external perturbations and excluding the influence of dormant months. Phenology is derived from MODIS Land Cover Dynamics yearly data at 500m⁷³. For each 500m pixel, mean greenup and dormancy dates are obtained for the time period 2001–2021. The resulting growing season product is aggregated to 0.05^∘ resolution, after accounting only for at least 50% forest cover at 0.005^∘ resolution.

Climate data

In order to explore the relationship between resilience and forest structural diversity (FSD), the effects of climate are taken into account. The climate variables considered build on those identified as relevant in the prediction of resilience by Forzieri et al.⁸. These are the total precipitation (Precip.), 2m air temperature (2m temp), surface net solar radiation (SSR), and Vapour Pressure Deficit (VPD). VPD is calculated as the difference between the saturated vapour pressure and the actual vapour pressure, and is an important factor in the regulation of stomatal conductance in plants⁷⁴.

These variables are obtained from the ERA5-Land reanalysis dataset^75,76. The ERA5-Land dataset is derived from ERA5 data, the fifth generation ECMWF reanalysis global climate and weather dataset at 0.25^∘ spatial resolution where the land component of the reanalysis is replayed to provide land variables at an enhanced resolution at 0.1^∘. For total precipitation, 2m temperature, and net solar radiation, the ERA5-Land hourly values for the period 2003-2021 are averaged to a daily mean across 8-day periods in synchronicity with the MODIS data. For VPD, the daily saturated vapour pressure, e_s, is calculated by⁷⁷

3

$$e_s=\left[e^{\circ }\left({\text{T}}_{\max }\right)+e^{\circ }\left({\text{T}}_{\min }\right)\right]/2\text{where:}{e}^{\circ }\left(\text{T}\right)=0.6108\text{x}{\exp }^{17.27\text{T}/\left(\text{T}+237.3\right)}$$ where ${\text{T}}_{\max }$ and ( ${\text{T}}_{\min }$ ) are the mean daily maximum and minimum air temperatures over each 8-day period. The latter formula is also used in the calculation of the actual vapour pressure from the dewpoint temperature. The ERA5-Land hourly reanalysis 2m air temperature and 2m dewpoint temperature data are used in the computation.

The ERA5-Land variables are interpolated from 0.1^∘ to 0.05^∘ spatial resolution, via bilinear interpolation in the case of the air and dewpoint temperature, and first-order conservative remapping for the precipitation and surface solar radiation. As with the kNDVI, the climate time series is only considered within the previously defined growing season. The final time series consists of 0.05^∘ resolution climate metrics averaged over 8 days over the period 2003 to 2021.

Auxiliary data

Four auxiliary datasets are also considered: digital elevation data, soil organic carbon content, nitrogen deposition data, and biogeographical regions. The digital elevation data are taken from the Shuttle Radar Topography Mission (SRTM) digital elevation data at 30m resolution and the standard deviation of elevation is computed at 0.05^∘ spatial resolution. This accounts for variation in the topography within a pixel, and in particular the effect this might have on the mean climate conditions⁷⁸. Soil organic carbon content is taken from the OpenLandMap Soil Organic Carbon Content dataset which has 250m spatial resolution and a depth of 30cm, and, after aggregating to 0.05^∘, accounts for differing soil quality between pixels⁷⁹. Nitrogen deposition data are taken from the EMEP Transboundary particulate matter, photo-oxidants, acidifying and eutrophying components 2024 status report and are available at 0.1^∘ resolution. The nitrogen deposition is taken as the total annual sum of reduced (wet and dry) nitrogen deposition, interpolated to 0.05^∘ and averaged over the period 2013-2018⁸⁰. Finally, the European Environment Agency dataset on biogeographical regions⁸¹ is processed in order to obtain three main biogeographic regions covering the area under analysis: alpine, temperate and mediterranean. This categorisation is visualised in Fig. 5b.

Forest structural diversity metrics

Forest structural diversity (FSD) metrics for horizontal (FSDH), vertical (FSDV) and combined horizontal and vertical (FSDH+V) diversity are derived from data provided by NASA’s GEDI instrument aboard the International Space Station. The native GEDI data provides approximately 25m-wide LiDAR measurements on forests between  ± ~ 50^∘ latitude³⁸. The derived metrics include the canopy cover fraction (CC), the percentage of the GEDI sample covered by a vertical projection of the tree crown, and forest Relative Height (RH) percentiles of the Vertical Profile (VP) of vegetation distribution (see Fig. 3. ref. ³⁸). RH98, for example, represents the forest height at which 98% of the reflected GEDI light is returned, relative to the ground, and is often used as a proxy for the canopy height⁸².

The GEDI metrics, collected between April 2019 and January 2023, are converted into FSD metrics via the methodology described in ref. ³⁹. In this study, we use the FSD dataset derived directly from GEDI data (i.e. the predictors rather than the machine learning-completed version). In summary, ’good’ quality flagged night-time, growing-season, GEDI observations are considered only on gentle (≤20^∘) slopes, with at least 10% tree and canopy cover and non-urbanised (≤5%) footprints without persistent water coverage. It is also required that the GEDI observations contain ground measurements. For consistency with the MODIS kNDVI data, a forest cover mask is also applied to GEDI observations during the FSD computation. The FSD metrics considered in this paper, motivated by metrics established in structural diversity mapping studies, are summarised in Table 1 and detailed below^64,83,84. Each metric is computed at a resolution of 0.05^∘, where each pixel is required to contain at least 40 valid GEDI observations.

The horizontal diversity metric, FSDH (τ_CH in ref. ³⁹), is calculated as the standard deviation of the canopy height (RH98) of the N GEDI observations within each 0.05^∘ pixel:

4

$$\text{FSDH}=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left[{\text{RH98}}_i-\mu \left(\text{RH98}\right)\right]}^2}$$ where μ is the mean. This metric describes the horizontal spatial diversity in the height of the canopy crown. Uncertainty in FSDH estimated using the standard error (SE) formula:

5

$$\text{SE}\left(\text{FSDH}\right)=\text{FSDH}/\sqrt{2N-2}$$

The vertical structural diversity, FSDV ( − τ_KU in ref. ³⁹), is defined as the inverse of the mean excess kurtosis in the pixel:

6

$$\text{FSDV}=-\frac{1}{N}\sum _i^N\left(\frac{\mathbb{E}\left[{\left({\text{RH}}_i-\mu \left({\text{RH}}_i\right)\right)}^4\right]}{{\left(\mathbb{E}\left[{\left({\text{RH}}_i-\mu \left({\text{RH}}_i\right)\right)}^2\right]\right)}^2}-3\right)$$ where RH_i is the above-ground Relative Height return energy distribution within the i^th GEDI sample shot, across a pixel with N GEDI samples. FSDV is a good approximation of the within-pixel inverse mean excess kurtosis of the above-ground vertical distribution of vegetation elements. A factor of 3 accounts for the excess kurtosis relative to a normal distribution. Excess kurtosis is a metric that captures the shape of the vegetation VP independently of tree height, unlike metrics such as the VP standard deviation. Higher positive FSDV values indicate a VP with greater vegetation presence in the upper and lower canopy, suggesting a more even vertical distribution of vegetation. On the other hand, lower or negative FSDV values signify a concentration of vegetation in the central part of the canopy.

The combined horizontal and vertical structural diversity, FSDH+V (τ_SW in ref. ³⁹), is defined from the Shannon Entropy⁸⁵. It is calculated by first representing each GEDI observation as a vector in a 4D Cartesian space, defined by the basis {RH50, RH75, RH98, CC(x10)}. This 4D space is then divided into 1-unit size bins, and entropy is computed based on the fraction of observations in the pixel that fall into each bin, p_α, capturing the spread and diversity of GEDI observations:

7

$$\text{FSDV+H}=-\sum _{\alpha }{p}_{\alpha }\log p_{\alpha }$$ where the summation is running over all of the bins in the 4D space. FSDV+H quantifies the variation among VPs within the pixel, with values ranging from zero, indicating VPs falling within a single bin, to increasingly positive values as the diversity among VPs within the pixel increases. It is therefore a measure of complexity and diversity in the substructure of the canopy across the pixel.

The uncertainties associated with FSDV and FSDH+V are estimated via a boostrapping approach: the metrics are re-calculated 100 times after (with-replacement) resampling in order to find the variability of each pixel’s excess kurtosis (FSDV) and shannon entropy (FSDH+V). The standard deviation is taken as the sampling uncertainty on each metric.

Maps showing the spatial distribution of the structural diversity metrics across Europe can be found in Fig. 1 and Supplementary Information A.

Data processing

The resulting datasets encompass a mix of vegetation and climate time series data over a temporal range from 2003 to 2021, and from which we extract a series of statistics, as well as soil, forest cover and topography data. They cover the spatial extent of 33.0 ≤ latitude ≤ 50, −10.7 ≤ longitude ≤ 44.6, which is the European area observed by GEDI, with a spatial resolution of 0.05^∘. At mid-latitudes, the selected European forests ensure comparison of similar forest species and are particularly suitable for study due to the challenges of inferring resilience at high and tropical latitudes, the limits of spaceborne lidar in measuring the high-latitudes, and as areas of high forest management²⁵.

From the time series data, we compute the mean and the coefficient of variation (C.V.), the latter expressed as the ratio between the standard deviation and the mean. These respectively provide information on the climate and vegetation background state and climate variability.

The temporal autocorrelation of perturbations from the mean state is then computed. In order to remove the dominant seasonal component of the vegetation cycle, both the kNDVI and the climate time series are deseasonalised by subtracting from each timestep the long-term mean (2003-2021) of each 8 days timestep. Long-term trends, potentially driven by climate change, are also removed from the deseasonalised time series. A linear regression model is fitted to each pixel’s time series and the linear trend is subtracted from individual kNDVI and climate observations. The deseasonalised and detrended vegetation and climate time series are used to compute the lag-1 temporal auto-correlation (AC1) and (in the case of the kNDVI time series only) the variance.

In order to remove the dominant seasonal component of the vegetation cycle, both the kNDVI and the climate time series are deseasonalised by subtracting from each timestep the long-term mean (2003-2021) of each 8 days timestep. Long-term trends, potentially driven by climate change, are also removed from the deseasonalised time series. A linear regression model is fitted to each pixel’s time series and the linear trend is subtracted from individual kNDVI and climate observations. The deseasonalised and detrended vegetation and climate time series are used to compute the lag-1 temporal auto-correlation (AC1) and (in the case of the kNDVI time series only) the variance.

Resilience metrics

Engineering resilience is interpreted as the rate at which a system is restored following a perturbation¹⁶. It can be characterised by exponential decay of the form x(t) ≈ x₀e^λt, where x(t) is the perturbation of the system about an equilibrium value at a given time t, and λ is the decay rate or restoration rate (here defined such that negative values are equilibrium-restoring)¹⁰.

If we include a stochastic forcing term, then the system physics describes an equilibrium-reverting random drift, known as an Ornstein-Uhlenbeck process^86,87. In a discrete time framework, the system displacement from equilibrium, between consecutive timesteps (t_n and t_n+1; defined such that t_n+1 − t_n = Δt = 1) can be written as^21,25

8

$$x\left(t_{n+1}\right)=\alpha x\left(t_n\right)+\sigma {ϵ}_n$$ where α is the autoregressive coefficient, which determines the strength of the dependence on the previous value and is directly related with the restoration rate, λ, and ϵ refers to the stochastic displacement with standard deviation σ. From this, it follows that the 1-lag autocorrelation (AC1) between consecutive timesteps is^86,88

9

$$\text{AC1}\left[x\left(t_{n+1}\right)\right]=\alpha =e^{\lambda }$$ and the variance is^86,88

10

$$V\left[x\left(t_{n+1}\right)\right]=\alpha E\left[x{\left(t_n\right)}^2\right]+E\left[\sigma ϵ\left(t\right)\right]={\sigma }^2/\left(1-{\alpha }^2\right)$$ If restoration rate approaches zero, the autocorrelation and variance tend towards unity and infinity respectively. This slowness in the return of the system to its equilibriums state is defined as critical slowing down. The two indicators act as early-warning signals for critical transitions as the system approaches a bifurcation point between two states^18,19.

Two comparative estimations of the recovery or restoration rate, λ, can then be defined from the AC1 and the variance^25,40:

11

$$\text{Rest. Rate AC1}=\mid ln\left(\alpha \right)\mid$$

12

$$\text{Rest}.\mathrm{Rate}\mathrm{Variance}=∣\frac{1}{2}ln\left(1-\frac{{\sigma }^2}{\text{V}\left[x\right]}\right)∣$$ The 1-lag autocorrelation, α, and standard deviation of the noise term, σ, can be computed by regressing the time series of perturbations in formula (8). The small number of pixels where the logarithms cannot be defined (α < 0 or $\frac{{\sigma }^2}{V\left[x\right]}>1$ ) are removed. These result from a mixture of high noise pixels and (particularly at high latitudes) high σ pixels, likely due to short growing seasons²⁵. Note that the restoration rates in this analysis are defined in terms of absolute values in equations (11) and (12), such that the rate is interpreted as the magnitude of the restoring force from the perturbation: the larger the restoration rate, the faster the system recovers.

In this paper, x(t) represents the vegetation productivity about its seasonal mean. This is defined by the pixel-wise deseasonalised, detrended kNDVI time series, with the climate and stochastic processes driving perturbations about this equilibrium. The first restoration rate, equation (11), is defined in terms of the autocorrelation of consecutive kNDVI displacements about the mean, and the second, equation (12), in terms of the variance of displacements about the mean and the variance of the residuals of a linear model between consecutive displacements.

Decreasing autocorrelation and variance of each forest pixel about the seasonal mean corresponds to an increasing restoration rate magnitude and therefore an increasing resilience. Through the comparison of previous time steps, the restoration rate AC1 metric relates more to the system memory, whilst the variance-derived restoration rate metric is associated more with the stability of forest productivity²¹. The two separately calculated metrics act as a cross check on the other, corroborating or refuting potential trends and patterns. Maps showing the distribution of the restoration rates calculated in the analysis can be found in Fig. 1 and supplementary information A.

The Random Forest model

A number of environmental factors have been found to drive spatial differences in the temporal autocorrelation and variance of kNDVI time series⁸. Aridity, in particular, is observed as a key driver of differences in the recovery of vegetation from perturbations between regions^29,30. The complex interactions between multiple environmental drivers, makes comparing the autocorrelation and variance of spatially separated regions more complicated than temporal changes across environmentally coherent regions.

In order to explore potential drivers of the resilience metrics across a geographically diverse region, such as Europe, a random forest (RF) regression model is chosen, a non-linear, non-parametric machine learning model which does not assume any form of interaction between drivers and response^89,90. The response variable of the RF models are the resilience metrics, the AC1-derived and variance-derived restoration rates, with separate models for each. Each model takes as predictors the previously described variables of mean background climate and kNDVI, the C.V. and the AC1 of the climate time series, as well as the stationary maps of soil carbon content, nitrogen deposition, forest cover, topographical variation and one of the three FSD metrics per model. The 18 drivers, summarised in Supplementary Information H, therefore incorporate vegetation properties, soil quality, and climate average, variability and memory. The paper shows the results of 6 separate RF models optimised for the different combinations of two resilience metrics and three FSD metrics.

The dataset of 77,797 pixels at 0.05 spatial resolution with a full set of driver and response values is divided into the same 70% training and 30% testing set split for each model. The Random Forest model comprises 500 trees, with 7 variables sampled at each split, as determined by a 10-fold cross-validation optimisation process. Uncertainty in the models is estimated by bootstrapping a further 100 samples from the training dataset, with the FSD metric values in each model randomly sampled from a normal distribution with standard deviation as the estimated associated error. The uncertainties associated with the limited sampling of the FSD metrics are expected to dominate the overall uncertainty, and therefore uncertainties associated with non-remotely sensed inputs (such as the ERA5 climate model variables) are not included in the analysis. Supplementary Information G shows how the RF model accounts for aridity, one of the dominant environmental drivers of vegetation resilience, demonstrating that the model residuals show no aridity gradient.

Supplementary Information B shows metrics for the performance of the different models in replicating the observed restoration rates. The section also displays the importance of the predictor variables for each model, as measured through the increase in the MSE of the model when each predictor is permuted. Supplementary information B also displays equivalent simple linear and multivariate regressions produced in support of the RF models. Due to the stochastic nature of RF models, some variability in the predictor importance is expected, therefore, the figures show the mean and 95% confidence interval of the bootstrapped models. Whilst the RF model does not provide information on the underlying physics and mechanisms relating the variables, the outputs may give insights to the relationships.

In order to support the analysis and assess the impact of resolution on the model relationships, the random forest inputs are also recomputed at 0.1^∘ and the analysis re-run (without computation of the corresponding uncertainties). The results are shown in Supplementary Information E.

The relationship between the resilience metrics and the structural diversity metrics

The relationship between the metrics of resilience and the metrics for structural diversity is analysed through Partial Dependence Plots (PDPs) and individual conditional expectation (ICE)⁹¹. These methods allow for interpretation of complex non-parametric models by relating the model output with predictors of interest.

PDPs show the global relationship between one or two predictors and the predicted variable, and therefore display the marginal effect of the variation of the input structural diversity metrics on the predicted value of the metrics of resilience, whilst controlling for the other variables. The R package, pdp⁹², is used to compute the partial dependence, f_s(x_s), of the output on a given set of input variables x_s. This is defined in terms of the prediction function ${\widehat{f}}_s\left(x\right)$ of the RF input variables x (with x_c the controlled variables)⁹³:

13

$$f_s\left(x_s\right)=E_{x_c}\left[\widehat{f}\left(x_s,x_c\right)\right]=\int \widehat{f}\left(x_s,x_c\right)p_c\left(x_c\right)\mathrm{d}{x}_c$$ where the marginal probability density of x_c is p_c(x_c) = ∫p(x) dx_s. Providing the relationships between the variables is not too strong (Supplementary Information I), the partial dependence can be estimated from the training dataset by averaging out over the controlled variables x_c such that

14

$$f_s\left(x_s\right)=\frac{1}{n}\sum _{i=1}^n\widehat{f}\left(x_s,x_{i,c}\right)$$ where n is the size, and i is the i^th entry of the training data. In effect, we obtain a function over the entire series of x_i,s values in the training dataset, where the value at each position, j, along the variable of interest, is the average of the prediction function applied to a duplicated training dataset except with the values of x_s replaced by x_j,s (i.e. the values of x_s in the duplicated dataset are the same). The calculation of partial dependence implies a causal effect within the model through the adjustment of predictors and measurement of the effect; however, does not necessarily require that the real-world effect is causational⁹⁴.

In the analysis presented, the partial dependence of the response variables, the resilience metrics, are first shown as a function of the metrics of structural diversity, controlling for the other RF model inputs, and the results are shown in Fig. 3. The response of the metrics of resilience are then shown as a function of one of each of the metrics of structural diversity and the mean daily growing season temperature, in Fig. 5, in order to show the combined relationship between resilience, temperature and structural diversity.

Whilst PDPs take global averages of the relationship between metrics, ICE figures focus on the local level relationship between predictor and response, visualising the local relationship at each data point or pixel⁹⁵. In summary, the partial dependence is therefore the average of the local level ICE relationships, whilst the ICE plots are disaggregations of the PDPs. The R package, ICEbox⁹⁶, is used to estimate the response function curve, f_s(x_s), for each data point or pixel, in a similar methodology to that described for the PDP. At each point, f_s(x_s) is estimated from the prediction function, $\widehat{f}\left(x_s,x_c\right)$ , applied to the values of the controlled variables, x_c, at that point, and the full range of the variable of interest across the dataset x_s as the dependent variable.

In order to improve the visualisation of the ICE curve for each pixel, the partial derivative of each curve with respect to the variable of interest is taken, after the application of a smoothing function. This enables the production of a map, Fig. 4, showing the derivative ICE (or dICE) value, for the relationship between the different resilience metrics and the different structural diversity metrics. The sign of the derivative shows the direction of the relationship, whilst the absolute value shows the magnitude of the relationship relative to other pixels. If there is no interaction between the feature of interest and the other features, then the partial derivative is approximately the same for all points, and thus variability in the dICE values informs us about pixel-level deviations from the global relationship.

Acknowledgements

The study was co-funded by the Exploratory Project ForBioRes of the European Commission, Joint Research Centre and by the EU-H2020 project FORGENIUS (Grant Agreement 86221). The views expressed are purely those of the writers and may in no circumstance be regarded as stating an official position of the European Commission. This is ClimTip contribution #88; the ClimTip project has received funding from the European Union's Horizon Europe research and innovation programme under grant agreement No. 101137601.

Author contributions

M. Girardello, A. Cescatti and M. Pickering conceived the original idea. M. Pickering, A. Elia, M. Girardello, G. Oton, G. Forzieri, M. Migliavacca, and A. Cescatti designed the study. M. Pickering, A. Elia, M. Girardello, G. Oton, G. Ceccherini, and M. Piccardo collected the data and conducted the analyses. M. Pickering and A. Elia led the writing with contributions from all authors. All authors contributed to discussions regarding the interpretation of the results and revisions of the manuscript.

Peer review

Peer review information

Communications Earth & Environment thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editors: Mengjie Wang. A peer review file is available.

Data availability

The datasets used in this study are publicly available^{69, 70, 71, 72–73,76,78, 79, 80–81}. The forest structural diversity metrics³⁹ are available here: https://figshare.com/s/daa9b652c12beb42e518.

Code availability

The custom code written to download and analyse the data and generate figures, is available here: https://github.com/markpickering1/diversity_resilience_analysis.git.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s43247-025-02592-8.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Friedlingstein, P et al. Global carbon budget 2021. Earth Syst. Sci. Data; 2022; 14, pp. 1917-2005. [DOI: https://dx.doi.org/10.5194/essd-14-1917-2022]

2. Shukla, P. R. et al. (eds.) Climate Change and Land: An IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems (IPCC, 2019).

3. Pielke, RA et al. Interactions between the atmosphere and terrestrial ecosystems: influence on weather and climate. Glob. Change Biol.; 1998; 4, pp. 461-475. [DOI: https://dx.doi.org/10.1046/j.1365-2486.1998.t01-1-00176.x]

4. Korosuo, A. et al. The role of forests in the EU climate policy: are we on the right track? Carbon Balance Manag.18 (2023).

5. Commission, E. & for Communication, D.-G.The European Green Deal—Delivering the EU’s 2030 climate targets (Publications Office of the European Union, 2023).

6. Parliament, E., for Internal Policies of the Union, D.-G. & Negre, F. The EU 2030 biodiversity strategy (European Parliament, 2020).

7. Roebroek, CT et al. Climate policies for carbon neutrality should not rely on the uncertain increase of carbon stocks in existing forests. Environ. Res. Lett.; 2024; 19, 044050. [DOI: https://dx.doi.org/10.1088/1748-9326/ad34e8]

8. Forzieri, G; Dakos, V; McDowell, N; Alkama, R; Cescatti, A. Emerging signals of declining forest resilience under climate change. Nature; 2022; 608, pp. 1-6. [DOI: https://dx.doi.org/10.1038/s41586-022-04959-9]

9. Boulton, C., Lenton, T. & Boers, N. Pronounced loss of Amazon rainforest resilience since the early 2000s. Nat. Clim. Change12 (2022).

10. Smith, T., Traxl, D. & Boers, N. Empirical evidence for recent global shifts in vegetation resilience. Nat. Clim. Change12 (2022).

11. Patacca, M et al. Significant increase in natural disturbance impacts on European forests since 1950. Glob. Change Biol.; 2023; 29, pp. 1359-1376. [DOI: https://dx.doi.org/10.1111/gcb.16531] 1:CAS:528:DC%2BB38XjtFWqsrbP

12. Seidl, R; Schelhaas, M-J; Rammer, W; Verkerk, H. Increasing forest disturbances in Europe and their impact on carbon storage (vol 4, pg 806, 2014). Nat. Clim. Change; 2014; 4, pp. 930-930. [DOI: https://dx.doi.org/10.1038/nclimate2393]

13. Forzieri, G et al. Emergent vulnerability to climate-driven disturbances in European forests. Nat. Commun.; 2021; 12, [DOI: https://dx.doi.org/10.1038/s41467-021-21399-7] 1:CAS:528:DC%2BB3MXltFOmsb4%3D 1081.

14. Forzieri, G. et al. Ecosystem heterogeneity is key to limiting the increasing climate-driven risks to European forests. One Earth (2024).

15. Bengtsson, J; Nilsson, SG; Franc, A; Menozzi, P. Biodiversity, disturbances, ecosystem function and management of european forests. For. Ecol. Manag.; 2000; 132, pp. 39-50. [DOI: https://dx.doi.org/10.1016/S0378-1127(00)00378-9]

16. Pimm, S. The complexity and stability of ecosystems. Nature; 1984; 307, pp. 321-326. [DOI: https://dx.doi.org/10.1038/307321a0]

17. of Engineering, N. A. Engineering Within Ecological Constraints (The National Academies Press, 1996).

18. Dakos, V et al. Tipping point detection and early warnings in climate, ecological, and human systems. Earth Syst. Dyn.; 2024; 15, pp. 1117-1135. [DOI: https://dx.doi.org/10.5194/esd-15-1117-2024]

19. Bathiany, S. et al. Ecosystem resilience monitoring and early warning using earth observation data: challenges and outlook. Surv. Geophys.46, 265–301 (2024).

20. Dakos, V et al. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PLoS ONE; 2012; 7, pp. 1-20. [DOI: https://dx.doi.org/10.1371/journal.pone.0041010]

21. Scheffer, M et al. Early-warning signals for critical transitions. Nature; 2009; 461, pp. 53-9. [DOI: https://dx.doi.org/10.1038/nature08227] 1:CAS:528:DC%2BD1MXhtVygsbjI

22. Dakos, V; Carpenter, SR; van Nes, EH; Scheffer, M. Resilience indicators: prospects and limitations for early warnings of regime shifts. Philos. Trans. R. Soc. B: Biol. Sci.; 2015; 370, 20130263. [DOI: https://dx.doi.org/10.1098/rstb.2013.0263]

23. Verbesselt, J. et al. Remotely sensed resilience of tropical forests. Nat. Clim. Change6 (2016).

24. Hirota, M; Holmgren, M; Nes, EHV; Scheffer, M. Global resilience of tropical forest and savanna to critical transitions. Science; 2011; 334, pp. 232-235. [DOI: https://dx.doi.org/10.1126/science.1210657] 1:CAS:528:DC%2BC3MXht1yktr7O

25. Smith, T. & Boers, N. Reliability of vegetation resilience estimates depends on biomass density. Nat. Ecol. Evol.7, 1799–1808 (2023).

26. Camps-Valls, G et al. A unified vegetation index for quantifying the terrestrial biosphere. Sci. Adv.; 2021; 7, eabc7447. [DOI: https://dx.doi.org/10.1126/sciadv.abc7447] 1:CAS:528:DC%2BB3MXpsVartrc%3D

27. Feng, Y. et al. Reduced resilience of terrestrial ecosystems locally is not reflected on a global scale. Commun. Earth Environ.2 (2021).

28. Liu, Y; Kumar, M; Katul, G; Porporato, A. Reduced resilience as an early warning signal of forest mortality. Nat. Clim. Change; 2019; 9, pp. 1-6. [DOI: https://dx.doi.org/10.1038/s41558-019-0583-9]

29. De Keersmaecker, W et al. A model quantifying global vegetation resistance and resilience to short-term climate anomalies and their relationship with vegetation cover. Glob. Ecol. Biogeogr.; 2015; 24, pp. 539-548. [DOI: https://dx.doi.org/10.1111/geb.12279]

30. Smith, T. & Boers, N. Global vegetation resilience linked to water availability and variability. Nat. Commun.14 (2023).

31. LaRue, EA et al. Structural diversity as a reliable and novel predictor for ecosystem productivity. Front. Ecol. Environ.; 2023; 21, pp. 33-39. [DOI: https://dx.doi.org/10.1002/fee.2586]

32. Hooper, DU et al. Effects of biodiversity on ecosystem functioning: a consensus of current knowledge. Ecol. Monogr.; 2005; 75, pp. 3-35. [DOI: https://dx.doi.org/10.1890/04-0922]

33. Mahecha, MD et al. Biodiversity and climate extremes: Known interactions and research gaps. Earth’s. Future; 2024; 12, e2023EF003963. [DOI: https://dx.doi.org/10.1029/2023EF003963]

34. LaRue, E et al. A theoretical framework for the ecological role of three dimensional structural diversity. Front. Ecol. Environ.; 2023; 21, pp. 4-13. [DOI: https://dx.doi.org/10.1002/fee.2587]

35. Guo, X et al. Regional mapping of vegetation structure for biodiversity monitoring using airborne lidar data. Ecol. Inform.; 2017; 38, pp. 50-61. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2017.01.005]

36. Ibáñez, I et al. Forest resilience under global environmental change: Do we have the information we need? a systematic review. PLoS ONE; 2019; 14, e0222207. [DOI: https://dx.doi.org/10.1371/journal.pone.0222207]

37. Thompson, I., Mackey, B., Mcnulty, S. & Mosseler, A. Forest Resilience, Biodiversity, and Climate Change. A Synthesis of the Biodiversity/Resilience/Stability Relationship in Forest Ecosystems (UNEP, 2009).

38. Dubayah, R et al. The global ecosystem dynamics investigation: High-resolution laser ranging of the earth’s forests and topography. Sci. Remote Sens.; 2020; 1, 100002. [DOI: https://dx.doi.org/10.1016/j.srs.2020.100002]

39. Girardello, M et al. A dataset on the structural diversity of european forests. Earth Syst. Sci. Data Discuss.; 2025; 2025, pp. 1-30.

40. Dakos, V; van Nes, EH; D’Odorico, P; Scheffer, M. Robustness of variance and autocorrelation as indicators of critical slowing down. Ecology; 2012; 93, pp. 264-271. [DOI: https://dx.doi.org/10.1890/11-0889.1]

41. Ehbrecht, M et al. Global patterns and climatic controls of forest structural complexity. Nat. Commun.; 2021; 12, [DOI: https://dx.doi.org/10.1038/s41467-020-20767-z] 1:CAS:528:DC%2BB3MXit1alsLY%3D 519.

42. García-Tejero, S; Spence, JR; O’Halloran, J; Bourassa, S; Oxbrough, A. Natural succession and clearcutting as drivers of environmental heterogeneity and beta diversity in north american boreal forests. PLoS ONE; 2018; 13, pp. 1-16. [DOI: https://dx.doi.org/10.1371/journal.pone.0206931]

43. (C3S), C. C. C. S. (ed.) European State of the Climate 2023.

44. Chen, N., Jayaprakash, C., Yu, K. & Guttal, V. Rising variability, not slowing down, as a leading indicator of a stochastically driven abrupt transition in a dryland ecosystem. Am. Naturalist191, E1–E14 (2018).

45. Neijnens, F; Siteur, K; van de Koppel, J; Rietkerk, M. Early warning signals for rate-induced critical transitions in salt marsh ecosystems. Ecosystems; 2021; 24, pp. 1-12. [DOI: https://dx.doi.org/10.1007/s10021-021-00610-2]

46. Huang, Y; Bathiany, S; Boers, N. Deep learning for predicting rate-induced tipping. Nat. Mach. Intell.; 2024; 6, pp. 1556-1565. [DOI: https://dx.doi.org/10.1038/s42256-024-00937-0]

47. Lloret, F; Calvo, E; Pons, X; Diaz-Delgado, R. Wildfires and landscape patterns in the eastern iberian peninsula. Landsc. Ecol.; 2002; 17, pp. 745-759. [DOI: https://dx.doi.org/10.1023/A:1022966930861]

48. Weir, J., Johnson, E. & Miyanishi, K. Fire frequency and the spatial age mosaic of the mixed-wood boreal forest in western canada. Ecol. Appl.10 (2000).

49. Seidl, R; SCHELHAAS, M-J; Lexer, M. Unraveling the drivers of intensifying forest disturbance regimes in europe. Glob. Change Biol.; 2011; 17, pp. 2842-2852. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2011.02452.x]

50. Mitchell, S. Wind as a natural disturbance agent in forests: a synthesis. Forestry: Int. J. For. Res.; 2012; 86, pp. 147-157. [DOI: https://dx.doi.org/10.1093/forestry/cps058]

51. Marini, L., Ayres, M. & Jactel, H. Impact of stand and landscape management on forest pest damage. Annu. Rev. Entomol. 67 (2022).

52. Rosenfeld, JS. Functional redundancy in ecology and conservation. Oikos; 2002; 98, pp. 156-162. [DOI: https://dx.doi.org/10.1034/j.1600-0706.2002.980116.x]

53. Schindler, DE; Armstrong, JB; Reed, TE. The portfolio concept in ecology and evolution. Front. Ecol. Environ.; 2015; 13, pp. 257-263. [DOI: https://dx.doi.org/10.1890/140275]

54. Ratcliffe, S et al. Biodiversity and ecosystem functioning relations in european forests depend on environmental context. Ecol. Lett.; 2017; 20, pp. 1414-1426. [DOI: https://dx.doi.org/10.1111/ele.12849]

55. Anderegg, WRL et al. Hydraulic diversity of forests regulates ecosystem resilience during drought. Nature; 2018; 561, pp. 538-541. [DOI: https://dx.doi.org/10.1038/s41586-018-0539-7] 1:CAS:528:DC%2BC1cXhslKjsLzP

56. Liu, D., Wang, T., Penuelas, J. & Piao, S. Drought resistance enhanced by tree species diversity in global forests. Nat. Geosci.15 (2022).

57. Meinen, C; Hertel, D; Leuschner, C. Root growth and recovery in temperate broad-leaved forest stands differing in tree species diversity. Ecosystems; 2009; 12, pp. 1103-1116. [DOI: https://dx.doi.org/10.1007/s10021-009-9271-3]

58. Pardos, M et al. The greater resilience of mixed forests to drought mainly depends on their composition: analysis along a climate gradient across europe. For. Ecol. Manag.; 2021; 481, 118687. [DOI: https://dx.doi.org/10.1016/j.foreco.2020.118687]

59. Jucker, T; Bouriaud, O; Avacaritei, D; Coomes, DA. Stabilizing effects of diversity on aboveground wood production in forest ecosystems: linking patterns and processes. Ecol. Lett.; 2014; 17, pp. 1560-1569. [DOI: https://dx.doi.org/10.1111/ele.12382]

60. Schmitt, S et al. Functional diversity improves tropical forest resilience: Insights from a long-term virtual experiment. J. Ecol.; 2020; 108, pp. 831-843. [DOI: https://dx.doi.org/10.1111/1365-2745.13320]

61. Sakschewski, B et al. Resilience of amazon forests emerges from plant trait diversity. Nat. Clim. Change; 2016; 6, pp. 1032-1036. [DOI: https://dx.doi.org/10.1038/nclimate3109]

62. Mori, AS; Furukawa, T; Sasaki, T. Response diversity determines the resilience of ecosystems to environmental change. Biol. Rev.; 2013; 88, pp. 349-364. [DOI: https://dx.doi.org/10.1111/brv.12004]

63. Bachofen, C et al. Tree water uptake patterns across the globe. N. Phytologist; 2024; 242, pp. 1891-1910. [DOI: https://dx.doi.org/10.1111/nph.19762]

64. Whittaker, RH. Evolution and measurement of species diversity. Taxon; 1972; 21, pp. 213-251. [DOI: https://dx.doi.org/10.2307/1218190]

65. Hakkenberg, CR et al. Inferring alpha, beta, and gamma plant diversity across biomes with gedi spaceborne lidar. Environ. Res.: Ecol.; 2023; 2, 035005.

66. Potterf, M et al. Contrasting Norway spruce disturbance dynamics in managed forests and strict forest reserves in Slovakia. Forestry: Int. J. For. Res.; 2022; 96, pp. 387-398. [DOI: https://dx.doi.org/10.1093/forestry/cpac045]

67. Pausas, JG; Keeley, JE. Wildfires and global change. Front. Ecol. Environ.; 2021; 19, pp. 387-395. [DOI: https://dx.doi.org/10.1002/fee.2359]

68. Thom, D. Natural disturbances as drivers of tipping points in forest ecosystems under climate change—implications for adaptive management. Forestry: Int. J. For. Res.; 2023; 96, pp. 305-315. [DOI: https://dx.doi.org/10.1093/forestry/cpad011]

69. Vermote, E. & Wolfe, R. Modis/terra surface reflectance daily L2G global 1km and 500m sin grid v061 [data set]. NASA EOSDIS Land Processes Distributed Active Archive Center (2021). https://doi.org/10.5067/MODIS/MOD09GA.061 (Accessed 4 April 2024).

70. Vermote, E. & Wolfe, R. MODIS/Aqua Surface Reflectance Daily L2G Global 1km and 500m SIN Grid V061. NASA EOSDIS Land Processes Distributed Active Archive Center (2021). https://doi.org/10.5067/MODIS/MYD09GA.061 (Accessed 3 May 2024).

71. Gorelick, N et al. Google earth engine: planetary-scale geospatial analysis for everyone. Remote Sens. Environ.; 2017; 202, pp. 18-27. [DOI: https://dx.doi.org/10.1016/j.rse.2017.06.031] Big Remotely Sensed Data: tools, applications and experiences.

72. Hansen, MC et al. High-resolution global maps of 21st-century forest cover change. Science; 2013; 342, pp. 850-853. [DOI: https://dx.doi.org/10.1126/science.1244693] 1:CAS:528:DC%2BC3sXhslCrsrbO

73. Friedl, M., Gray, J. & Sulla-Menashe, D. Modis/terra+aqua land cover dynamics yearly l3 global 500m sin grid v061 [data set]. NASA EOSDIS Land Processes Distributed Active Archive Center (2022). https://doi.org/10.5067/MODIS/MCD12Q2.061 (Accessed 4 April 2024).

74. Grossiord, C et al. Plant responses to rising vapor pressure deficit. N. Phytologist; 2020; 226, pp. 1550-1566. [DOI: https://dx.doi.org/10.1111/nph.16485]

75. Muñoz Sabater, J et al. Era5-land: a state-of-the-art global reanalysis dataset for land applications. Earth Syst. Sci. Data; 2021; 13, pp. 4349-4383. [DOI: https://dx.doi.org/10.5194/essd-13-4349-2021]

76. Muñoz Sabater, J. ERA5-Land hourly data from 1950 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). https://doi.org/10.24381/cds.e2161bac (2019). (Accessed 13 September 2022).

77. Allan, R. & Pereira, L.Crop evapotranspiration—Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56, Vol. 56 (FAO, 1998).

78. Farr, T. G. et al. The shuttle radar topography mission. Rev. Geophys.45https://doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 (2007).

79. Hengl, T. & Wheeler, I. Soil organic carbon content in x 5 g / kg at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution (version v02) [data set]. Zenodo (2018).

80. MSC-W, J., CCC, CEIP & Report, C. Emep status report 1/2024: Transboundary particulate matter, photo-oxidants, acidifying and eutrophying components. https://www.emep.int/mscw/mscw_moddata.html (2024).

81. (EEA), E. E. A. Biogeographical regions, europe 2016, version 1 (2016) (Accessed 28 May 2024).

82. Ngo, Y.-N., Ho Tong Minh, D., Baghdadi, N. & Fayad, I. Tropical forest top height by GEDI: from sparse coverage to continuous data. Remote Sensing15 (2023).

83. Hakkenberg, CR; Goetz, SJ. Climate mediates the relationship between plant biodiversity and forest structure across the united states. Glob. Ecol. Biogeogr.; 2021; 30, pp. 2245-2258. [DOI: https://dx.doi.org/10.1111/geb.13380]

84. Ma, Q. et al. The coordinated impact of forest internal structural complexity and tree species diversity on forest productivity across forest biomes. Fundamental Research4, 1185–1195 (2024).

85. Ortiz-Burgos, S. Shannon-Weaver Diversity Index 572–573 (Springer Netherlands, 2016).

86. Dijkstra, H. A. Nonlinear climate dynamics. (Cambridge University Press, 2013) https://doi.org/10.1017/CBO9781139034135.

87. Gardiner, C. Handbook of stochastic methods for physics, chemistry, and the natural sciences (1985) https://doi.org/10.1007/978-3-662-02452-2.

88. Ives, AR. Measuring resilience in stochastic systems. Ecol. Monogr.; 1995; 65, pp. 217-233. [DOI: https://dx.doi.org/10.2307/2937138]

89. Breiman, L. Random forests. Machine Learning45 (2001).

90. Liaw, A; Wiener, M. Classification and regression by randomforest. R. N.; 2002; 2, pp. 18-22.

91. Molnar, C. Interpretable Machine Learning 2nd edn (Shroff Publishers, 2022).

92. Greenwell, BM. pdp: An R package for constructing partial dependence plots. R. J.; 2017; 9, pp. 421-436. [DOI: https://dx.doi.org/10.32614/RJ-2017-016]

93. Friedman, JH. Greedy function approximation: A gradient boosting machine. Ann. Stat.; 2001; 29, pp. 1189-1232. [DOI: https://dx.doi.org/10.1214/aos/1013203451]

94. Zhao, Q. & Hastie, T. Causal interpretations of black-box models. J Bus Econ Stat39 (2019).

95. Alex Goldstein, JB; Kapelner, A; Pitkin, E. Peeking inside the black box: Visualizing statistical learning with plots of individual conditional expectation. J. Computational Graph. Stat.; 2015; 24, pp. 44-65. [DOI: https://dx.doi.org/10.1080/10618600.2014.907095]

96. Goldstein, A; Kapelner, A; Bleich, J; Pitkin, E. Peeking inside the black box: Visualizing statistical learning with plots of individual conditional expectation. J. Computational Graph. Stat.; 2015; 24, pp. 44-65. [DOI: https://dx.doi.org/10.1080/10618600.2014.907095]