Main

Natural ecosystems are severely threatened by climate change and biodiversity loss; the Amazon, African and southeast Asian rainforests are key examples that have attracted substantial recent attention1,2,3. These tropical vegetation systems have been inferred to exhibit multistability for broad ranges of mean annual precipitation4,5; within the same precipitation ranges, both the rainforest state and an alternative savannah state are simultaneously stable. This implies that, even absent long-term changes in local or regional precipitation, transitions from the current rainforest state to the savannah state are possible and may be triggered by external perturbations such as droughts, forest fires and deforestation6. Although ecosystem transitions in tropical rainforests have received widespread attention, the risk of transitions to alternative ecosystem states appears to be a global characteristic that extends to high-latitude7,8 and dryland ecosystems9. Given that ecosystem transitions could turn net carbon sinks into carbon sources3 and the tremendous potential of vegetation to reduce atmospheric carbon dioxide concentrations10, the mitigation of anthropogenic climate change and the maintenance of global biodiversity are strongly dependent on the resilience of vegetation systems worldwide.

Ecosystem resilience is typically defined as the capacity to resist and recover from external disturbances11,12,13. Unfortunately, this definition only allows for the empirical measurement of resilience either in controlled experiments (by applying an artificial disturbance) or by waiting for occurrences of large external disturbances to natural vegetation systems. Due to the scarcity of suitably strong external perturbations, it is difficult to quantify the resilience of natural ecosystems at a global scale, and in particular to investigate resilience changes over time.

Theoretically, the fluctuation–dissipation theorem (FDT) from statistical mechanics14,15,16,17 suggests that for specific classes of systems, the response to external perturbations can be expressed in terms of the characteristics of natural fluctuations around the equilibrium state. In other words, the FDT states that the rate at which a system will return to equilibrium following an external disturbance can be determined from its internal natural fluctuations. The tremendous practical value of the FDT comes from the fact that, if it can be shown to hold for a given system, the response to external perturbations can be predicted on the basis of the internal variability of the system in question. Evidence that the FDT holds has been revealed in several real-world systems17, ranging from financial market data18,19 to atmospheric and climate dynamics20,21.

Several studies have suggested that the lag-one autocorrelation (AC1)—a measure of how strongly correlated neighbouring time spans of a given time series are—and variance of a system can be used as measures of vegetation resilience1,22,23,24,25,26,27. The variability of natural fluctuations can be estimated in terms of the variance22,27,28, while the strength of the system’s memory can be measured using the AC11,23,24,25,28. Low-dimensional dynamical system frameworks and designed experiments justify this choice by showing that variance and AC1 increase as the system approaches a critical threshold beyond which a bifurcation-induced transition—a jump to an alternative stable state—occurs, which is interpreted as a loss of resilience29,30. The increase in AC1 together with a corresponding increase in variance have been termed early-warning signals for critical transitions; the underlying change in dynamics is referred to as ‘critical slowing down’22,28. It has been shown that early-warning signals can be identified before abrupt climate transitions evidenced in palaeoclimate records31,32,33 as well as in ecosystem28 and climate34,35 model simulations. However, although the AC1 and variance have been used to quantify the stability or resilience of different systems, their actual suitability as measures of ecosystem, and in particular vegetation, resilience has not been confirmed outside of controlled and model-based experiments36,37, and in particular not based on empirical evidence.

In this article, we use empirical remotely sensed vegetation data to test for the correspondence between theoretical vegetation resilience—AC1 and variance—and the rates of recovery from perturbations. We first use large perturbations to derive empirical recovery rates for diverse landscapes, vegetation types and climate zones using two independent vegetation datasets based on optical (advanced very-high-resolution radiometer (AVHRR) normalized difference vegetation index (NDVI), 1981–201538) and passive microwave (vegetation optical depth (VOD), 1992–201739) data; these data measure changes in vegetation with different methods and thus provide complementary information for our analysis. We then show that for VOD, the empirically estimated recovery rates from large external perturbations are indeed closely related to the continuously measurable response to small natural fluctuations, quantified here by AC1 and variance. We further show that the AC1 and variance of NDVI are not well matched to empirically estimated recovery rates from large disturbances and conclude that VOD is a more suitable basis for measuring vegetation resilience. We emphasize that while both AC1 and variance have previously been used to estimate vegetation resilience1, their theoretically expected relationships with recovery rates from perturbations, and thus with resilience, have yet to be confirmed empirically for vegetation systems. Moreover, temporal changes in AC1 and variance of remotely sensed vegetation indices, as we investigate here, have rarely been studied40,41.

