Ecosystems are showing symptoms of resilience loss

Ecosystems are prone to non-linear dynamics that can shift their function and structure from one configuration to another [1, 2]. Examples of such regime shifts include the transitions from forest to savannas [3], the collapse of coral reefs [4], kelp forest to urchin barrens [5], peatland transitions [6], or the emergence of hypoxic dead zones in coastal systems [7, 8]. Over 30 different types of regime shifts at the ecosystem scale have been reported in the literature [9]. Because regime shifts can be hard to reverse [10] and can impact the services that society gets from ecosystems [11], predicting regime shifts is a policy relevant endeavor. Even imperfect predictions can help society to prepare for navigating some of these transitions, or avoid them when possible. Yet predicting when and where they will occur remains challenging for most ecosystems [12]. Understanding this is deeply related to our ability to observe and measure resilience.

Resilience is the ability of a system to withstand disturbances without losing its function, structure, and hence its identity [13, 14]. Formally, it is the size of the system's basin of attraction [13, 15, 16]. Several metrics have been used to approximate resilience including the depth of the basin, slope, distance to the threshold, probability of tipping, resistance, elasticity, among others [16]; yet the simplest and cross-system indicators used are based on recovery time [13, 17, 18]. Complex systems when close to critical transitions leave statistical signatures in the time series of its observables known as critical slowing down [1, 19, 20]. It means that the system takes longer to recover after a small disturbance, which translates into increases in variance, autocorrelation, and skewness or flickering [1, 12]. Similar indicators exist for spatial data which includes spatial correlations, discrete Fourier transforms, spatial variance, skewness, power spectrums, and patch-size distributions [21]. These methods however have some limitations. They require long time series to detect useful signals [12], and they can fail when regime shifts are driven by stochastic processes [10, 22].

Recent theoretical and empirical developments have addressed some of these limitations. On the theoretical front, critical speeding up has been proposed as a suitable alternative to detect stochastically driven critical transitions [23]. While critical slowing down relies on the assumption that resilience loss is driven by a widening and depth loss of the current basin of attraction, critical speeding up assumes that the basin shrinks by narrowing the basin [23] (figure 1). Both techniques pick up resilience loss by measuring changes in the higher moments of the time series distribution, detecting sudden increase (decrease) of variance, autocorrelation, skewness or kurtosis. Another proposed proxy of resilience is the fractal dimension [24, 25], which is an indication of self-similarity across scales. The fractal dimension is related to how adaptable a system is to perturbations, or how easily it finds modes to deal with disturbances, applications of this are found in the diagnosis of cardiac disorders [24, 25] and engineering [26]. Exit time has also been proposed as a resilience indicator, it does however requires high resolution time series with multiple shifts to render useful insights [18]. Autoregressive state-space models and dynamic linear models are also proposed as alternatives to generic early warnings which could help circumvent some of their limitations by calculating changes on the leading eigenvalue in time series, if λ > 1 the time series signals a loss of stability and resilience [27–29]. On the empirical realm, recent studies have pointed to remote sensing products [30] and climate change simulations [31] as suitable high dimensional datasets for testing some of these tools in quantifying resilience.

Zoom In Zoom Out Reset image size

This paper aims to identify where regime shifts are likely to occur by detecting signals of resilience loss in terrestrial and marine ecosystems primary productivity. Primary productivity is a proxy of how much living organisms are able to transform solar energy and nutrients into biomass through photosynthesis. Seeing from space, primary productivity is the pulse of life in the biosphere and is typically measured as carbon stored in plants and phytoplankton. This paper applies the traditional early warning signals based on critical slowing down; and adapts the methods to include critical speeding up metrics, and fractal dimension (see section 3). To that end, gross primary productivity, terrestrial ecosystem respiration, and chlorophyll-a concentration were used as proxies of primary productivity of marine and terrestrial ecosystems. These variables have been harmonized by the Earth System Data Lab, so all data layers share the same time (weekly) and spatial (0.25 degree) resolution [32]. To quantify and compare the change in resilience indicators, the absolute difference between the maximum and minimum values per indicator was used (Δ, figure 2), and a segmented regression was used to detect changes in slope and break points in the time series (section 3). A series of logistic and random forest regressions were used to gain insights into what is driving resilience loss in ecosystems world-wide.

