Ecosystem Resilience Monitoring and Early Warning Using Earth Observation Data: Challenges and Outlook

Abstract

As the Earth system is exposed to large anthropogenic interferences, it becomes ever more important to assess the resilience of natural systems, i.e., their ability to recover from natural and human-induced perturbations. Several, often related, measures of resilience have been proposed and applied to modeled and observed data, often by different scientific communities. Focusing on terrestrial ecosystems as a key component of the Earth system, we review methods that can detect large perturbations (temporary excursions from a reference state as well as abrupt shifts to a new reference state) in spatio-temporal datasets, estimate the recovery rate after such perturbations, or assess resilience changes indirectly from stationary time series via indicators of critical slowing down. We present here a sequence of ideal methodological steps in the field of resilience science, and argue how to obtain a consistent and multi-faceted view on ecosystem or climate resilience from Earth observation (EO) data. While EO data offers unique potential to study ecosystem resilience globally at high spatial and temporal scale, we emphasize some important limitations, which are associated with the theoretical assumptions behind diagnostic methods and with the measurement process and pre-processing steps of EO data. The latter class of limitations include gaps in time series, the disparity of scales, and issues arising from aggregating time series from multiple sensors. Based on this assessment, we formulate specific recommendations to the EO community in order to improve the observational basis for ecosystem resilience research.

Article Highlights

• Substantial differences exist in defining vegetation “resilience” and its estimators, creating sometimes unnecessary uncertainty

• Progress has been made in constructing indicators of resilience, but there are underappreciated caveats associated with the measurement process of Earth observations

• Novel developments and increasing diversity in datasets from Earth observations, in combination with theoretical insights, can help to improve resilience monitoring from space

1 Introduction

Concerns have been raised that different terrestrial ecosystems might respond to ongoing anthropogenic climate and land-use change with abrupt and potentially catastrophic regime shifts. Indeed, several vegetation ecosystems, such as the Amazon and other tropical rainforests, but also drylands or boreal forests, have been suggested to be potential Earth system tipping elements (Lenton et al. 2008; Armstrong McKay 2022; Boers et al. 2022). To reliably project the future fate of terrestrial ecosystems and to assess the risk of local- to regional-scale regime transitions it is crucial to be able to measure their resilience based on empirical spatiotemporal data.

The resilience of an ecosystem is often qualitatively defined as the ability of the system to cope with disturbances, where a disturbance is “any relatively discrete event in time that disrupts ecosystem, community, or population structure and changes resource pools, substrate availability, or the physical environment” (White and Pickett, 1985). Disturbances are usually induced by anthropogenic or natural discrete events which act to perturb a system away from its dynamic equilibrium, after which the system gradually recovers. In principle, disturbances do not necessarily have to be extreme or rare, though practically they need to be large and rare enough to be distinguished from the normal variability of the system.

We note that, in this review, “equilibrium” refers to any state of the system that does not show significant changes in its properties over time, except when external forces impose a change. A special case of an equilibrium is a fixed point x* where the evolution operator f of an autonomous dynamical system dx/dt = f(x) becomes 0, i.e., f(x*) = 0 (for example, see Kuznetsov 2004). Under more complicated conditions with time varying forcing, one can define a “pullback attractor” to which the system converges when the initial condition is applied at a time point shifted infinitely into the past (Chekroun et al. 2011). This perspective of equilibrium is purely mathematical and differs from the concept of thermodynamic equilibrium in physical systems, which refers to a state of maximum entropy for a given set of constraints.

From a dynamical systems perspective, a local and a global notion of resilience can be distinguished (Holling 1973; Krakovská et al 2023). The local perspective applies close to an equilibrium state and defines resilience as the recovery rate after a perturbation (e.g., Smith et al. 2022). Such a perturbation could be a single disturbance event. If the perturbation is small, and the reference state is a stable equilibrium in the sense of dynamical systems theory, this definition of resilience is the fixed point’s stability, expressed by the eigenvalue of the slowest mode of recovery (Holling 1973). This equilibrium-based notion of resilience has been coined as “engineering resilience” (Pimm 1984).

In contrast, the global perspective defines resilience by how large a perturbation must be for the system to not recover from it at all. In dynamical systems terms, this non-equilibrium based definition refers to the width of a basin of attraction, or “ecological resilience” (Holling 1973) and is sometimes also called “latitude” (Walker et al. 2004). Figure 1 illustrates these two perspectives of resilience with a widely used analogy, where the system’s state is represented by a ball in a potential landscape.

Fig. 1

Conceptual difference between engineering resilience (local linear stability, return rate) and ecological resilience (basin width). The numbers in the lower panel correspond to the numbers in the upper panel and indicate where the system is in the depicted potential landscape over time

These two views of resilience are not unrelated. A system can only display engineering resilience if it also has ecological resilience regarding a disturbance (Fig. 2, case 1). Environmental change (captured as a slowly changing parameter in models) can reduce the width of the basin and at the same time make it less steep around the equilibrium. As a consequence, the system is more vulnerable to disturbances, first in the sense that it shows a smaller recovery rate, and second in the sense that these disturbances can tip the system to an alternative state (Fig. 2, case 2). In the extreme case, the whole basin and with it the stable equilibrium disappears at a so-called bifurcation point (Fig. 2, case 3). At this point, both ecological resilience (basin width) and “engineering” resilience (recovery rate) would decrease to zero. Consequently, indicators measuring “engineering resilience”, which are closely related to the concept of critical slowing down (CSD), are often referred to as “early warning signals” of abrupt shifts associated with a loss of ecological resilience (Scheffer et al. 2009). It is, however, a common misconception to believe that resilience monitoring (in the sense of measuring changes in engineering resilience) would only be informative if alternative equilibria existed; even when alternative equilibria do not exist, the engineering resilience measures the recovery capacities of the system. This also means that CSD indicators alone do not provide any informative value on the existence of alternative equilibria (Kéfi et al. 2013).

Fig. 2

Idealized cases of abrupt events and how they are related to different kinds of resilience, and resilience indicators. case I (blue): A large external state perturbation happens, and thereafter the system recovers gradually. The recovery time measures “engineering resilience”. case II (yellow): A similar event occurs but this time the external shock is larger (or the resistance of the system smaller), so that the system does not show “ecological resilience”, i.e.,, it does not return to the previous state. case III (red): A gradual parameter change drives rapid change in the system because it crosses a bifurcation point where the old equilibrium disappears

Another quantifier of resilience which to some degree integrates the above two views can be defined based on mean first exit times from a given basin of attraction to the basin of attraction of an alternative equilibrium (Ciemer et al. 2019; Arani et al. 2021). A given state is naturally more resilient in this sense if the expected first passage time to the alternative basin is longer. To first approximation, the resilience in terms of mean first exit times depends on the width and depth of the basin of attraction, as well as on the strength of the noise. This quantifier of resilience can be called global in distinction to the local concept of engineering resilience. For an excellent in-depth review of different concepts of ecosystem resilience we refer to Dakos and Kéfi (2022).

In this review, we focus on terrestrial ecosystems and how to measure changes in their resilience to disturbances under increasing anthropogenic pressures, as an important field of application for Earth observation (EO) systems. The motivation for this article is to work out major conceptual as well as practical caveats that need to be addressed before the value of novel and better resolved EO data can be manifested in monitoring resilience (and changes thereof) in the Earth system. While recent studies have argued for enhanced and improved resilience monitoring, highlighting specific examples (Lenton et al. 2022; Lenton et al. 2024), we here discuss in more detail regarding the conceptual, practical and computational challenges of data analysis, in light of the most recent literature.

Some of the discussed caveats are not exclusive to EO data, but all of them are important to consider in EO-based resilience monitoring. Our intention is (i) to make “data users” (e.g., Earth system modelers and dynamical systems experts) aware of these challenges and of arising opportunities, and (ii) to make “data producers” (e.g., remote sensing experts) more aware of the specific requirements that resilience monitoring implies for the data. We believe that both would contribute to more meaningful resilience assessments by future studies. The majority of work discussed in this article focuses on engineering resilience and the closely related concept of CSD; first, because there are arguably more observed examples of disturbance and recovery than examples of crashes without recovery, and second because it is less clear how to accurately define/measure ecological resilience mathematically and in practice.

We will first (Sect. 2) discuss the conceptual challenges associated with resilience quantification. After that, we focus on the practical problems related to estimating ecosystem resilience from EO data. To this end, we review how engineering resilience has been measured in practice using two different approaches.

The first approach (Sect. 3) is by computing properties of the recovery process after single large perturbations. This includes methods that detect disturbances in time series (Sect. 3.1) and methods that quantify recovery rates after such disturbances (Sect. 3.2). This can be considered the direct way to estimate the recovery rate and thus engineering resilience.