By comparing with the empirical rates of recovery from external perturbations, we demonstrate using VOD that both AC1 and variance provide robust, empirically verified global resilience measures. On the basis of this relationship, we further quantify global-scale changes in vegetation resilience since 1992 and find coherent resilience loss across land-cover types that has accelerated in the past two decades.

Quantifying vegetation recovery from external perturbations

Vegetation in the natural world is constantly subject to disturbances that vary greatly in frequency and intensity. Many of these signals are subtle, and identifying minor and short-term disturbances is difficult. Large excursions from the typical vegetation state of an ecosystem can, however, be identified by abrupt transitions in time series of vegetation indices. The empirical local recovery rate can then be estimated after each abrupt negative transition by fitting an exponential function to the time series as it recovers towards its previous state (Fig. 1, see Methods for details).

Fig. 1: Global vegetation data.
figure 1

a, Global long-term mean of VOD39 (1992–2017). b, VOD time series for a given location in the Brazilian Amazon (8.375° S, 50.875° W). Raw time series in black, with deseasoned and detrended time-series residual in blue (see Methods for details). c, Recovery of the exemplary time series to the previous mean state after a rapid transition, with commensurate exponential fits. Rare large disturbances, such as those in 2007 and 2010, can be used to track the recovery of vegetation and assign a recovery rate using an exponential fit. See Extended Data Fig. 1 for a corresponding figure based on the NDVI. Exp., exponential.

Full size image

Both VOD (Fig. 1) and NDVI (Extended Data Fig. 1) are subject to the same types of major external disturbances (for example, droughts or fires) that can rapidly reduce both vegetation density (VOD) and vegetation productivity or greenness (NDVI). It is important to note that while both datasets measure vegetation, the data do not describe the same vegetation parameters and hence do not respond identically to external shocks; this can in some cases mean that the number of detected transitions differs between the two vegetation datasets over the same period. In addition, while vegetation recovery is measurable in both data, the time frame of those recoveries, and hence the fitted exponential function, can be dramatically different for the same perturbation (Fig. 1 and Extended Data Fig. 1). Further discussion of the limitations of the disturbance detection procedure can be found in Methods.

Estimating resilience from intrinsic variability

We find globally well-distributed recovery rates from diverse external shocks (Fig. 2 and Extended Data Fig. 2). Not all landscapes have experienced rapid and drastic changes in vegetation over the satellite measurement period; for such regions, it is impossible to directly measure vegetation resilience in terms of recovery from an external shock. Even in regions where perturbations are relatively frequent, they are too sparsely distributed to allow for an estimation of changes in the recovery rate, and thus resilience changes, through time (Fig. 2a).

Fig. 2: Global distribution of recovery rates.
figure 2

a, Recovery rate (for well-determined exponential fits, R2 > 0.2) for VOD (n = 11,538 perturbations for 10,620 unique locations). b, Theoretical estimate of the recovery rate computed via rAC1 = log[AC1] (Methods) from the AC1 of the detrended and deseasoned VOD time series at each location. c, Theoretical estimate of the recovery rate computed via rVar = σ2/(2〈x2〉) (Methods) from the variance 〈x2〉 of the detrended and deseasoned VOD time series at each location. Bare earth, snow and anthropogenic land covers are excluded from the analysis42 (Methods). Note the sparsity of grid cells where there have been abrupt shocks that can be exploited to estimate the recovery rate (a), as opposed to theoretical measures (b,c) that can be computed for all grid cells with vegetation. Also note the similarity of the spatial patterns in b and c and their resemblance to the spatial pattern shown in a as far as there are values for the recovery rate available. d,e, Relative deviation d of theoretical recovery rate estimated from AC1 (d) and variance (e) (for example, d = (r − rAC1)/r). Clear patterns of over- and underestimation of recovery rate indicate that the theoretical framework does not perform equally in all locations. See Extended Data Fig. 2 for a corresponding figure based on the NDVI.

Full size image