Zoom In Zoom Out Reset image size

The generic resilience indicators do not necessarily align with critical slowing down or speeding up theories. Figure 2 illustrates the analysis for one pixel where all indicators show signals of resilience loss. The time series is first pre-processed to remove confounding factors such as seasonality or long term oscillations, log-transformed to reduce the influence of outliers or shock events, and centered to zero mean and unit variance (section 3). All time series were first-differenced and passed a unit root test that guarantees stationarity; in other words, if a signal is detected it is not the product of residual long-term trends or seasonal variation on the data. The generic indicators are then calculated by using a rolling window with length equivalent to 50% of the data. To detect changes in trends of the indicators we used three approaches. First, Δ captures the difference between the minimum and maximum values of the indicator over time. By itself Δ is not very informative unless one has its distribution for the entire planet to compare against (figure 3). In the pixel described in figure 2 the observed Δ is indeed on the tail of the distribution for all indicators. Second, a linear regression was fitted to each indicator to test whether the slopes were different from zero, but it tends to treat sharp increases or decreases of the resilience indicators as outliers. Hence, the third approach was a segmented regression, enabling the detection of a point in time where there is a significant change on the slope of two linear regressions. The summary of the analysis for all pixels across all datasets shows that only in few cases autocorrelation and variance increase or decrease in tandem, the signatures of critical slowing down or speeding up. In contrast, large values of Δ, interpreted here as symptoms of resilience loss, are detected in several different arrangements, most commonly in the agreement between kurtosis and skewness.

Zoom In Zoom Out Reset image size

Ecosystems world-wide are showing symptoms of resilience loss. The absolute difference in resilience indicators (Δ) emphasizes jumps in the time series and enables comparison with normal variation adjusted to each biome type. Arctic ecosystems such as boreal forests, taiga and tundra show the strongest signals of resilience loss globally (figures 4 and S1 (available online at stacks.iop.org/ERL/17/065013/mmedia)). However, the extremes of the distributions (5% and 95% percentiles) of each resilience proxy reveal that all ecosystems are losing resilience, for some of them up to 30% of their global area using the gross primary productivity or terrestrial ecosystem respiration data sets (figures S2 and S3). Despite data incompleteness for marine ecosystems at high latitudes, some signals of resilience loss are detected in Arctic marine systems and the Southern Ocean (figures 5 and S4). The Easter Indo-Pacific and Tropical Eastern Pacific Oceans are the marine realms with larger areas showing symptoms of resilience loss (figure 5). The high oceans (grey area in figure 5), however, show by far the larger areas affected with hot spots outside the Caribbean basin in the Tropical Atlantic, the Tropical Pacific and southeast Madagascar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Symptoms of resilience loss are coherent in space and time. Although the analysis was done independently for each time series and variable, spatial aggregation and coincidence of break points in time suggest that the signals are not artifacts of the data used (figures S5–S7). On the contrary, it supports the idea that there are edges in the three-dimensional space (longitude, latitude and time) that enclose volumes whose dynamics can shift in tandem [31]. Resilience indicators are remarkably consistent across metrics for marine systems both in space (figure S4) and time (figure S7). There is high agreement between kurtosis and skewness across all data sets (figures S2–S4); they can signal the possible shifting of basin of attraction or dynamic transients [1, 12].