The second approach (Sect. 4) is to infer engineering resilience indirectly, from the normal variability of the system as it is permanently adjusting to many small perturbations. The fact that the recovery rate can be inferred from natural variability is suggested by theoretical insights, such as those from the theory of stochastic dynamical systems (Kuehn 2011) and the fluctuation–dissipation theorem (Kubo 1966). These theoretical relationships have only recently been empirically validated in vegetation ecosystems (Smith et al. 2022; Smith and Boers 2023a, 2023b). The major advantage of this second approach is that it allows to measure resilience everywhere and not only where large perturbations have occurred, and continuously in time. We discuss metrics to measure resilience indicators from permanent fluctuations (Sect. 4.1), and ways to cope with long-term trends and seasonality (Sect. 4.2) and spatially and mechanistically complex ecosystems (Sect. 4.3).

In Sect. 5, we discuss future opportunities and challenges for resilience monitoring using EO data. As more and more data with higher and higher resolutions from different sensors is becoming readily available, the vision to monitor vegetation resilience from EO data can be seen as a Big Data challenge. However, this requires both meaningful measures of resilience in high-dimensional systems and a way to automate their computation as well as the detection of significant changes, and there might be limitations associated with EO data. Finally, Sect. 6 provides a brief summary and our most important conclusions.

2 Challenges in Defining Resilience

2.1 Diversity of Resilience Metrics

The concept of resilience has been a central topic in the field of ecology for approximately 50 years (Yi and Jackson 2021). Defining resilience is challenging, not only because of differences in definitions and semantics, but also because of the difficulties in finding appropriate metrics to quantify resilience. While ecological and engineering resilience as outlined in Sect. 1 are the most common categories, other definitions of resilience exist, and many variants are used when quantifying resilience (Ingrisch and Bahn 2018; Yi and Jackson 2021).

The term ‘resilience’ is also used differently in different scientific disciplines and contexts, and even where the qualitative definitions appear to be the same, the mathematically defined metrics can differ. In principle, there are two kinds of challenges in trying to define resilience. The first is a conceptual challenge of how to precisely define resilience and distinguish it from other related properties. Often, different studies focus on different properties of disturbances and recovery processes, but nonetheless use the same terminology (Ingrisch and Bahn 2018; van Meerbeek et al. 2021), which can create unnecessary confusion, adding to the actual epistemic uncertainties. The second kind of challenge arises when such concepts are applied to observations of geoscientific or ecological systems, where practical estimation problems occur (e.g., arising from uncertainties and gaps in the measurement process or biases induced by other effects). We will deal with conceptual challenges in this section, while later sections focus on more practical problems associated with applying the definitions to EO data.

Regarding conceptual challenges, an important factor is that the recovery rate generally varies over time in systems with more than one dimension (e.g., an ecological model describing several coexisting species), even if the linearisation around the fixed point is valid. In the following, we assume that the one-dimensional simplification is justified. A more complete formal definition and comparison of resilience properties in higher dimensions is given in Arnoldi et al. (2016).

Even in a one-dimensional setting, a great deal of ambiguity arises from the fact that there are several properties that characterize and influence the recovery process (as shown in Fig. 2, case 1), and that one can measure: the magnitude of the disturbance, the ratio between this magnitude and the size of the external perturbation (resistance or robustness, also called Lloret’s resilience, see Lloret et al. 2011), the total duration of the recovery process, the ratio between disturbance magnitude and duration (abruptness), the initial pulse rate, the rate of recovery (engineering resilience), the area under the pulse curve (see Jentzsch and White 2019, and Fig. 1 therein), the magnitude of the perturbation needed to change an ecosystem to a different set of processes and structures (Angeler and Allen 2016), as well as combinations of multiple properties (Mitra, Kurths and Donner, 2015). In systems more complex than the simple dynamical systems shown above, it can also happen that the system does recover to some extent, but not fully (i.e.,, an intermediate type between the blue and yellow trajectory in Fig. 2), which is referred to as the degree of recovery (e.g., see Tao et al. 2022). For a comprehensive review of resilience metrics, see Yi and Jackson (2021) and Dakos and Kefi (2022).

Among the mentioned properties of disturbance and recovery phenomena, resistance is arguably among the most popular quantifiers (Pimm 1984; Ingrisch and Bahn 2018). Resistance has been defined as the ability of individuals or ecosystems to persist or stay essentially unchanged under a disturbance (Grimm and Calabrese 2011; Nimmo et al. 2015). The higher the resistance, the smaller is the impact of a given disturbance on ecosystem state variables. Therefore, resistance is typically defined based on the concurrent impact of a disturbance, quantified either by the (relative) difference or ratio of the state variable before and following a disturbance (Yi and Jackson 2021), and is closely related to disturbance severity used as characteristics of disturbance regimes in ecology (Turner 2010).

Resistance has also been included as one component of engineering resilience, next to recovery rate (Holling 1996). However, resistance is a property that is different from the recovery rate from a perturbation (Ingrisch and Bahn 2018) (though in a simple one-dimensional dynamical system setting, both should be proportional), and it can be insightful to consider both when assessing resilience. Particular care has to be taken when controlling for one of them to estimate the other. Even when the same qualitative definition of resilience is used (e.g., recovery rate), it matters how the quantified resilience is normalized in order to obtain unitless quantities. For example, the same absolute recovery time will lead to different recovery rates depending on whether one normalizes the time series by the baseline (reference state) value, or by the impact (amplitude) of the perturbation (Ingrisch and Bahn 2018), which is determined by the system’s resistance.

In the following, if not specified further, we refer to resilience as the engineering resilience in the sense of a recovery rate that is not normalized with the resistance. We will discuss some of the associated conceptual challenges in Sect. 2.2 and the practical problems of its estimation from Earth observations in Sect.  3 and 4.

2.2 Conceptual Challenges

Real ecosystems are much more complex than the idealized dynamical systems shown in Figs. 1 and 2. Hence, a number of additional caveats have to be kept in mind when assessing their resilience.

• (i) Equilibrium assumption.

  The notions of equilibrium (a stable fixed point in terms of dynamical systems) and out-of-equilibrium is questionable in complex ecosystems that might be better described with chaotic dynamics even under constant environmental conditions (Hastings and Wysham 2010). Another caveat is how to define a system’s equilibrium in practice. Examples for a predefined state can be the empirically observed long-term average (Ingrisch and Bahn 2018), but typically covering only a few years given observational constraints, or it can be defined in absolute terms such as the FAO (Food and Agriculture Organization of the United Nations) definition of forests (> 40% tree cover, 5 m height), yet how different definitions impact the outcome of resilience assessments remains largely not well quantified.

• (ii) High dimensionality.

  Even in case one-dimensional concepts apply mathematically (as mentioned above), there is typically a variety of variables to work with. Yi and Jackson (2021) consequently argue that it is challenging to calculate ecosystem resilience by a single universally applicable measure or method. For example, multiple state variables have been used to estimate resilience, some related with a system’s functioning (e.g., tree growth, assessed by basal area increment or tree ring width), others related with a system’s state (e.g., leaf area or vegetation indices related to photosynthetic activity, assessed by optical remote sensing) or a system’s composition (e.g., species richness). These reflect different facets of ecological systems that do not necessarily correspond to the same measures of “resilience”. As an example, when observing a deviation of the normalized Difference Vegetation Index (NDVI), the post-recovery NDVI, even if similar to pre-disturbance NDVI, may represent a different species composition or even another vegetation state (e.g., shrubs, grasslands instead of trees; Buma 2012). It is hence important to be aware of and clearly state the variables and mechanisms one aims to resolve, which also affects the sampling time (see below). For example, the generation time of a tree species may reflect the time scale of recovery for a forest after a severe dieback event, but may not matter for the recovery of leaves or biomass after less extreme events. Furthermore, measurable variables typically correspond to proxies of state variables, rather than the state variables themselves (e.g., spectral vegetation indices as proxies for photosynthetic activity/productivity). Given this complexity, it is generally of advantage to have knowledge on the drivers and main ecosystem processes involved, in order to make meaningful choices on relevant variables to assess the resilience of the system. This also applies to the right time period and sampling intervals, as we will discuss below.