The FDT suggests that the rate of a system’s recovery from large external perturbations is related to the variability (quantified by variance22,27,28) and memory timescale (quantified by AC11,23,24,25,28) of natural fluctuations around the equilibrium16. Theory predicts an exponential relationship between the AC1 and the negative recovery rate r, that is, AC1 = erΔt, and a power-law relationship between the variance of the VOD time series x and the recovery rate r, that is, 〈x2〉 = −σ2/2rΔt, where σ is the standard deviation of the driving noise, r < 0, and we set the time steps to Δt = 1 (see Methods for details). For the set of locations where empirical recovery rates can be estimated (Fig. 2a), both AC1 and variance can be derived directly from the corresponding time series. For areas where it was possible to empirically estimate the recovery rate from large perturbations (Fig. 2a), there is broad spatial agreement with the AC1 (Fig. 2b) and variance (Fig. 2c) estimates (see also zoomed-in maps, Supplementary Figs. 1 and 2). Moreover, the two theoretical recovery rate estimates themselves, which are available for all vegetated grid cells, exhibit similar spatial distributions (compare Fig. 2b,c), especially if the relative order of values is considered (see the rank comparison in Supplementary Fig. 3). Note that the AC1-based estimate for the recovery rate r mostly underestimates the recovery rate (Fig. 2d), especially in parts of North America, central Europe and Southern Africa, while the variance-based estimate for the recovery rate mostly overestimates the recovery rate (Fig. 2e).

To more concisely compare the empirical (Fig. 2a) with the two theoretical recovery estimates (Fig. 2b,c), we compare them on a point-by-point basis (Fig. 3). For the VOD, the expected relationships hold remarkably well; for NDVI, the link between empirical and theoretical resilience metrics is much weaker (see Extended Data Fig. 3).

Fig. 3: Empirical confirmation of recovery rates.
figure 3

Comparison between empirically measured recovery rates and theoretical resilience metrics calculated over the five years preceding each transition, for VOD data39. a, AC1 versus recovery rates r from exponential fits to recovering time series with R2 > 0.3; the magenta (blue) line shows binned medians (means), which are close to the exponential fit of the empirical relationship between recovery rate and AC1 values (red line). Grey shading shows data interquartile range. The AC1 thus shows the expected exponential relationship with the recovery rate, but quantitatively, some deviations from the theoretically expected AC1 = er (black line) are apparent. b, Same as a but for the variance. The variance indeed shows the expected power-law relationship with the recovery rate, but as for the AC1, there are some deviations from the theoretically expected \(\langle {x}^{2}\rangle =-{\sigma }_{r}^{2}/2r\) relationship (black line), where we use the spatial mean of the driving noise σr. The mean variance and corresponding interquartile range are also shown for the case where the individual σr values for each grid cell are used to compute the variance (orange line, with shaded interquartile range). c, Binned medians of AC1 as a function of the empirically measured recovery rate r, for increasing thresholds on R2 of the exponential fit to the recovering time series after abrupt transitions, as indicated in the legend. d, Same as c but for the variance. Note that the match between empirical and theoretical estimates of the recovery rate improves the more restrictive the empirical estimation of the recovery rate is; low R2 variance medians in d plot on top of each other until R2 > 0.3. Bare earth, snow and anthropogenic land covers are excluded from the analysis42 (Methods). See Extended Data Fig. 3 for a corresponding figure based on the NDVI and Extended Data Fig. 4 for a corresponding figure using the whole time series to compute AC1 and variance with VOD. See Supplementary Fig. 4 for alternative measures of theoretical variance based on different σ estimates.

Full size image

When considering the AC1 and variance values directly as functions of the recovery rates for all available grid cells together, the theoretically expected relationships are overall corroborated by the observational data, although differences between geographical regions are neglected when investigating the relationship in this way. As expected, some differences are therefore visible (compare Fig. 2). We note that the correspondence between theoretical and empirical estimates becomes substantially better if only recovery rates from exponential fits with R2 > 0.5 are considered, compared with recovery rates from all fits with R2 > 0.1 (Fig. 3). This indicates that the poor exponential fits to the recovering time series after transitions are a key reason for the differences between measurement and theory and suggests in turn that the more reliable the recovery rate estimate, the closer is the match between empirical and theoretical estimates of the recovery rates. We also note that for the variance, uncertainties in estimating the standard deviation of the driving noise σ also probably play a role. Estimating the variance from the empirically determined recovery rate via 〈x2〉 = −σ2/2rΔt requires an estimate of σ. We calculate each individual σi and then bin the resulting data points to obtain the orange curve in Fig. 3b, while we use the globally averaged σ to obtain the black curve.

Global shifts in vegetation resilience