An alternative interpretation for coherence of signals is spatial and temporal autocorrelation, or that they are merely responses to shocks on environmental variables. Potential issues of temporal autocorrelation were dealt with by processing the data to the point that the remaining time series were stationary [33]. Issues of spatial autocorrelation are dealt with (below) by fitting explanatory models and subsampling in space, while introducing fixed effects for geography related variables (biome, latitude, longitude). To further test the robustness of the signals, a subsample of time series of detected places (N = 500) were compared against the same time series but with the order of observations changed (10 permutations per time series, 5000 permutations). When the ordering of observations is lost, the signal of resilience loss disappears, meaning it has e.g. a significantly lower autocorrelation and standard deviation on the absolute value scale (figure 6). A Wilcoxon two-sided test, comparing the distributions of permutations against observed data, suggests that the means are indeed different ($p \ll 0.05$ for all comparisons). Thus, the results are not merely noise induced by outliers (e.g. shocks) because the same time series with the same outliers but different ordering does not show the early warnings to the same extent as the observed data. On the contrary, for some statistics like autocorrelation and standard deviation, their Δ values go back to what is expected to be normal (the global medium range, see i.e. figure 3, or the insets in S2–S4 for all variables).

Zoom In Zoom Out Reset image size

An additional robustness test is comparing the generic indicators against model-based early warnings. The former are statistical trends that do not necessarily account for mechanisms and can be vulnerable to false positives and false negatives. The latter are auto regressive dynamic linear models that enable fitting coefficients that change over time, in particular the approximation of the leading eigenvalue of the system λ [27–29]. The leading eigenvalue is a much stronger indicator that the system is stable or not, if λ > 1 it is a signature that the system is going away from equilibrium, possibly an alternative state or a dynamic transient. For each data set a sample of a hundred pixels was randomly drawn for detected places with at least three different resilience indicators, and compared against places that showed no early warning (figure 7). Model-based indicators showed resilience loss for at least half of the sample of chlorophyll A and a third of the sample of gross primary productivity. For places where λ < 0 the distribution still shows a displacement towards one when compared to places with no signals. For terrestrial ecosystem respiration however, model-based early warnings do not support the patterns found by the generic indicators.

Zoom In Zoom Out Reset image size

If resilience is approximated as the size of the basin of attraction, critical slowing down is sensible to reductions in depth while critical speeding up is sensible to reductions in width (figure 1). However, these theories when applied to ecological problems are typically approached as low dimensional models with perhaps one controlling factor, one driver. Recent experimental evidence suggests that when ecosystems are subject to multiple drivers, early warning indicators can fail or provide contradictory signals [34]. Lack of consistency between signals suggests that there is no one preferred theory at place. But instead, multiple drivers are interacting in pushing ecosystems outside their realm of stability, some of them through increasing stochasticity while others through slow forcing.

To test that hypothesis, I used logistic regressions and random forests to investigate what is driving the detection of early warnings. Detection by at least two proxies of resilience loss was used as a response variable. Explanatory variables for terrestrial ecosystems included time series of temperature, precipitation, burned area, and land cover. For marine systems, the explanatory variables were restricted to sea surface temperature and sea surface salinity. Variables with long and frequent time series (temperature, salinity and precipitation) were filtered with a Fourier transform to test whether the signals are predicted by long term variation, annual cycles, fast oscillations, or the linear (slow) trends. All regressions used a stratified sampling design, balanced on detection, and with fixed effects per biomes or marine realm respectively (see method). The subsampling is necessary to control for two sources of bias. First, it breaks biases induced by spatial autocorrelation and the potential inflation of effects or small standard errors [33]. Second, it reduces bias by a relatively larger amount of pixels were resilience loss is not detected, or by biomes that naturally have larger areas.