• (iii) Interactions between disturbance and recovery time-scales.

  Ecosystems with low resilience will only slowly regain their pre-disturbance state, or may not ever recover completely. The essential metric used to assess resilience is thus ‘recovery time’ or ‘recovery interval’, which needs to be interpreted in light of the disturbance interval. As an example, a system taking 100 years to recover is likely to fully recover if disturbances only occur every say 500 years, but if same system is disturbed every 50 years, it will never reach its pre-disturbance state (see Senf and Seidl 2021a, b), and its response to later perturbations might change due to nonlinear or even irreversible legacy effects (Hansen et al. 2022; Anderegg et al. 2020).

  This is an important point to consider, as changing disturbance frequencies (i.e.,, due to climate change) might change the resilience of a system even though recovery times stay similar. To judge whether a system is actually at risk of experiencing a permanent state change, we hence need to consider both changing recovery times but also changing disturbance intervals. An example of this is increasing fire frequency observed in some regions, where ecosystems will be stuck in successional stages because they reburn before late-successional stages can be reached (a so-called fire-trap; see Freeman et al. 2017 and Hoffmann et al., 2020). In contrast, changes in ecosystem composition in response to changes in fire frequency can also promote resilience, depending on the ecosystem property under consideration, e.g., increased carbon storage (Mack et al., 2021).

  Additionally, the recovery to a given state depends on changes in the forcing, e.g., drought recovery might be faster if coinciding with increased water availability. An alternative way to assess recovery is to evaluate how ecosystem function is changed following a disturbance. This can be quantified by evaluating departures of the sensitivity of the target variable(s) to forcing drivers relative to a long-term sensitivity (functioning) baseline (Yu et al., 2022; Wu et al. 2022).

• (iv) Memory effects and system boundaries.

  As explained in point (iii), the time scale separation of an abrupt disturbance and the consequential slow recovery is an idealization. Many ecological systems exhibit memory (Ogle et al. 2015; Cranko Page et al. 2022) and specifically stress memory (Walter et al., 2013; Ruehr et al. 2019). Ecological memory might lead to altered responses to perturbations, either by conferring increased or reduced resistance (Gessler et al. 2020; Wu et al. 2022). Increased resistance might result from positive legacy effects of responses to disturbance. For example, investment in root growth during drought (McDowell et al. 2008) might attenuate impacts of a subsequent drought, or reduced competition due to mortality might increase resources available for surviving individuals (Gessler et al. 2020). Reduced resilience following disturbance might result from functional impairment (e.g., damage to the plants’ hydraulic system; see e.g., McDowell et al., 2016; Ruehr et al. 2019; Wagner et al. 2023). Further, some disturbance resistance strategies might increase vulnerability to subsequent disturbances. For example, drought avoidance strategies like stomatal closure or defoliation might increase the risk of carbon starvation and long-term metabolic decline (Poyatos et al. 2013; McDowell et al., 2016), increasing tree vulnerability to biotic attacks. In addition, these effects likely depend on the magnitude of the disturbance impact (Ruehr et al. 2019).

  Hence, one needs to keep in mind that concurrent impacts might not only reflect the direct disturbance effect, but also the legacy effects from similar or other types of past disturbances due to stress memory (Gessler et al. 2020; Bastos et al. 2021; Anderegg et al. 2020). It has therefore been proposed to quantify resistance by the concurrent sensitivity to a given disturbance driver (ΔV/ΔF, with state variable V and forcing F, Fu et al., (2020); Xiao et al. (2023)), which better allows to estimate “the persistence of relationships within a system” (Holling 1973), and allows to compare shorter periods of time (Fu et al. 2020). It should be noted, however, that these are also likely to be time-scale dependent (Linscheid et al. 2021).

  Besides these challenges related to time scales, resistance and resilience can also be hard to disentangle practically because of the presence of climate-ecosystem and ecosystem-biogeochemical feedbacks. For example, the precipitation-evaporation feedback (Eltahir and Bras 1996) leads to the fact that precipitation at a certain location is not an external driver, but already affected by the response of the ecosystem to other drivers, or the same drivers at other locations. Another example is that increased carbon storage under more frequent fires reported by Mack et al. (2021) was linked to higher nutrient availability, hence, the direct effect of fires in isolation (as might be possible to be studied in models) is not well separable in observations. It is hence not always possible to observe the internal ecosystem dynamics and external stressors separately from each other.

• (v) Value judgements and other subjective choices.

  In highly complex (socio-)ecological systems, it also matters whether the restored ecosystems fulfill the same functions (or provide similar ecosystem services) as prior to the perturbation, i.e., if they can adapt to change in a way that is judged as beneficial (Helfgott 2015). In this sense, “resilience” has more facets than just the recovery to a previous state, or the persistence in remaining in a certain regime. It has been argued that it is a feature of resilience for systems to be self-organized (rather than just changing when being forced); and to permanently change in order to adapt to the environment by learning. For instance, species and hence ecosystems permanently change due to evolution. It is then particularly relevant to decide which function or property of the system shall be the focus (“resilience of what to what?”; Carpenter et al. 2001; Cutter 2016), and, in case of socio-ecological systems, “resilience for whom?” (Cutter 2016). If the structure and/or composition changes, but from a human perspective the ecosystem provides the same function, we might still call it resilient. In this sense, there is also a subjectivity that depends on the selected functions (or ecosystem services) and how we value them (e.g., timber production rate, species diversity, habitat for particular species). As different functions can have distinctive resilience, it matters which ones we focus on, which is necessarily a strongly human-based perspective. In this article, we however ignore these value-related issues and focus on the conceptual and practical challenges when defining and measuring resilience for a given ecosystem property or function.

3 Recovery from Individual Perturbations

As outlined above, when referring to resilience as the return rate to equilibrium after disturbances, one can hope to measure resilience by identifying large state perturbations and then estimating the recovery time empirically (van Nes and Scheffer 2007). In principle, this involves two steps: (i) detecting sufficiently large and sufficiently abrupt shifts in order to identify a clear starting time for the recovery process, and (ii) measure the time the system needs to recover. Both aspects need a conceptual choice and quantification of a reference state before the perturbation and after recovery. Closely related to measuring the recovery time, one can also on the same basis measure the recovery rate directly by fitting an exponential function to the part of the time series between the perturbation and the time when the system has fully recovered (Smith et al. 2022). As explained in Arnoldi et al. (2016), this method assumes that the recovery rate does not change during the recovery process, i.e., that the initial and asymptotic recovery are identical. This is only guaranteed in one-dimensional systems, or in case the perturbation happens in the direction of the eigenvector with the largest negative eigenvalue.

3.1 Methods for Detecting Abrupt Shifts

Sudden and large perturbations of a state variable can result from (i) disturbances from which the system can recover (Fig. 2 case 1, blue trajectory), (ii) large disturbances that push the system into an alternative state (Fig. 2 case 2, yellow trajectory) or (iii) “tipping points” originating from the internal nonlinear dynamics of the system, which can be triggered by only small and gradual changes in external conditions (Fig. 2 case 3, red trajectory), and can be practically irreversible. The first kind is useful to probe the linear stability (engineering resilience) of the system, whereas the second and third are associated with ecological resilience and result in a transition to an alternative state. In observed complex systems, it is often not clear into which of these (idealized) classes an abrupt shift falls, in particular shortly after the shift has occurred. However, all three are important events to be monitored before they can be assessed further.

Interestingly, methods to detect irreversible shifts and methods to detect large anomalies (temporary disturbances) are closely related (Flach et al. 2017; Bathiany et al. 2020). For instance, integrating a time series (i.e., computing the cumulative sum) will convert an extreme event into an abrupt persistent shift. Vice versa, differentiating (computing the difference between subsequent data points) turns an abrupt transition back into an extreme event. Consequently, a very fundamental and successful detection approach (depending on the length of the time series) works via computing the cumulative sum of a time series and then detecting extreme values (Basseville and Nikiforov 1993).

It is important to be able to detect extreme or abrupt shifts fully automatically, i.e., without the need for visual inspection to identify candidates of events. First, automatic detection allows more objective and reproducible results, and second, EO datasets are increasing in number, length, and resolution, thus increasing the amount of data that needs to be monitored and processed, which would be hardly possible by eye. A similar development is occurring for climate models of climate and parts of the Earth system.

Several detection strategies with an increasingly large number of algorithms have been developed and applied in the Earth sciences to detect abrupt shifts automatically (see, e.g., Aminikhanghahi and Cook 2017). Historically, an important application of such methods was to remove artificial offsets in observational time series, introduced by changes in the measurement process (e.g., when replacing one instrument by another). To this end, so-called change-point detection methods were introduced which detect the points in time where a sudden change in the statistical properties (usually the mean) of a time series occurs. These change-point detection frameworks involve fitting specific parametric models to time series (e.g., using ramp regression as a piecewise linear model; see Mudelsee 2000), and testing for the likelihood that different models are needed to explain the data in each segment before and after a suspected change-point (Reeves et al. 2007; Bai and Perron 2003; Zeileis et al. 2003). Such a task can also be framed as a model selection problem and information criteria can be applied (Jamali et al. 2014; Reeves et al. 2007).

Originally, several assumptions had to be made when applying these methods, which are too restrictive in the case of time series from climate or ecological systems. Consequently, change-point detection methods have been extended to be able to cope with large autocorrelation (Mudelsee 2000; Rodionov 2004, 2006; Wang 2008; Ducré-Robitaille et al. 2003), trends and non-stationary/non-iid residuals (Matyasovszky 2011; Gallagher et al. 2013), seasonality (Verbesselt et al. 2010a), non-Gaussian distributions (Lazante 1996; Chu 2002), as well as measurement and dating uncertainties (Goswami et al., 2018).