Rapid large-scale perturbations are not evenly distributed in space and time (compare Fig. 2), which renders a reliable estimation of temporal resilience changes in terms of empirical recovery rates impossible. As justified by the relationship between recovery rates and theoretical resilience metrics (compare Fig. 3), we instead calculate resilience in terms of both the AC1 and variance in rolling five-year windows over all vegetated areas (Fig. 4). In the following, we define resilience loss (gain) if at least one of the two indicators (AC1 or variance) shows a statistically significantly increasing (decreasing) trend while the other indicator does not exhibit a significant trend in the other direction (Fig. 4).

Fig. 4: Global resilience trends.
figure 4

a–c, Direction (+/–) of global resilience trends for AC1 and variance using VOD data39 for 1992–2017 (a), 1992–2004 (b) and 2004–2017 (c). Bare earth, snow and anthropogenic land covers are excluded from the analysis42 (white areas) (Methods). Linear trends are calculated on the basis of five-year rolling-window AC1 and variance estimates; only trends with P < 0.05 in either AC1 or variance are shown in colour (see Methods for details on significance testing). Pixels with mixed significant trends (for example, AC1 positive, variance negative) are shown in grey. See Supplementary Table 1 for raw pixel counts. See Extended Data Fig. 5 for global AC1 and variance trends for all three periods, and Supplementary Fig. 5 for latitude-aggregated trends. Note the increases in the strength of resilience loss since the 2000s, especially in the tropics (Extended Data Figs. 5–7).

Full size image

Over the period 1992–2017, the spatial pattern of resilience trends in terms of AC1 and variance is complex (Fig. 4a and Extended Data Fig. 5) but follows consistent latitudinal patterns where equatorial (for example, Amazon and Congo basins) and monsoon-driven (for example, southeast Asia) areas show generally increasing resilience (negative trends in both indicators), and high-latitude areas typically show decreasing resilience trends, especially for the Northern Hemisphere. The global picture for long-term resilience trends is thus mixed; there is only a slight majority of grid cells with resilience losses (54.2%) compared with the number of grid cells with resilience gains (41.6%) over the whole period 1992–2017. The given percentages refer to the set of grid cells that have at least one statistically significant trend in either of the two indicators; the unconfined class with significant yet opposing trends contributes the remaining ~4%. When we restrict the analysis to the first half of our study period (1992–2004), trends are again mixed, with increasing resilience in the tropics and decreasing resilience at higher latitudes (Fig. 4b); these trends are stronger for variance than for AC1 (Extended Data Fig. 5).

From the early 2000s onward, however, we observe a marked increase in resilience loss in terms of both indicators (that is, significantly positive trends in AC1 and variance; Fig. 4c and Extended Data Figs. 5–7). We observe an increase from 28.2% to 59.4% of pixels with resilience loss between the periods 1992–2004 and 2004–2017; the percentage of pixels showing resilience gains decreased from 37.9% to 33.8%. Areas with significant yet opposing trends contribute the remaining 33.8% and 6.8%, respectively; many regions with opposing significant trends until 2004 show coherent resilience loss in both indicators for the period since 2000, 2002 or 2004 (Fig. 4c and Extended Data Fig. 7). Some regions, such as the high northern latitudes, southern Africa and parts of Australia, show consistent resilience losses throughout the study period, which broadly agrees with previous findings based on alternative resilience metrics and AVHRR NDVI data41. Many regions, in particular the equatorial rainforest belt, have reversed from gaining resilience (blue regions, 1992–2004) to losing resilience (orange and red regions, 2000s onwards). Long-term (1992–2017) trends thus conceal a strong reversal from gains to losses in resilience in many regions.

When changes in AC1 and variance are aggregated by land cover42, we infer that evergreen broadleaf forests show overall lower AC1 and variance (higher resilience) than other land-cover types (Extended Data Fig. 6); nevertheless, the global tendency is towards aggregate decreases in resilience (in terms of AC1) across all land-cover classes. These trends maintain a similar form if a three-year or seven-year rolling window is used to calculate continuous changes in resilience (Extended Data Fig. 6). It should be noted, however, that this approach conceals considerable spatial trend variability (Extended Data Fig. 5) and will (although confined to single land-cover types) smooth over vast differences in biomes worldwide; hence, these aggregated time series should be carefully interpreted in the context of the global trend maps (Fig. 4 and Extended Data Fig. 5). Variance presents a more mixed picture when aggregated by land-cover class, with losses of resilience being expressed more strongly since the early 2000s (Extended Data Fig. 6). Previous work proposed that AC1 will always increase towards a critical transition, but variance can in some cases decrease27; the two metrics are also not guaranteed to change at the same rate. This is also to some degree expressed in our global trends (Extended Data Fig. 5), where variance trends, particularly for the tropics, are more strongly negative than for AC1 over both the whole period 1992–2017 and the early period 1992–2004. Both AC1 and variance trends, however, are majority positive for the recent period ~2000–2017 (Fig. 4c and Extended Data Figs. 5–7).