The strongest predictors of resilience loss are indeed a combination of slow forcing and stochasticity in environmental variables such as temperature, precipitation or sea surface salinity. For terrestrial biomes, mean temperature, mean precipitation and variability in fast oscillations and annual cycles are the strongest predictors, while marine realms are predicted by mean temperature, mean salinity and their variability at different time scales (figure 8). The logistic regression facilitates a relatively straightforward interpretation [35]. Once geography and biome type has been controlled for, the reminder coefficients indicate what increases the odds of detecting resilience loss. However, this linear approach suffers from correlations in the predictors, namely the different time scales at which the hypothesis needs to be tested, whether it is slow forcing or stochasticity at different time scales what drives the signals. In fact, the predicting accuracy for the logistic regression is poor (area under the receiver operating characteristic curve ROC 0.71, 0.65 and 0.59 for gross primary productivity, terrestrial ecosystem respiration and chlorophyll A respectively). The random forest approach leads to higher predictive power (ROC 0.893, 0.874, 0.827 respectively). It is robust to potential correlations and required less pre-processing for feature engineering but is less amenable to interpretation [35]. It reveals which variables improve predictions but not necessarily in which direction they are affecting resilience loss. The results of both approaches confirm the hypothesis that resilience loss is driven by a combination of slow forcing and stochasticity on potential drivers.

Zoom In Zoom Out Reset image size

Resilience is the ability of any system to deal with change while keeping its structure and functions, thus its identity [13, 14]. The dynamic indicators used here approximate resilience as recovery time [1, 23], and for the fractal dimension a collection of behaviors available to deal with disturbance [24]. These metrics, however, do not directly inform about other dimensions of resilience such as the size of the basin of attraction, the amount of disturbance that the system can stand, the mean return time, the distance to tipping points and thresholds, or adaptive and transformative capacities [14, 16]. The indicators here used can be seen as symptoms of instabilities being developed on the time series of primary productivity for global ecosystems over the past two decades. The symptoms do not imply imminent regime shifts, they help identify places where their probability is increasing given the limitations of the data. The data is relatively short to account for critical transitions in the biosphere, yet it is one of the best observational records to approximate resilience globally. These symptoms are necessary but not sufficient evidence of unfolding regime shifts. Yet, the analysis here presented enables the identification of areas vulnerable to regime shifts at a scale relevant to decision makers.

Globally 29% of terrestrial biomes and 24% of marine realms show signals of resilience loss in at least one of the indicators used (figures 4, 5 and S1). The overall patterns here reported agree with recent reports documenting ecosystems' degradation worldwide. Forests are becoming more vulnerable to droughts, and the combined effects with increasing fire frequency are exposing them to major diebacks expected by mid-century [36]. Temperature thresholds for terrestrial primary productivity have been identified [37, 38] where carbon uptake is potentially degraded (sink to source transition). Less than 10% of the terrestrial biosphere has already crossed the threshold, and under business-as-usual scenarios, half of the biosphere is expected to cross these thresholds by the end of the century, with the most affected areas being Canadian and Russian boreal taigas as well as the Amazon and South East Asia rainforests [37]. These are biomes where the strongest signals of resilience loss were detected with a right shift on the distribution of $|\Delta|$. Other studies quantifying terrestrial ecosystems' resilience with NDVI data also show strong signals from tundra and boreal forests [39]. A recent quantification of aridity thresholds shows that up to 28.6% of current dryland area can cross these thresholds by 2100 in the most drastic climate scenarios [40]. The results here presented confirm early warning signals of resilience loss in drylands as well, particularly with the fractal dimension, skewness and kurtosis (figures S2 and S3).

The marine patterns presented here also align with previous studies. Deepening of the ocean's mixed layer can decrease light conditions near the surface, decreasing nutrient exchange in the water column and consequently primary productivity [41]. The area outside the Caribbean basin reported here coincides with an area where salinity has contributed to ocean stratification [41]. Upwelling systems are also hotspots where resilience loss is identified. Upwelling strength is expected to change with climate change, with strengthening already reported in the California and Benguela currents, while weakening in the Iberian-Canary system [42]. Upwelling weakening can limit nutrients in marine food webs, while strengthening can over enrich nutrients and facilitate the onset of oxygen minimum zones [7]. This study provides an additional line of support that these systems are being destabilized.