Despite these improvements, almost all methods still are restricted to univariate time series. EO data (like simulated data), in contrast, provides spatio-temporal information, allowing for detection methods that focus on a certain spatial scale of interest, and take spatial connections into account in order to improve statistical significance and to classify abrupt shifts in terms of their spatio-temporal structure (Bathiany et al. 2020). Table 1 summarizes some typical examples of applications in the climate sciences, which share the particular target of detecting real abrupt climate shifts (as opposed to artificial shifts in the measurement process).

Table 1 Methods to detect abrupt shifts automatically in Earth system data

Inspired by ecological disturbances, many techniques for detecting abrupt changes in annual time series of optical medium resolution remote sensing data have been developed. Such techniques assume that a sudden shift in the vegetation community, e.g., by fire, will result in a measurable change in the annual spectral reflectance properties. Detecting changes is mostly done by classical change-point detection algorithms, but those methods are increasingly combined with supervised classification models to filter out falsely detected breaks. Examples of such products are available for Canada (White et al. 2017) and Europe (Senf and Seidl 2021a, b).

3.2 Quantifying Recovery After Perturbations

After having identified a perturbation of the system, for instance with the methods discussed in Sect. 3.1, and under the assumption that the system is ecologically resilient (can recover to its old state), quantitative properties of the recovery can be measured in order to quantify the system’s “engineering resilience”. As explained in Sect. 1, in this framework, more resilient vegetation or ecosystems will return to their previous state (or a predefined state) more quickly. In practice, the disturbance can be thought of as a drought, fire, storm, or any other event that leads to an abrupt change in structure and frees resources for the development of new vegetation. In principle, it is the recovery rate which quantifies what we introduced as “engineering resilience” above, i.e., the linear stability of the system’s state in terms of dynamical systems theory. However, we will discuss in this section that it matters how this property is defined precisely (in particular, if and how the time series are normalized), and how it is measured.

The simplest way to measure recovery time is to define a reference state, and measure how long it takes for a given ecosystem to return to that reference after disturbance (Senf and Seidl 2021a, b; Wu et al. 2022; Bastos et al. 2011). This provides a time scale for recovery, which can be interpreted as a recovery rate. The main practical difficulties here are defining the reference state, and estimating the recovery time in the presence of other perturbations, which can make it difficult to estimate when the recovery process has ended. This stands in contrast to the detection methods introduced in Sect. 3.1, which typically do not require defining any baseline state a priori. Another practical difficulty in EO applications is that time series might not be long enough to robustly define a long-term reference state. Circumventing this caveat is the motivation of space-for-time approaches, discussed in Sects. 4.3 and 5.8.

To quantify the reference state the system recovers to (baseline value), one typically uses either the pre-disturbance state, or a baseline value from an undisturbed control treatment (if available). As mentioned in Sect. 2, it is a subjective choice whether the recovery is formulated in terms of absolute values, or normalized with regard to the baseline or the disturbance impact (Ingrisch and Bahn 2018). In any case, it matters whether the reference state can be well quantified using the data sources available. For instance, the recovery in NDVI or similar spectral indices will reach pre-disturbance values quite quickly after disturbances, even though this will be driven by understory plants and early vegetation (for example, see Bastos et al. 2011). Using reference states with more clearly defined ecological metrics (e.g., common definitions of forest based on tree cover and height) might overcome such problems, but require more advanced data. Here, novel developments in remote sensing might help close this gap (e.g., using LiDAR data to look specifically at structural recovery; see Stritih et al. 2023), but these are still experimental and not widely applicable to assess resilience.

In practice, a pre-defined time period after the disturbance is often selected, and the regrowth capacity is evaluated. A wide variety of measures have been applied that compare the conditions after the disturbance to a reference condition and potentially also the condition at the end of the disturbance (e.g., Frazier et al. 2018; Tao et al. 2022; Lloret et al. 2011; Van Meerbeek et al. 2021). These measures express the regrowth relative to the disturbance impact (e.g., the recovery magnitude relative to the disturbance magnitude) or to the reference state itself. Essentially, measuring regrowth in this way corresponds to measuring recovery, but with different definitions of what the reference state is.