Note that many regions where we observe strong vegetation resilience loss are also fire prone (for example, Siberia, Canada and western North America); increasing fire frequencies due to drier conditions in these regions could explain some of the observed recent vegetation resilience loss43. Increases in temperature, alongside changes in precipitation and weather extremes, could also be a potential driver of changing vegetation resilience; we emphasize, however, that a detailed analysis of the different potential causes for the inferred resilience loss (and in particular its acceleration during the past two decades) is still lacking and is an important topic for future research.

Discussion

Our results provide empirical evidence that both AC1 and variance are directly related to vegetation resilience, defined as the recovery rate from external perturbations. The AC1 and variance can hence be used to estimate resilience in situations where controlled experiments are not possible and external perturbations are rare. Our findings, therefore, justify the usage of AC1 and variance as vegetation resilience metrics1,41,44 and provide an empirical basis for future studies based on these theoretical resilience metrics. However, our results also show that the resilience estimates derived from the common AC1 and variance metrics directly using theoretical relationships may be slightly biased, and instead the modified empirical relationships revealed in Fig. 3 should be used to translate AC1 and variance into the recovery rate as a measure of resilience. On the basis of the thus empirically confirmed relationship between AC1/variance and vegetation resilience, we infer a heterogeneous spatial pattern of resilience gains and losses; resilience losses in the high northern latitudes are consistent since the early 1990s, but in the tropics, we detect gains during the 1990s and pronounced resilience losses since around the year 2000. While the directions of AC1 and variance trends broadly agree (Fig. 4 and Extended Data Figs. 5–7), there remains considerable spatial heterogeneity.

We find marked differences in our results when using the NDVI instead of the VOD data. While we cannot say with complete certainty what drives this disparity, it is likely that differences in the parameters measured by the satellites play a critical role. VOD is primarily sensitive to vegetation density and, thus, will respond to changes in both leafy and woody biomass39. NDVI, however, is sensitive to ‘greenness’, which is often interpreted as vegetation productivity or chlorophyll content; it is well known that NDVI is a poor estimator of biomass45. Recovery in NDVI after a disturbance can thus be rapid, even if a completely new species mix accounts for the post-disturbance vegetation growth (for example, forest replaced by grass). VOD, however, will remain suppressed until vegetation density (for example, leaves and stems) returns. It is thus likely that the empirically derived recovery rates for NDVI contain much higher levels of noise and that some recoveries to previous NDVI values represent a transition to a new vegetation mix rather than a return to the actual previous vegetation state. The relatively poor constraint on vegetation type provided by NDVI is a major barrier to its use in assessing ecosystem state and stability; we therefore propose to rather employ VOD data for such purposes.

A few potential caveats should be kept in mind when interpreting our results. (1) We do not have a strong constraint on the type and cause of the vegetation perturbations used to calculate recovery rates. Sufficient data on all types of disturbances, their spatial extent and their magnitude do not exist; we thus rely on a data-driven approach to estimate the timing and magnitude of a given disturbance. We note, however, that since we determine the empirical recovery rates using only parts of the time series following an abrupt transition, we can estimate a recovery rate without knowing what kind of event (for example, fire, drought) caused the abrupt transition. (2) Possible spurious or missed time-series transitions are carried forward into our analysis of the global relationship between empirical and theoretical vegetation resilience; this probably accounts for some of the scatter seen in Fig. 3 (see Methods for further details). (3) Some changes in variance and autocorrelation are not necessarily related directly to vegetation resilience, for example, in the case of time-lagged vegetation response to water deficits46 that could modify the measured AC1. At the global scale of our analysis, however, we posit that our empirical confirmation of resilience metrics and long-term trends remain robust. (4) We are limited by the mathematical framework to studying only systems that return to the previous state and therefore probably miss many important ecosystem transitions from which there has been no recovery to the original state. Finally, (5) it is important to note that we cannot say for certain whether the acceleration of resilience loss observed in the past decades (Fig. 4) will continue into the future; indeed, it is possible that global vegetation resilience is responding to a (multi-)decadal climate variability mode (compare Extended Data Fig. 6), which could in principle drive a global-scale reversal towards renewed resilience gains. Theoretically, a critical transition will occur when the AC1 reaches a value of one, corresponding to a vanishing recovery rate; in practice, however, extrapolating AC1 trends into the future is not feasible. Our results are based on empirical data and are thus not predictive; they show only how vegetation resilience has changed in recent decades. We have also not assessed changes in the magnitude or frequency of external disturbances (for example, droughts47), which also play a key role in controlling global vegetation health; a comparison between vegetation resilience and contemporaneous changes in external disturbances would provide key context for the attribution of observed resilience changes to explicit drivers. Despite these caveats, our work represents the first empirical confirmation of previously proposed vegetation resilience metrics in terms of variance and AC1 and thus provides the basis for further investigations.