The results are limited by the temporal and spatial scale of the available data. If the grain of the data is not frequent enough to match fast dynamics, long enough to capture change in slow processes [1, 12, 18], or the spatial resolution is too coarse, local transitions in space and time can be missed. This is the case with hypoxic areas, where large oxygen minimum zones are identified, but not the diverse range of smaller local cases that have been previously reported [7]. Only regime shifts that are drastic enough to change primary productivity as observed from remote sensing products can be identified. Regime shifts that impact specific populations or community assemblies without changing primary productivity are missed. Such is the case of coral transitions in the Great Barrier Reef where ~50% of the reef community collapsed following the heatwave events of 2016–17 [43].

Because of these limitations, the results likely are an underestimation. The estimates are also conservative, with an arbitrary 5% and 95% quantile of the Δ distribution as detection threshold. A lower cutoff would enlarge the areas where resilience loss is detected, but also increase the risk of false positives. A similar study using NDVI data estimated up to ~65% of terrestrial ecosystems show early warning of critical transitions, with strong bias towards boreal forest and taiga [39]. The estimates presented here complement previous efforts [30, 39] in taking into consideration fixed effects by biome and a pre-processing technique that removes seasonality and long-term variations that can lead to errors, or bias towards high variable environments (higher latitudes).

Another limitation of this study is the lack of ground truthing. The results provide a spatial prediction of where ecosystems might be losing their resilience. But since critical transitions have not happened yet in many of these places, it is not possible to contrast the models presented here with their real predicting power. Current databases that track such transitions are biased to studies and observations on the global north and coastal ecosystems [7, 9], their coverage is not yet sufficient to be used as ground truth across all ecosystems and biomes. Nevertheless, the robustness of the signals detected were tested against null permutation models and model-based early warnings. The first test showed the signals are not an artifact of the data, potential temporal autocorrelation was dealt with stringent data pre-processing and randomizing the ordering of observations. The second test supported our findings for chlorophyll A and gross primary productivity, but failed in supporting the findings of terrestrial ecosystem respiration.

Model-based approaches such as autoregressive time varying models can help interpreting the results and robustness of generic indicators [27–29]. In this study, these methods gave stronger confidence to the results from gross primary productivity and less support to terrestrial ecosystem respiration derived signals, although both results qualitatively point to the same patterns. The combination of methods thus helped distinguish what are optimal observables of resilience loss, an open question and possibly a fruitful avenue for future research [12]. Unfortunately, model-based methods are computationally expensive and do not scale up to the data requirements of remote sensing products, on the order of 10⁵ to 10⁷ spatial pixels (depending on spatial resolution) times time series length. Generic indicators are still useful in narrowing down the scope at which computationally expensive methods become practical.

Computing the probability density function of Δ helps interpreting signals of resilience loss. Previous work using generic indicators typically have ground truth in the form of modeling simulations with some noise or actual experiments. Because the ground truth is known, there was not concern whether the magnitude of the increase (decrease) in early warnings was big enough to actually be considered a warning. Here we do not have ground truth like experiments, but we can compute the distribution of each indicator for the entire planet enabling the comparison of what constitutes a big jump on an indicator as opposed to its normal expected variability for each biome (figure 3). The magnitude depends on the indicator and data pre-processing choices, but the position of the pixel in the distribution is relatively unaffected. Thus, Δ is less sensitive to pre-processing choices. While in theory increases or decreases on a particular statistic can be interpreted as resilience loss (an apparent contradiction), the distribution of Δ helps interpreting whether an increase or decrease is towards the tails of the distribution (resilience loss) as opposed to an increase or decrease towards the center of the distribution (a resilience gain). The fractal dimension here introduced is unambiguous and do not depend on scale, for time series it is bounded between 1 and 2 and higher values are always more resilient [24, 44]. Future studies can use these two facts to track ecological recovery.