The choice of reference pre-disturbance values varies widely across studies, some considering only the single year before disturbance (Thurm et al. 2016; Sohn et al. 2016) as reference, but most using the mean values of the 2–3 years prior to disturbance (Yi and Jackson (2021). In such short periods, it is likely that the state variable is influenced by the natural interannual variability of many ecological systems (Zhu et al. 2017; Li et al. 2022), risking to lose representativeness. Therefore, in more recent studies, longer-term reference periods have been used to quantify impacts (Flach et al. 2018; van Schaik et al. 2018; Bastos et al. 2021), yet long time series prior to disturbance might not always be available.

Another problem with using recovery time as an indicator of resilience is that it assumes that recovery will be complete after some fixed time. However, systems might take decades to centuries to reach pre-disturbance values, especially for systems composed of long-living organisms such as forests. The observation length of many data series is hence shorter than the recovery interval. This makes it especially difficult to assess changes in recovery over time, which would be of interest to detect changes in resilience caused by, e.g., climate change. Early indicators of recovery might already serve as a proxy for the long-term developmental trajectory, which yet remains to be empirically validated (Seidl and Turner 2022).

Instead of measuring the time span until a “recovery” criterion is met, it might be advantageous to directly quantify the recovery rate by fitting an exponential function to the section of a time series after a disturbance (Fig. 3) (Smith et al. 2022; Bastos et al. 2011). Based on an assumed linear dynamical systems framework, the recovery of a system after a perturbation should be exponential. Hence, the exponent of an exponential function fitted to the time series of restoration directly quantifies the recovery rate of an ecosystem.

Fig. 3

Observed perturbations and recovery in vegetation optical depth (65.625°N, 155.875°W) as measured by passive microwave data (Moesinger et al. 2020). Top: full time series from 1992–2016, raw time series in black, with deseasonalized and detrended time-series residual in blue. bottom: period from 1999–2012; Large disturbances as those in 1999 and 2008 are an opportunity to estimate the recovery rate using an exponential fit

4 Recovery from Continuous Perturbations

The problem with estimating resilience via direct estimation of the recovery rate is that one needs time series that exhibit large perturbations and that it will practically be close to impossible to measure resilience changes over time in this way. In this Section, we will discuss time series that show no large disturbances but permanent variability. Though the resilience indicators obtained from such permanent small perturbations are arguably more indirect than the recovery from single large perturbations, theory predicts a relationship between them and the advantage is that we do not need to rely on the presence and detection of large perturbations as discussed in Sect. 3.

4.1 Statistical Resilience Indicators

In the absence of large externally-induced perturbations, the internal variability of the system can also be used to gain insight in the resilience of a system to external influences. That is because the internal ‘fluctuations’ of a system typically contain information about the system’s response to external influences, as stated by the fluctuation–dissipation theorem (Kubo 1966) in statistical mechanics. It is worth noting that this link has also been observed in climate models (Cionni et al. 2004; Gritsun and Branstator 2007; Leith 1975). Hence, interpreting a system’s variance and autocorrelation (in short, AC1, using a time lag of one resolved time step) as indicators of linear stability has become a widely adopted and straightforward approach across various fields.

Such an empirical approach is particularly relevant when monitoring systems that might be losing resilience, or even approaching a bifurcation point. When moving toward such a point, the dynamics of the system will typically exhibit critical slowing down (CSD). The initial discovery of this phenomenon can be traced back to statistical physics (Gaspard et al. 1995; Solé et al. 1996; van Kampen 1981). Subsequently, it has been investigated in the field of ecology (Wissel 1984; Oborny et al. 2005) and climate research (Kleinen et al. 2003; Held and Kleinen 2004). Critical slowing down refers to the situation where a system takes longer to return to its equilibrium state following a perturbation, indicating a decrease in its recovery rate. In this regard, indicators of CSD aim to detect a loss in “engineering resilience” (decreasing linear stability) as a proxy of closeness to a bifurcation point. It is the crossing of such a point that finally causes a catastrophic transition to an alternative state (loss of Ecological resilience), and therefore typically is seen as an impact-relevant property.

When a system approaches a bifurcation-induced tipping point, its dynamics tend to resemble a red noise process, wherein each state becomes more similar to its past state, leading to an increase in autocorrelation (Wiesenfeld 1985). Simultaneously, the slowness of the system causes a rise in variance as subsequent random fluctuations accumulate due to the slower and slower decay of these fluctuations (Wiesenfeld and McNamara 1986). Consequently, both variance and autocorrelation progressively increase until the system reaches a bifurcation point. This theory has prompted the use of rising variance and autocorrelation as “early warning signals” for impending bifurcations (Scheffer et al. 2009). Additionally, other related indicators have been proposed, such as higher statistical moments (Guttal and Jayaprakash 2008a) and spatial correlations in spatially extended systems (Dakos et al. 2010; Donangelo et al. 2010; Guttal and Jayaprakash 2008b).

There are several practical and conceptual limitations associated with these statistical resilience indicators (Bathiany et al. 2016a). First, practical limits exist because these methods require long time series with sufficiently slow and gradual change in the external forcing (i.e.,, a time-scale separation between external forcing and internal system dynamics). Second, on a conceptual level, these indicators might lack predictive power, because they might detect false positives due to, e.g., multiplicative noise (Bathiany et al. 2013), negative feedbacks in most of the measured regime, or inertial effects (Wagner and Eisenman 2015; Boulton and Lenton 2015; Bathiany et al. 2016b), the confounding effect of colored or non-stationary noise (Dakos et al. 2012), and the large spatial complexity of real-world systems (Bathiany et al. 2013; Held and Kleinen 2004). As an example, if a given climatic driver shows changes in its AC1, then the ecological metric will too (assuming sensitivity remains), which would be reflected in resilience change, while functioning was actually maintained.

There have been some recent advances however to mediate some of these problems. In particular, Boers (2021) introduced a new restoring rate parameter λ as a resilience indicator that is not prone to false alarms based on changes in the properties of the driving noise when using a more general regression method than autocorrelation (see also Boettner and Boers 2022), while Clarke et al. (2023) analyze the ratio of spectra (ROSA) of the system state variable relative to a forcing variable. Moreover, spatial information from EO records can be increasingly beneficial by allowing (i) to generate more data about systems than we would have in case of univariate time series (often by making a “space-for-time” assumption), and (ii) to assess spatial resilience indicators that can complement univariate indicators like autocorrelation, variance or feedback parameter λ (Dakos et al. 2010; Boulton et al. 2022).

Recent work has further linked the statistical measures of resilience to the recovery rate discussed in Sect. 3, over a globally distributed set of ecosystem recoveries. Smith et al. (2022) found that the theoretically expected relationship between AC1/Variance and recovery rate broadly held in empirical vegetation EO data, indicating that AC1 and variance can indeed be related to the recovery rate after a disturbance. This implies that AC1 and variance can be used to estimate how fast a system will recover after a shock, and provide another potential avenue to estimate recovery times and, hence, changes in resilience.

4.2 Pre-processing Time Series to Derive Resilience Indicators

Estimating the resilience of a system using Critical Slowing Down indicators (e.g., the restoring rate λ, lag1-autocorrelation (AC1) or variance) generally requires that the time series describing the system is stationary—that is, without seasonality or trends. Stationarity helps to avoid conflating normal seasonal variations with changes in short-term oscillations, which are signals of changes in system resilience, and avoids biases in the resilience indicators induced by trends in the time series. Therefore, data has to be pre-processed before statistical resilience indicators can be computed from it. It has been shown that the methodology to detect so-called anomalies is typically very important for the detection of extremes (Flach et al. 2017) as well as the estimation of resilience indicators and their trends (Smith and Boers 2023b). We note that these procedures to remove e.g., seasonality and trends are not necessary for estimating recovery times after large and rare perturbations as those discussed in Sect. 3, which focuses on the actual recovery after a perturbation (Verbesselt et al. 2010a).

A wide variety of approaches exists to remove seasonal and long-term trends from data (Smith and Boers 2023b; Verbesselt et al. 2016); the most suitable method is driven by the data set in question. For example, highly seasonal data such as NDVI or other vegetation time series can be decomposed into a trend and seasonal component with harmonic models treating seasonality as a combination of sine waves (e.g., BFAST, Verbesselt et al. 2010a; Verbesselt et al. 2010b, Masiliunas, 2021). There is no one optimal deseasoning method; the efficacy of different procedures depends strongly on the underlying signal that is to be isolated. In addition to the deseasoning method itself, additional parameters such as the length of any seasonal or annual smoother, the choice of a stable reference period with which to construct anomalies, and the number of overlapping sine functions used to represent annual seasonality can have a strong influence on the resultant deseasoned residual.

The type of long-term trend removed is also important to consider—removing a linear trend from non-linear data can leave undesirable traces in the time series that will influence estimates of changing resilience (Smith and Boers 2023b; Verbesselt et al. 2010a; Verbesselt et al. 2016). If, after detrending, the time series still contains a trend, estimating AC1 and variance (or other resilience metrics) based on a rolling window will artificially inflate local autocorrelation, as compared to a locally flat time series. Thus, a trend that fluctuates through time can induce changes in AC1 which might be confused with AC1 changes driven by high-frequency fluctuations or changes in the system dynamics. The impacts of these combined deseasoning and detrending methods can be seen in Fig. 4; not all methods perform equally well at reconstructing a stationary residual, and each leave distinct traces in the inferred changes in AC1 through time. It is important to carefully consider both the means of removing seasonality and trend, and how that procedure can influence resulting estimates of system resilience.

Fig. 4

Impacts of deseasoning/detrending procedures on a simple toy model (from Smith and Boers 2023b, Fig. 1) for a system moving Toward a bifurcation induced transition. Left column: seasonality modeled using a regular sine curve, right column: seasonality varies in amplitude and timing between years. A, B Time series with (black) and without (blue) noise. C, D Perfect deseasoning via removing the a-priori known seasonality curve. E, F Deseasoned time series via a harmonic fit and rolling mean (red), via long-term daily means and a linear detrender (black) and via STL (purple). Gray shading covers periods where rolling mean and STL detrending is uncertain due to incomplete time windows. G, H AC1 computed over a five-year rolling window for each method of removing seasonality and trend. While all models show an increase in AC1 as a transition approaches, not all methods preserve the monotonic increase in AC1 from (A, B). Dashed lines indicate less reliable AC1 values due to incomplete moving time windows, shaded areas cover one STD from the mean. Figure license: https://creativecommons.org/licenses/by/4.0/

There also exist methods which remove both trend and seasonality in tandem, such as Seasonal Trend Decomposition by Loess (STL) (Cleveland et al. 1990). STL has the benefit of removing more complex seasonal signals (i.e.,, not sine waves) and dynamic trends (i.e.,, changing trends through time). However, overfitting, e.g., by fitting high-order splines, can remove portions of the dynamics of the system that are actually part of changing system dynamics—and hence potentially resilience, and not of the seasonality or trend (Fig. 4).

It should be emphasized that no method of deseasoning and detrending is perfect—each will have limitations that need to be considered (Smith et al. 2023; Verbesselt et al. 2016). The non-stationarity of the resilience indicators that one desires to measure can often not be separated from other non-stationary features because they occur on similar time scales. For example, if seasonality strongly increases through time, a single harmonic model fit to the whole time series will leave traces in the residual, as it is focused on removing the average signal through time. This effect could be removed by locally fitting harmonic models (e.g., every three years), with the potential downside of leaving step-changes in the residual year-on-year. In short, the potential limitations of deseasoning methods—and their knock-on effects on resilience estimates—need to be carefully considered. Therefore, careful attention is needed when fitting time series in order to explore the residuals.

4.3 Multivariate Approaches and Spatial Resilience Indicators

In the previous sections, we focused on the analysis of time series of a single state variable, so-called univariate time series. However, in practice, often more data is available, be it time series of multiple state variables or time series of spatial evolution.

Multivariate data is useful because natural systems are complex systems that have many processes acting on, and interacting across, different time scales and affecting different state variables and/or locations differently. For instance, the growth rate varies from one plant type to the other, and soil water dynamics play out on shorter time scales compared to plant dynamics. Multivariate data therefore ‘samples’ the full variety of processes in a system better compared to univariate data, which can lead to better resilience estimates. Mathematically, this means that multivariate data enables the estimation of multiple modes (“eigenvalues and eigenvectors”) compared to univariate methods that are typically limited to one mode only.

To use multivariate data, often the univariate methods from the previous sections are extended to a multivariate setting using mode decomposition methods to derive eigenvalue–eigenvector pairs from a time series. One of the most widely-used techniques for this is a principal component analysis, for which a singular value decomposition (SVD) of the multivariate (i.e.,vector-valued) time series $X$ forms the basis. The downside of the use of such methods is that the time ordering of the original time series is completely disregarded when constructing the eigenvectors. It is, however, possible to explicitly account for this time ordering, for example, using dynamic mode decomposition (DMD) (Kutz et al. 2016). A detailed and practically focused methodology can be found in e.g., Weinans et al. (2019).

Once such a mode decomposition has been obtained, the system’s resilience can be derived from it. A naive implementation would be to look at the least stable eigenvalue–eigenvector pair (i.e., the eigenvalue with the highest real part) as this mode corresponds to the slowest recovering dynamics. However, this mode might not always correspond to the relevant dynamics. For example, when the eigenvalue has imaginary parts it could correspond to oscillatory dynamics, such as seasonality (see also Sect. 4.2). If this is not the case, it could also happen that the system has very slow dynamics that are always present but do not contain information on the system’s resilience. For instance, in semi-arid areas, vegetation patterns migrate slowly over a landscape (DeBlauwe 2012; Bastiaansen et al. 2020) but that dynamics has little to do with the restorative dynamics of the vegetation. Hence, looking at all modes might be necessary to obtain the proper resilience signal.

Proper analysis of these time series using multivariate methods can lead to more accurate resilience estimates and/or reduce the required time span of the data at hand (Weinans et al. 2019; Weinans et al. 2021; see also e.g., Bastiaansen et al. 2021 for exposition of this in the context of global climate models). A number of methods have recently also been further developed and assessed that allow for the automatic detection of abrupt shifts (Bathiany et al. 2020), extreme events (Flach et al. 2017) and changing resilience indicators (Weinans et al. 2021) in multidimensional systems. There is particularly high potential to look at multivariate recovery indicators from remote sensing, with some recent evidence in the literature that multivariate indicators could indeed improve the assessment of post-disturbance recovery (Morresi et al. 2019). When monitoring forests, data on vegetation coverage, tree density, tree height, above-ground biomass, soil moisture and many more might be present (Lechner et al. 2020), and satellite data nowadays typically has high enough resolution to capture spatial patterns, such as those appearing in semi-arid regions (Bastiaansen et al. 2018; DeBlauwe et al. 2008; Veldhuis et al. 2022) or in peatlands (Eppinga et al. 2008). Also, Viana-Soto et al. (2022a) derive fractions of trees, non-treed vegetation and bare ground from remote sensing data, which allows one to study recovery trajectories in a 3D space spanned by those three land covers. A particularly promising approach in this direction is to apply machine learning methods in order to identify optimal indicators (Bury et al., 2021).

When spatial data is available, it is also possible to use the form of the eigenvector instead of the eigenvalue (or in combination with) to estimate the system’s resilience. Here, the idea is that the destabilizing eigenvector of spatially extended systems has some striking characteristics that can be used as a metric (Held and Kleinen 2004; Dakos et al. 2010; Guttal and Jayaprakash 2008b). That is, when a system approaches a bifurcation-induced tipping point, it might do so in a very specific way, thus displaying specific spatial characteristics. Precisely which characteristics to look for varies from system to system. For example, in spatially homogeneous systems the spatial correlation tends to increase when a tipping point is approached (Dakos et al. 2010). In a similar vein, it is also possible to look at the spatial variation or decorrelation length. However, for spatially patterned systems, typically the spatial correlation does not precede critical shifts, and one should look for different signs (Bastiaansen et al 2020; Kéfi et al 2014). For instance, the patchiness of some dryland ecosystems, known for exhibiting spatially patterned states, seems to increase as they approach a critical shift (Kéfi et al 2014).

The strength of these spatial early warning signs is that they do not require lengthy time series to be available. However, they do require a dataset with enough spatial detail. Therefore, in particular for systems where constant monitoring is hard, or for systems that have strong seasonality, infrequent analysis of the spatial structure of a snapshot of the system might be a good alternative.

5 Challenges and Opportunities of EO-Based Resilience Monitoring

5.1 Technological Advances in Observing Terrestrial Vegetation

The challenges of resilience monitoring discussed above in principle apply to all Earth observation records. Before we highlight how these challenges manifest themselves in current datasets and applications, we here provide a very brief summary of recent developments in the retrieval of terrestrial vegetation datasets from Earth observations.

In principle, satellites measure radiative fluxes which by a retrieval algorithm are translated into products that are meant to represent a property on the ground, like vegetation coverage or biomass. The measurement can be passive or active. Passive remote sensing measures the properties of solar radiance reflected by the Earth surface, or radiance emitted by the Earth itself. Active remote sensing systems send artificially generated radiation Toward the Earth and measure its reflection properties, such as intensity and polarization. The latter category includes observation systems such as radar and lidar.

A common and largely available passive EO technology is multispectral imaging, where satellites collect data over multiple wavelength regions, such as the visible, near-infrared, and shortwave infrared part of the electromagnetic spectrum. Multispectral data form the basis of widely used vegetation indices like the Normalized Difference Vegetation Index (NDVI) or the Enhanced Vegetation Index (EVI). They are closely (though not linearly) related to photosynthesis and thus plant vigor and productivity. Yet, NDVI as well as EVI saturate in regions with high biomass (Sellers 1985; Liu et al. 2011), and the short wavelengths render them tainted by atmospheric conditions such as clouds and aerosols (Samanta et al. 2012). In addition, seasonal changes in sun-sensor geometry can lead to serious artifacts (Morton et al. 2014).

Radiometry, and its passive-microwave products such as Vegetation Optical Depth (VOD), overcome these obstacles, as biomass density drives the natural emission of radiance from the vegetation itself and the attenuation of natural emissions from the soil (Liu et al. 2011; Konings et al. 2019; Li et al. 2021; Schmidt et al. 2023). Also, improved plant surface temperature measurements from ECOSTRESS (Fisher et al., 2020) allow insights into physiological mechanisms like stomatal control and hence into the water and carbon cycle and land–atmosphere coupling. Radiometry’s drawback is its low spatial resolution. In contrast, the measurement of backscatter from actively sent microwave signals, as done by synthetic aperture radar (SAR) technology allows for resolutions at meter-scale (Konings et al. 2019). SAR provides insights into terrain structure, vegetation canopy structure, and soil moisture. Besides classical multispectral and SAR-based indices, novel indices such as solar-induced chlorophyll fluorescence (SIF) can even better reveal photosynthetic activity and phenology (Zhang et al. 2022).

Finally, recent spaceborne LiDAR missions, such as GEDI and ICESat-2, provide near-global and global sampling of 3D forest information. Initial exploration of those data focused on combining GEDI samples of forest heights (relative height observations) with Landsat and Sentinel-2 images to derive global, wall-to-wall 30 m and 10 m canopy height maps, respectively (Potapov et al., 2021). These global products will enable a global analysis of canopy height time series in the future, and in turn engineering resilience studies on forest regrowth. Nevertheless, the canopy heights should be used critically as potential systematic effects in satellite LiDAR height observations may lead to an underestimation of forest regrowth (Milenkovic et al. 2022). ICESat-2, on the other hand, already provides time series of forest heights as this mission, in contrast to GEDI, samples along the same tracks and seasonally since 2018. Besides forest height, both GEDI and ICESat-2 provide information on vertical vegetation profile metrics such as leaf area index and canopy gaps (Dubayah et al. 2021), but the exploration of those variables in resilience studies is yet to come (see Stritih et al. 2023 for an example).

On the one hand, these novel developments provide us with a richer and much better resolved database on which resilience monitoring can rely. We provide a timely review of recent developments in the next section, with focus on forest monitoring. On the other hand, a number of conceptual and practical caveats need to be resolved in order to harvest the full potential of these developments. To raise awareness of the gap between current practices in the generation of EO data and the requirements of resilience monitoring, we discuss the most apparent and urgent limitations in Sect.  5.4–5.9.

5.2 Monitoring Resilience in Aquatic and Marine Ecosystems

Although examples in the rest of this review focus on terrestrial ecosystems, the notion of ecosystem resilience in general, the potential to monitor it with remote sensing, and the above conceptual and practical challenges all apply to aquatic and marine ecosystems in a similar way. For example, shallow lakes have been a prototypical example for the potential existence of alternative stable states (Scheffer and Jeppesen 2007), which differ in the turbidity of the water. In principle, earth observation by satellite-borne sensors can be useful for probing the resilience of such aquatic systems (Oezkundakci and Lehmann, 2019). A large body of research also focuses on the resilience of coastal wetlands, e.g., by observing the recovery after droughts, hurricanes and other storms, heavy precipitation and fluvial flooding. Such studies make use of a large number of remote sensing sources such as MODIS and Landsat (Tahsin et al. 2018), vegetation-based indices (Balogun et al. 2020) and sometimes use the classical statistical resilience indicators as outlined in Sect. 4.1 (Zhao and Li 2022). Moreover, similar detection methods to automatically detect abrupt shifts as outlined in Sect. 3 are being developed and applied in the marine sciences (Yang et al. 2022). Other ecosystems under investigation are giant kelp forests (Cavanaugh et al. 2019), coral reefs (Knudby et al. 2014) and seagrass ecosystems (Clemente et al., 2023), with particular concern of tipping in these systems (van de Leemput et al. 2016; van der Heide et al., 2011). While returning to terrestrial applications in the remaining sections, we point out that some of the same challenges as discussed throughout this review also apply to aquatic and marine ecosystems as well.

5.3 Developments in Earth Observations with Regard to Terrestrial Ecosystem Resilience Monitoring

In recent years, Earth observation has become a major tool in assessing resilience (Nikinmaa et al. 2020; De Keersmaecker 2022; Senf 2022), often involving the technological advances outlined in the previous section. Special focus has been put on forest ecosystems, because forests globally are changing in response to climate change and increasing land use pressures (McDowell et al. 2020). For instance, recovery from disturbance such as windthrow, bark beetle outbreaks or fire have been estimated from optical remote sensing time series in order to assess the forest’s capacity to return to similar spectral conditions than before disturbance (White et al. 2019; White et al. 2022; Chirici et al. 2020). Most of those studies made use of free image archives such as Landsat (Wulder et al. 2019), which provides long enough time series to estimate recovery signals in forests. The use of Landsat also allows for the quantification of historic disturbance events (Senf and Seidl 2021a, b; Hansen et al. 2013), which in turn can be used to estimate recovery by space-for-time substitution. More specifically, combining the historic information on disturbances with recent data on forest structure, such as derived from airborne or spaceborne Light Detection and Ranging (LiDAR) data, allows to create chronosequences of structural attributes (e.g., canopy cover, height) over variable years post-disturbance (Viana-Soto et al. 2022b; Stritih et al. 2023; Senf et al. 2019). From those chronosequences, it can be estimated how quickly forests return to certain structural states or recover certain functions, which substantially increases the applicability of the result in ecosystem management (but see also discussion in Sect. 5.4). Similar results can be achieved also without the pitfalls of space-for-time substitution (Johnson and Miyanishi 2008) by recent developments in Earth observation, including machine learning regression (Senf and Seidl 2021a, b), synthetic spectral unmixing analysis (Senf et al., 2020) or artificial intelligence (Liu et al., 2023). Those novel approaches convert the spectral signals into interpretable units, such as canopy cover or biomass, allowing for example for the assessment of community change after fire disturbances (White et al. 2023; Viana-Soto et al. 2022a). Novel sensor systems, including spaceborne hyperspectral data or Radar, will likely add important information to ecosystem resilience, but their full potential hasn’t been explored yet.

5.4 Signal Versus Functional Recovery

Many remote sensing based recovery studies use signal recovery as a measure of resilience. Signal recovery is derived from a time series of a single remote sensing observable, being it a spectral index, e.g., NDVI or RVI (the radar vegetation index, Saatchi 2019a, b), or a biogeophysical variable derived from the radiance measurements, e.g., the fraction of absorbed photosynthetic radiation (fAPAR) derived from solar-reflective measurements, or vegetation optical depth (VOD) from microwave observations. Yet, by focusing on a single ecosystem property, signal recovery does not always reveal to what extent the same level of processes or structure are retained in the system. For instance, while after a disturbance the NDVI signal of a tropical forest may recover to its pre-disturbed state, the secondary forests are known to be structurally and functionally different (Caughlin et al., 2021). In addition, signal recovery does not provide information on the process behind the recovery. For instance, monocultures can be more sensitive to diseases but recover relatively fast when they are replanted. Yet, the fast signal recovery capacity of these monocultures does not necessarily reflect their resilience to disturbances. Consequently, it is important to find a proxy variable measurable by remote sensing which represents the in-situ property under investigation. If this is not directly possible, multivariate approaches based on several, independent remote sensing observables, and additional knowledge on the relationships of time scales may help constrain the recovery process quantitatively.

5.5 Spatial and Temporal Resolution

The scale of EO measurements rarely matches up with the system under study, neither in space nor in time. For example, Landsat has a theoretical revisit time of ~ 16 days (or much longer in areas with frequent cloud cover), which is too low to capture some of the dynamics of many ecosystems. To track phenological changes, or the response to extreme events like droughts or storms, daily resolution would be required. Conversely, MODIS with a daily repeat time lacks the spatial resolution to resolve small-scale features on the spatial scale of individual trees that can be picked out with Landsat data (30 m resolution). The trade-off between time and space becomes even more of an issue when very-high-resolution optical and Synthetic-aperture radar (SAR) are used, which, depending on the satellite constellation and priority of the ground target in the acquisition plan, may provide only a few acquisitions per year. It is not always clear what the optimal spatial or temporal resolution to study a given system is, and to what extent EO data captures the dynamics well enough to estimate changes in resilience through time. Multi-scale approaches, i.e., tackling a problem using various spatial and temporal resolutions, might make resilience assessments from EO data more robust. Also, it is strongly recommended to first define the scales necessary to understand a system and then choose the appropriate EO dataset, and not vice versa (Senf 2022).

5.6 Data Gaps

Gaps are common in EO data, e.g., due to cloud cover, smoke, snow and frozen ground, fog, and in the case of radar and very-high-resolution optical imagery, because of prioritization in the acquisition plan. For optical data, gaps are particularly a problem in the tropics, where cloud cover is frequent and persistent. For microwave data, large seasonal gaps may occur when the soil is frozen, covered with snow, or inundated. Gaps can cause problems both in theoretical resilience metrics (e.g., AC1, variance) and engineering metrics (recovery time). For example, gaps in a time series can bias AC1 and variance if the number of measurements over a given time period fluctuates, or if there are sporadic large jumps through time (Ben-Yami et al. 2023; Smith and Boers 2023b). Engineering resilience metrics can also be biased by missing data by either under- or over-estimating the length of a recovery due to missing data at the start or end of the perturbation or recovery period. More generally, deseasoning and detrending methods do not always take missing data into account, and can yield inconsistent results when data is sparse (Smith and Boers 2023b; Verbesselt et al. 2016). Many statistical methods can handle time series with data gaps without gap filling (e.g., by computing variance or autocorrelation from only the existing values). However, the gap properties (gap size, frequency etc.) should certainly be taken into account carefully when preprocessing time series to derive robust resilience indicators. For example, masking an artificial stationary time series with the same gap structure as in the observed data can give an estimate of gap-induced biases (Ben-Yami et al. 2023).

5.7 Time Series Length

Most resilience methods require long time series, i.e., at least long enough to cover the recovery of a single disturbance event or even a multiple of the decorrelation time scale when using statistical resilience metrics like the recovery rate λ, AC1 or variance. Especially for systems with slow internal dynamics, such as forests with long regeneration times, this might be problematic given the relatively short length of the satellite record in general, let alone individual instruments. Besides, the time series length should also allow for a joint analysis of the recovery and disturbance regime. For example, as outlined in Sect. 2.2, a system that takes 50 years to recover and is subject to a lower disturbance frequency will be able to restore most important properties, while a system with a higher disturbance frequency than the recovery period will likely have a higher chance of shifting into a state change (Senf and Seidl 2021a, b). Finally, a sufficiently long time series is required to separate cyclic changes from long-term change (see also Sect. 4.2).

However, because EO data acquisitions have become widespread and frequent only recently, many available data sets are relatively short. Very few of the early remote sensing instruments are still being used, and most missions have ended or been supplanted by newer satellites. Thus, for long-term time series analysis such as the quantification of resilience changes, data from multiple satellites need to be considered. This requires a thorough homogenization of random and systematic differences between instruments, either at the level of radiometric or spectral properties (radiances, reflectances, backscatter, or brightness temperatures; e.g., Xiong et al. 2020) or starting from biogeophysical variables derived from single sensors (e.g., Moesinger et al. 2020; Preimesberger et al. 2021).

5.8 Space-for-Time Substitution

Even though most EO products still have insufficient length to robustly assess ecosystem resilience, it is now getting easier to obtain high-resolution spatial data. Consequently, it is an appealing idea to use a space-for-time substitution. The idea to infer temporal information from spatial differences is not new; for instance, different succession stages in forest development can be identified by looking at species composition at different sites (Picket, 1989). Besides such qualitative comparisons, it is also suggested to use a similar space-for-time approach for quantitative studies (Blois et al 2013), which might also be suitable for resilience monitoring. Here, one could for instance study spatial data along an environmental gradient to get insight into the temporal evolution in accordance with the changing conditions (e.g., to assess Amazon rainforest regrowth; Milenkovic 2022). The spatial resilience indicators as outlined in Sect. 4.3, which are justified from the theory of stochastic dynamical systems, are one promising approach in this context—at least in case their validity to describe the resilience of complex ecosystems can be confirmed. Another example of a space-for-time approach is the diagnosis of multimodal distribution of tropical tree cover for a similar range of mean annual precipitation (Staver et al. 2011; Hirota et al. 2011). If the tree cover data can be understood as a representation of the same stochastic system (space-for-time replacement applies), this finding implies the existence of multiple stable states under the same annual rainfall conditions, and hence potential tipping points in tropical forests. Caveats related to the retrieval and classification of tree cover from satellites in this context have been largely resolved (Hanan et al., 2014; Staver and Hansen 2015).

However, there are also caveats and doubts regarding the validity of space-for-time replacements (Damgaard 2019). Most prominently, the fundamental assumption behind space-for-time methods is that temporal and spatial evolution are comparable and have the same drivers, and that other environmental factors, that differ spatially, can be neglected. This assumption is likely to not hold under present-day conditions where CO₂ fertilization, climate change and nutrient deposition play together, so that systems are not close to equilibrium. Also, even in equilibrium, there are usually more possible explanations for multimodality of the data distribution than the existence of multiple stable states. For example, a bimodality in tree cover states can also be imposed by bimodal environmental effects like human activity (Wuyts et al. 2017), other climate variables than annual mean precipitation (Good et al. 2016) or by nonlinear and transient growth dynamics (Yin et al. 2014). Also, the noise of a system might often have a higher spatial correlation compared to its temporal correlation. As far as we are aware, there are no comprehensive theoretical studies of this assumption in the context of resilience monitoring. Future research in this direction might reveal the merit of such approaches.

5.9 Impact of Data Artifacts on Resilience Metrics

The need to reduce gaps in satellite datasets (see Sect. 5.6) and to generate longer time series than the lifetime of a satellite, has motivated the construction of datasets from multiple instruments. In resilience research however, multi-instrument records pose serious challenges to the methodology used to detect change, e.g., how much of the change can be attributed to environmental shifts, and how much is simply due to the use of more precise instruments, different satellite retrieval schemes, changes in satellite overpass times, or other measurement changes? harmonization efforts (e.g., Moesinger et al. 2020; Preimesberger et al. 2021) usually focus on mean and variance to enable the detection of trends (e.g., changes in leaf area index through time). Unfortunately, harmonization of higher-order statistics such as autocorrelation has received less attention. This implies that satellite datasets may be biased even though they have been evaluated in terms of mean values, and extra care must be taken to ensure that changes in resilience indicators can really be attributed to environmental change (Smith et al. 2023). Further collaboration is needed between environmental scientists and the data curation and creation community to document the potential impacts of changing satellite mixes on resilience indicators, and to better document what statistical measures can be safely evaluated up to what precision.

Artificial data features that can remain after insufficient preprocessing can also be a cause of an erroneous resilience assessment. In principle, these are general issues with remote sensing products, but can be particularly relevant for resilience estimators even if the products might show unbiased mean values. For example, when aggregating overlapping images from neighboring orbits, the spectral properties of the same object may be different just because the images have different looking angles to the object (Dorigo 2012). This may cause artifacts in the time series, e.g., additional noise, and ultimately lead to a false detection of change. Thus, it is important to report what radiometric and geometric corrections were applied to the satellite data to ensure that the data is analysis ready. Sometimes, particularly in high resolution synthetic aperture radar data with strongly oblique view angles, a ground object itself causes artifacts. For example, when a patch of vegetation is completely removed, e.g., after a clear-cut event, pixels along the patch edge in the near-field of the radar satellite remain in the “shadow”, i.e., are occluded by the neighboring forests. Such radar image artifacts cause an abrupt drop of the signal, which facilitates the deforestation detection (Bouvet et al. 2018), but does not recover, because the pixel remains in the radar shadow of the neighboring forest. Including such pixels in the analysis may lead to a false interpretation of forest recovery.

6 Conclusions, Recommendations, and Outlook

We have outlined conceptual and practical challenges of monitoring ecosystem resilience and resilience changes from EO data. We first discussed the long-standing diversity and subtlety in defining resilience conceptually, beyond the general classes of ecological versus “engineering” resilience. Given the large variety of resilience quantifiers, particular care should be taken in future research articles when naming and describing the particular approach taken. In this regard, “resilience” should be considered a topic or a qualitative group of system properties, and distinguished from particular quantifiers, which should be named more precisely and with more care.

We also discussed methods to detect and quantify abrupt shifts (in ecological systems: disturbances), the recovery time and statistical indicators of linear stability (recovery rate, AC1, variance, spatial correlations and others). Besides the need to acknowledge the general limitations of such methods (like the validity of the stochastic system or statistical robustness), we emphasized several additional caveats arising from data uncertainties. Many of these uncertainties arise from combining different sensors (Smith et al. 2023), from gaps in the time series (Ben-Yami et al. 2023), and from detrending and deseasoning the data (Verbesselt et al. 2016; Smith and Boers 2023b).

It has consequently been shown that for some vegetation types resilience measures are not necessarily consistent with each other even given the same data (Smith and Boers 2023b), which points to potential data uncertainties or simply unsuitable EO data for quantifying the relevant aspects of vegetation dynamics, preventing the use of these measures.

Despite the growing availability of EO data, there remain several key limitations to its use in monitoring terrestrial ecosystem resilience and its changes, particularly in near-real-time applications, maintaining a gap between the development of novel remote sensing techniques and their applications for ecosystem resilience monitoring (also see Senf 2022). Many of these limitations are areas of active research, and can be ameliorated through improved models, data processing schemes, and longer time series.

Based on our assessment above, we make the following key recommendations:

• 1. There is need for more clarity in resilience studies on (i) the definitions of resilience used, (ii) what the ecological metric actually refers to (vegetation cover, biomass, composition, function, or other), and (iii) the careful choice and clear definition of reference states considered.

• 2. Innovations are needed in ways for dealing with inherent seasonal and phenological variations in EO time series (e.g., indices such as NDVI). Some systems have high inter-annual variability, and awareness and best practices are required to deal with it.

• 3. Large opportunities arise from the availability of higher spatial resolution imagery such as the Sentinels (1 and 2), and future Very High Spatial resolution satellite missions (e.g., Home—Copernicus Contributing Missions Online). Key here is to see which methods are able to measure resilience by analyzing spatial variability (heterogeneity), establish whether space-for-time assumptions can be justified, and further explore the possibility of estimating resilience via the potential of artificial intelligence, making use of the multitude of geo-spatial data, on different aspects of the ecosystem, that are currently available (Smith et al. 2023).

• 4. Time series are becoming longer, and are also available in near real time. They should be processed carefully regarding seasonality and the treatment of gaps and other uncertainties, and the limitations should be documented in a way that allows for quantitative uncertainty estimates of statistical indicators.

• 5. The suitability of remote-sensing datasets for resilience analyzes should be better verified, especially for high-biomass regions like the tropical forests. As biomass might be the most indicative parameter of rainforests, radar and radiometry in the very low frequency range (P- or L-band) should be desirable due to the deeper penetration from below the canopy. Here, upcoming satellite missions such as the European Copernicus Biomass and American NISAR mission are promising. Furthermore, future availability of in-situ measurements, for example based on terrestrial or airborne LiDAR campaigns, are recommended to evaluate the different data sets, especially in these regions.

• 6. Novel measurements such as land surface temperature (e.g., the future Copernicus Land Surface Temperature Monitoring satellite mission, NASA’s CIRiS CubeSat, and commercial space based alternatives), vegetation structure (e.g., vegetation height via GEDI, ICESAT-2, and future BIOMASS mission) and fluorescence (e.g., via the future ESA FLEX satellite mission) provide direct biophysical variables to monitor terrestrial vegetation dynamics of the ecosystem. Before directly utilizing them in describing ecosystem change and recovery, one should first understand the uncertainties and limitations associated with those novel products as they are often optimized for other applications than resilience monitoring, or the pre-processing is simplified to foster big-data handling and processing.

• 7. Knowledge transfer and closer collaborations between marine, aquatic and terrestrial ecologists and remote sensing experts may create synergies in monitoring ecosystem resilience in different realms.

Overall, we conclude that it is beneficial that there are more and more applications of different resilience indicators (“early warning signals”) to Earth observations (EO). These observations can provide advantages, e.g., higher spatial coverage and hence spatial indicators, potentially by leveraging space-for-time relationships. However, it is vital to be aware of the caveats associated with data generation.

We also note that the research community designing retrieval algorithms and generating EO datasets on the one hand, and the community interested in diagnosing and interpreting ecosystem resilience (often ecologists, climate scientists and applied mathematicians) on the other, are still not in sufficiently close scientific contact. This brings the risk of generating datasets without a well-defined scope of use, and the risk of (mis)using satellite datasets for resilience monitoring despite technical and statistical biases.

We hence see the need for closer collaborations that help close the gap between what EO data can offer and the requirements of EO users. Specific questions of this dialog should at least comprise: What variables are needed for what purpose? On what resolution? Which statistical properties should be optimized, and how are biases communicated? How should gaps and multi-sensor products be compiled by EO scientists? We specifically recommend better transparency and documentation in this regard. Only then can the value of long and frequent observations from space live up to their great potential, and help us get a clear view on the resilience of the Earth system, and help society navigate into a safe future.