Our study shows that the satellite-derived VOD data can be used to establish a global empirical manifestation of the FDT for vegetated ecosystems. Vegetation resilience, defined as the capacity to recover from external perturbations, can hence be approximated from the characteristics of natural internal variability in terms of AC1 and variance. On the basis of this correspondence, we identify a global loss in vegetation resilience over the course of the past decades, although the spatial pattern is heterogeneous and the inferred resilience changes depend on climate zones. The spatial pattern is complex for the full period for which reliable VOD data are available (1992–2017), with overall resilience gains in the tropical belts and losses in the higher northern and southern latitudes. From the 2000s onwards, however, we find globally almost coherent resilience loss; further work is required to constrain the causes of this loss and especially to investigate whether the observed resilience losses can be attributed to anthropogenic climate and land-use change. Our results establish a firm basis for a global, satellite-driven monitoring system of ecosystem resilience.

Methods

Data preparation

We use two vegetation datasets in our analysis to provide a holistic view of vegetation response to shocks and stresses. (1) VOD at 0.25° spatial resolution; specifically, we employ the Ku-band and use daily values for the period 1992–201739. Note that we do not use the entire VOD data record (1987–) as some pixels exhibit extreme discontinuities before 1992 (Extended Data Fig. 6). We posit that this is due to the change from Special Sensor Microwave/Imager satellite F08 to F11 in the VOD dataset39. While we observe these discontinuities only in the tropics, we choose to discard all data before 1992 for consistency; it should be noted, however, that our global-scale results are robust whether we use 1987 or 1992 as our first year of data (Extended Data Fig. 6). (2) NDVI (from AVHRR) at 1/12° spatial and 15-day temporal resolution for the period 1981–2015; specifically, we use GIMMSv3g38. We further use the moderate-resolution imaging spectroradiometer (MODIS) MCD12C1 land-cover database (2014 annual composite, resampled via the mode of land covers in each VOD/NDVI pixel)42 to break our analyses into distinct land-cover types (for example, Extended Data Fig. 6).

To limit the impact of anthropogenic land use on our results, we further use MODIS MCD12Q1 (500 m, annually 2001–2017) land-cover data to identify any pixels that were at any point during the period 2001–2017 subject to human land use (for example, urban, cropland). We then remove any NDVI/VOD pixels that had one or more anthropogenic land-cover pixels (at least one 500 m pixel) in at least one year between 2001 and 2017. This step helps to remove pixels that, for example, were once logged and then returned to grasslands; those pixels would not be classified as ‘anthropogenic’ for the entire period following the logging and thus might introduce spurious results. While this does not completely eliminate anthropogenic influence from our results (we do not have sufficient land-cover data before the MODIS sensing period), it conservatively removes all 0.25° (~25 × 25 km) regions where human use occurred. We thus cannot completely rule out the influence of human-driven land-cover change on our results at the global scale but have endeavoured to remove it to the furthest extent possible given data limitations. As a final robustness check, we have also used the ref. 48 global deforestation dataset to remove any pixels from our long-term trend data (Fig. 4) that suffered forest loss (Supplementary Figs. 6 and 7); as this dataset also includes non-anthropogenic forest loss—for example, due to natural fires—it serves as an even more conservative land-cover removal step. Removing these additional pixels does not substantially impact our reported long-term trend results or our inferred conclusions.