Despite its limitations, this first-order approximation to resilience loss can help outline priority areas for management. Russia, Canada, the US, and Australia are countries with the largest areas of resilience loss identified, yet by proportion of territory, small island states top the ranking. When accounting for the diversity of ecosystems showing signs of resilience loss, megadiverse countries like Brazil, India, Mexico, Indonesia, Australia, or Colombia are on the top 10 (figure 9). Australia has recently been reported as a hotspot for ecological collapses for both marine and terrestrial ecosystems [45]. The spatial resolution of the maps here presented enable countries, regions, and even municipalities to update land use planning and take the vulnerability of their ecosystems into consideration. Companies for example, can take such risk into account when deciding on relocation, resource outsourcing, or investments. Countries and municipalities can balance the trade-off between maximizing a particular ecosystem service (e.g. food production) in favor of multifunctional landscapes that include other values such as recreation, spiritual values or conservation. No matter the scale, we all have a role to play in caring for ecosystems resilience and maintaining their ability to provide the ecosystem services we all depend on.

Zoom In Zoom Out Reset image size

Future global resilience assessments could benefit from other data streams, particularly other anthropogenic drivers overlooked here. Answering the question of resilience of what to what will require observations in several ecological and social data streams. For example, longer recovery time might be related to more frequent and intense perturbations (rate induced tipping), not necessarily reductions on the basin of attraction size [46]. Recent qualitative reviews on regime shifts drivers diversity show that no regime shift is caused by single drivers. Multiple drivers often interact in non-random patterns where food production, nutrient cycling, climate change and urbanization co-occur strongly [47, 48]. Additional data would however require long time series coverage and high spatial resolution. The accuracy of the predictive models presented here can be improved by using non-linear approaches such as deep neural networks or other machine learning techniques [49]. However, these approaches have limited interpretability and often require manually annotated data sets to measure performance. Qualitative efforts such as the Global Ocean Oxygen Network [7] or the Regime Shifts Database [9] can provide the annotated examples to train such models. A spatially explicit map of resilience loss can also help quantify the risk of cascading effects in ecosystems previously identified [47]. As new Earth observations become available, these global maps can be updated and track how ecosystems resilience is evolving, where are they recovering and where they are becoming more vulnerable. This paper showcases the first steps toward an ecological resilience observatory.

3.1. Data

3.2. Pre-processing

3.3. Proxies of resilience

3.4. Regressions

All code used in this analysis is available at https://github.com/juanrocha/ESDL. The computations in the paper would have not been possible without open source software including the following packages: broom [57], corrgram [58], deSolve [59], diptest [60], facterize [61], forecast [62], fractaldim [63], furrr [64], future [65], GGally [66], ggridges [67], ggplot2 [68], here [69], janitor [70], MARSS [71], progress [72], patchwork[73], purrr [74], raster [75], rsample [76], segmented [77], sf [78], slider [79], tictoc [80], tidymodels [81], tidyverse [82], tsibble [83], and urca [84] in the R programming language [85], and ESDL, FFTW, NetCDF, and Statistics in Julia [86].

This work would have not been possible without the open data provided by Copernicus Climate Change and Atmosphere Monitoring Service, NASA, FLUXCOM initiative, and the Global fire emissions database. JCR would like to thank the European Space Agency and the Max Planck Institute of Biogeochemistry for an early adopter grant to use the Earth System Data Lab curated data sets and computational facilities. The manuscript has benefited from comments from Fabian Dablander, Steven Lade, Thorsten Blenckner, Megan Meacham, Giesela Rohr and Stephen Carpenter. JCR was supported by Formas Grant Nos. 942-2015-731, 2020-00198 and 2019-02316, the latter through the Belmont Forum.

All data used in this study is publicly available through the Copernicus Climate Change and Atmosphere Monitoring Service, NASA, FLUXCOM initiative, or the Global fire emissions database. Links to each data set are provided when introduced under the section data or regressions. No new data were created or analyzed in this study.