Cloud cover and other data artefacts are removed from the NDVI data using an upward-smoothing approach to gap filling49. VOD data are resampled to a twice-monthly time step to match the temporal resolution of the NDVI data by taking the median of each time window; this step ensures that divergent results between the two vegetation datasets are not due to spatial or temporal sampling differences. Using these cleaned and evenly sampled time series, we then deseason and detrend the data using seasonal trend decomposition by loess (STL50,51,52). We decompose the full-year signal using a period of 24 (one year at bi-monthly time sampling) and an adaptive loess filter. We use a value of 47 for the trend smoother (one point less than two years) and 25 for the low-pass filter (one point more than one year), according to the rules of thumb originally presented by ref. 50 (see code archive53 for details). We then maintain the residual (deseasoned and detrended) time-series term for analysis.

Note that the VOD dataset is a multi-satellite composite, with variable overlap between different input Ku-band datasets39. As multiple datasets are averaged in different configurations throughout the VOD period, there is the potential for changes in noise levels that could influence the computed AC1 and variance values if the underlying signal (for example, vegetation) changes on a slower timescale than the measurement noise. Stronger averaging associated with an increasing number of satellites would lead to step-wise increases in AC1 and step-wise decreases in variance.

For the period we consider, we do not see step changes in AC1 or variance as would be expected if the noise level or character changed with the introduction or removal of a new satellite; indeed, we see consistent resilience loss during long periods of constant satellite configurations (for example, 2002–2009, Extended Data Fig. 6). Furthermore, there are no contemporaneous jumps in the variance, which would also be expected to change with shifts in data averaging. We posit that the changes in AC1 and variance that we observe are highly unlikely to be driven by data aggregation and are instead representative of a global change in vegetation resilience.

Perturbation detection and recovery analysis

We use two methods to detect perturbations in our residual time series: (1) a moving-average54 and (2) a linear-fit approach55. For both methods, we use an 18-point (9 month) moving window over our residual time series and calculate either the simple mean difference between the first and second halves of our moving window (method 1) or a linear trend over the moving window (method 2). We then smooth these resultant derivative time series with a Savitzky–Golay filter (7 points, first order) to remove high-frequency noise56. Finally, we isolate any derivative values above the 99th percentile and label consecutive time steps as individual disturbance periods. We then use the highest peak within each disturbance as the perturbation date. Note that the results of our analysis are nearly identical whether we use method (1) or method (2) to detect perturbations; thus, we present here only data based on method (1). In our tests, a comparable set of disturbances was found using 12-, 24-, 36- and 72-point moving windows, which resulted in similar spatial (for example, Fig. 2) and global (for example, Fig. 3) patterns; for simplicity, we present only results using the 18-point moving window here.

As we use a percentile approach to delineate large perturbations, we will not always capture each perturbation for a given time series; our detected perturbations will be biased towards the largest excursion within each individual time series. We acknowledge that not all events will be equally represented in both the VOD and NDVI datasets; in the case where a much stronger response is engendered in one dataset than the other, the percentile threshold may not identify the same event in both time series. Furthermore, we will by construction detect some non-significant perturbations, in particular for the case where a given time series does not experience a strong disturbance. We thus impose the condition that the raw time series must descend more than 0.01 to be considered a valid perturbation. While we do not identify every perturbation over the entirety of both datasets, we generate a large and diverse set of recovery rates that are well distributed in space and time. To ensure that our estimated recovery rates represent a return to the previous state, and not a transition to a new vegetation regime, we further apply the condition that the five years of data before and after the disturbance must pass a two-sample Kolmogorov–Smirnov test (P < 0.05). We choose five years as our baseline to minimize the impact of long-term (for example, decadal) changes in vegetation state while maintaining enough data on both sides of the transition for a robust comparison.

For each detected time-series perturbation, we then find the local minimum of the residual time series with a two-month constraint to account for the fact that disturbances are often detected before the residual time series reaches its lowest point. We then take a period of five years after the local minimum and fit an exponential function, capturing both the exponent r and the coefficient of determination R2. To create the map for Fig. 2, if there is more than one transition at a given pixel location, we use the average recovery rate of all transitions. For Fig. 3, we maintain all recovery rates (for example, a single time series could contribute more than one recovery rate). We note that most locations studied have only one significant transition during the study period, and it is a relatively small number that have two or more. The computed transition points and recovery rates can be found in our data repository53.

Resilience estimation

Resilience is defined as the capacity to recover from external perturbations11,12. Quantitatively, it can be determined in terms of the recovery rate r after a perturbation to some value x0:

$$x(t)\approx {x}_{0}{e}^{rt}$$

where x(t) is the state of the system at time t after the perturbation. If r is negative, the system will recover to its equilibrium state at rate ∣r∣. The characteristic recovery time is given by ∣r∣−1. Note that for positive r, the initial perturbation would instead be amplified, indicating that the system is not resilient. Empirically, we estimate r for each perturbation in each residual NDVI and VOD time series as described in the previous section.

The AC1, a measure of how strongly correlated neighbouring time spans of a time series are, has been suggested as a measure for resilience1,23,24,25,57 and more generally as an early-warning indicator for forthcoming critical transitions28,31. Theoretically, this can be motivated from a linearization of the stochastic dynamics around a given equilibrium point x*. For the fluctuations \(\bar{x}=x-{x}^{* }\)

$$\frac{{\mathrm{d}}\bar{x}}{{\mathrm{d}}t}=\kappa \bar{x}+\sigma \eta \,,$$

which defines an Ornstein–Uhlenbeck process with linear damping rate κ < 0 and white-noise forcing with standard deviation σ > 0. It can be shown that the variance \(\langle {\bar{x}}^{2}\rangle\) and lag-n autocorrelation α(n) of the stochastic process obtained from a discretization of the Ornstein–Uhlenbeck process into time steps Δt are given by58

$$\langle {\bar{x}}^{2}\rangle =\frac{{\sigma }^{2}}{1-{e}^{2\kappa {{\Delta }}t}}\approx -\frac{{\sigma }^{2}}{2\kappa {{\Delta }}t}$$

and

$$\alpha (n)={e}^{n\kappa {{\Delta }}t}\,.$$

If the stability of an equilibrium state gradually decreases, κ will approach zero from below, and correspondingly, the variance \(\langle {\bar{x}}^{2}\rangle\) will diverge to positive infinity and the AC1 α(1) will increase towards + 1. These increases in the damping rate κ, as well as the variance of the fluctuations \(\langle {\bar{x}}^{2}\rangle\) and the AC1 α(1) can thus serve as precursor signals for a forthcoming critical transition and, in relative terms, as measures for stability or resilience changes.

The theoretical estimates for the recovery rates shown in Fig. 2b for AC1 and in Fig. 2c for the variance are given in terms of the damping rate κ, obtained by inverting the preceding equations. For the variance, an estimate of the driving noise σr is also needed, which we obtain from

$$\frac{{\mathrm{d}}\bar{x}}{{\mathrm{d}}t}=r\bar{x}+{\sigma }_{r}\eta \,,$$

where we used the empirically estimated recovery rate r rather than the damping rate κ on the right-hand side. Practically, we obtain very similar theoretical expressions for the variance when computed using the σκ obtained when putting κ instead of r into the preceding equation (Supplementary Fig. 4).

For an empirical confirmation of the FDT, we thus have to show that for the observed vegetation data, an exponential relationship between the AC1 and the recovery rate r, as well as a power-law (1/r) relationship between the variance and the recovery rate r, hold. It is important to note that for the comparison between empirical recovery rates and AC1 or variance, we consider only time series that eventually return to their pre-disturbance state, implying that the residual time series under study are, apart from infrequent large perturbations, approximately stationary.

Long-term trend estimation

To better understand temporal changes in vegetation resilience, we calculate the AC1 and variance on moving windows (with a size of 3, 5 and 7 years) over each entire residual time series. Using these windowed AC1 and variance measurements, we calculate Kendall–Tau59 statistics to check for increasing or decreasing trends. As our rolling-window data are by construction serially correlated, we test for statistical significance based on a set of 10,000 phase-shuffled surrogates, which preserve the variance and autocorrelation function of the original time series31,32,33. These phase surrogates are obtained by computing the Fourier transform of the original time series, uniformly randomly shuffling their phases and then applying an inverse Fourier transform to each of them. We then calculate the probability that our measured AC1 Kendall–Tau trends are significant using a threshold of P < 0.05. Finally, we discard six months of data at either end of each time series before calculating trends, as the variance and autocorrelation of the residual produced by the STL procedure are less reliable within one half of the length of the seasonal decomposition window. The python codes to replicate our trend estimation procedure can be found in our code repository53.