A review of measuring ecosystem resilience to disturbance

Resilience is the central concept for understanding how an ecosystem responds to a strong perturbation, and is related to other concepts used to analyze system properties in the face of change such as resistance, recovery, sustainability, vulnerability, stability, adaptive capacity, regime shift, and tipping point. It is extremely challenging to formulate resilience thinking into practice. The current state-of-art approaches of assessing ecosystem resilience may be useful for policy makers and ecosystem resource managers to minimize climatological or natural disaster related impacts. Here, we review the methods of assessing resilience and classify and limit them to three cases: (a) forest resilience based mainly on remote sensing and tree-ring data; (b) soil microbial community resilience based on laboratory and field studies; and (c) hydrological resilience of terrestrial biomes based on the Budyko framework and climate data.


Introduction
The magnitude and frequency of climate-related extreme events, such as heat waves, extreme drought, and flooding (Coumou and Rahmstorf 2012, McPhillips et al 2018, Kornhuber et al 2019, is increasing as CO 2 continues to rise and the climate continues to warm (Hansen and Sato 2016). Global warming is a fundamental cause of the increase in extreme climate events as predicted by the first and second laws of thermodynamics, as a consequence of warmer air that can hold more water vapor (Yi et al 2015). The sustainability of socioecological systems in the face of increasing extreme events has become a global concern (Steffen et al 2018). The resilience concept has been most widely used to describe manifold properties (or abilities, see box 1) of socioecological systems in the face of sudden changes, such as: (a) resistance-the ability to withstand disturbance; (b) recovery-the ability to return to the original state; (c) stability-the ability to retain the same function and structure; (d) vulnerability-inability to withstand disturbance, and (e) adaptive capacitythe ability to deal with change (Walker et al 2004). As a stable socioecological system becomes unstable, the system will shift into an alternative stable state, in which its resilience properties are functionally different from the old one (Ives and Carpenter 2007). Because of the synthetic nature of the concept of resilience, the number of papers relating to resilience has been increasing exponentially in a great variety of interdisciplinary, resilience-related journals (figure 1). The meaning domain of the resilience concept as presented in this literature has become larger but vague. Therefore, what resilience exactly means has become a hot topic for resilience scientists. Various journals have opened special issues to debate and review resilience (e.g. Carpenter et al 2001, Folke 2006, Janssen et al 2006, Millar and Stephenson 2015, Nimmo et al 2015, Reyer et al 2015, Trumbore et al 2015, Meerow and Newell 2016. Despite an evolving literature, resilience remains a useful concept for studying the ability of socioecological systems to cope with extreme events. At the same time, it seems impossible to quantify resilience by a universal index, formula, or analytic approach. Quantifying resilience is notoriously difficult, like measuring 'hard work' by a universal index. No doubt everyone knows what 'hard work' means in a particular context, but it is difficult to compare two or more different jobs by using a formula because 'work' does not hold the meaning of 'force × displacement' as defined in physics. Nevertheless, much effort has been made to develop planning and management approaches Resilience (a) Measure of the persistence of systems and of their ability to absorb change and disturbance and still maintain the same relationships between populations or state variables (Holling 1973). (b) The capacity of a system to manage disturbance by resisting change (robustness), recovering from change (stability), or adapting and benefiting (adaptability judged as beneficial) as a result of change (Helfgott 2015). (c) The ability of a system to recover to the original state upon a disturbance (Scheffer 2009). (d) Measure of the amount of change needed to change an ecosystem from one set of processes and structures to a different set of processes and structures (Angeler and Allen 2016).
Resistance (a) Staying essentially unchanged despite the presence of disturbance (Grimm and Calabrese 2011). (b) The ability to persist during the disturbance (Nimmo et al 2015).

Recovery
The process by means of which an ecosystem bounces back to its pre-disturbance status (Fraccascia et al 2018). Adaptability (a) The ability of a system to adjust, modify or change its characteristics or actions to moderate potential damage, take advantage of opportunities, or cope with the consequences of shock or stress (Brooks and Adger 2004). (b) The ability of a system to adjust to climate change, to moderate potential damages, to take advantage of opportunities, or to cope with the consequences (IPCC 2011).

Vulnerability
The state of susceptibility to harm from exposure to stresses associated with environmental and social change and from the absence of capacity to adapt (Adger 2006). Sustainability The ability to be maintained at a certain rate or level (Gasser et al 2019). Stability The ability of a system to return to an equilibrium state after a temporary disturbance; the more rapidly it returns and the less it fluctuates, the more stable it would be (Holling 1973). Equilibrium The entropy production of an adiabatically insulated system must be zero for equilibrium (reversible transformation) and positive for non-equilibrium (irreversible transformation) (De Groot and Mazur 2013).

Steady state
The state variables of a system are independent of time (Nicolis and Prigogine 1977).

Regime shift
A sudden jump from one dynamic regime to another one (Scheffer 2009). Tipping point Also called critical thresholds, at which the system shifts abruptly from one state to another . Perturbation An effect; the response of an ecological system to disturbance (Rykiel 1985). Disturbance A cause; a physical force, agent, or process, either abiotic or biotic, causing a perturbation in an ecological system (Rykiel 1985). and tools to measure resilience and thereby manage systems to minimize the impacts of disturbance. Generalizing rather than actually quantifying resilience of the overall system, these approaches measure the components of resilience within the system such as resistance, recovery and so on.
In this paper, we attempt to clarify what widely used resilience indicators exactly mean in stability theory, from which Holling (1973) adapted them. We limit the review of resilience measurement methods to three cases, all having to do with terrestrial ecosystem measurement: (a) forest resilience measured by remote sensing and tree-growth observational data; (b) soil microbial community resilience measured by laboratory and fieldwork; and (c) hydrological resilience of catchments as assessed by the Budyko framework and climate data. Here, measuring resilience refers to measuring resilience components, such as resistance and recovery or combinations of them, elasticity, etc. There is no single property of resilience, such as recovery, or a universal index or other variable that encompasses the totality of the concept of resilience. In this review, we only focus on resilience indices relating to resistance and recovery processes of a perturbed departure from a stable steady-state, and do not include those describing an unstable steadystate and the regime shift at the tipping point. The underlying reason for why we selected these three cases is that they are substantially different in measuring methods across the three large research communities of forest, soil, and water.
The selection of a variable to characterize a property of resilience is heavily dependent on data availability and the ability to at least partially quantify the level of ecosystem health. For example, ecosystem productivity is a good indicator of ecosystem health but it can be assessed using different methods or indices such as leaf area index (LAI), normalized difference vegetation index (NDVI), and normalized tree-ring width index (RWI). Even using the same index, different investigators might use different formulations to address their issues on specific ecosystems with specific spatial-temporal scales. Thus, although we will narrow our review to three cases, the diversity of approaches to measuring resilience is still high. In this review, we provide comprehensive lists of different resilience formulations for each category, summarize their calculation steps, and compare their advantages and disadvantages as listed in the tables.

Key definitions and concepts
To better understand what part of resilience we are able to measure by an index or analytical formulation, we must first clarify the derivation of the idea of resilience and related concepts. The term resilience has a history. It first depended on a simple formulation of 'bounce back' (Alexander 2013) before Holling (1973) introduced the concept of resilience to ecology. Holling (1973) defines resilience as a 'measure of the persistence of systems and of their ability to absorb change and disturbance and still maintain the same relationships between populations or state variables.' His definition has been extended and used by various scientific disciplines. Consequently, the concept of resilience has become rich, and the conceptual deviation (or vagueness) is increasing. Many resilience review papers have been dedicated to clarifying the confusion among resilience definitions and meanings (e.g. Carpenter et al 2001, Brand and Jax 2007, Fraccascia et al 2018. Holling's descriptive resilience concept was mainly based on nonlinear system theory (May 1972a(May ,1972b. Many terminologies of nonlinear dynamics such as stability, equilibrium, domain of attraction, and phase space appear throughout Holling's benchmark paper (Holling 1973) (figure 2).
The key contribution of Holling to resilience concepts is the combination of resilience with nonlinear stability theory.
To better understand Holling's descriptive resilience concept with the language of nonlinear dynamics, let x denote state variable of a nonlinear ecosystem. Its change rate dx/dt is governed by where f (x, λ) is a nonlinear function of state variable x, and λ is an environmental control parameter. At steady state x s , the equation (1) becomes Multiple steady states can be obtained by solving the algebraic equation f (x s , λ) = 0. The number of steady states depends on the nonlinearity caused by internal feedback processes. For illustration, let f be simplified as When λ < |λ * | = 1.09, there are three steady states (two are stable indicated by circles and one is unstable indicated by filled circle in figure 3), while if λ > |λ * | = 1.09, a single stable steady state exists. At the threshold value λ * = 1.09, the number of steady states is switched, at what is called a bifurcation point or tipping point. Steady states are fixed points in the phase space (figure 3). The phase space for this nonlinear system is one-dimensional because the system has a single state variable. These steady states are also called singular points (time-independent) and all other points in the phase space are called regular points. The stable singular points are usually called attractors and unstable ones are called repellers. Holling (1973) illustrates phase space in his figure 1 and singular points in his figure 2. The stability of these steady states can be expressed with a potential function. The two stable steady states are the two minima of the double-well potential, and the unstable steady state is the maximum between the two minima (figure 3). The depth of the potential well between a minimum and the maximum represents the resilience of the steady state to perturbation. The deeper a potential well, the greater the resilience of the steady state, and hence Potential Phase space Bifurcation Figure 3. Illustration of potential, phase space, and bifurcation based on the equation (3). The empty circles indicate stable steady-states, filled circles are unstable steady-states, and half-filled circles are the steady-states at the tipping point. The parameter value dictates steady-state number, stability, and switch condition. At λ = 0, the double-well potential is symmetric, which means that two stable steady-states have same resilience (ability) to absorb disturbance. As λ = ±0.5, the biased double-well potential illustrates that the steady state with a deeper well has stronger resilience than the shallower one. For parameter values beyond threshold λ * = 1.09, only a single minimum exists with no potential barrier. The system's evolution can be shown with a phase line that includes all possible states. The convergent arrows indicate the destination (attractor) of system evolving state over time (a decay process of a perturbed state, i.e. recovery process), while divergent arrows indicate the growing process of a perturbation from a repeller (unstable steady-state). As the control parameter increases or decreases to the threshold values, a repeller and an attractor meet and disappear together and the system switches onto a new attractor. In the bifurcation diagram, as the parameter increases and passes the positive threshold (λ * = 1.09), the system becomes unstable and then jumps onto a new stable steady-state (right potential well). Then, the control parameter λ is decreased slowly until it reaches the negative threshold (λ * = −1.09), and the system runs down to the left potential well. More discussions can be found in the main content.
resistance and recovery ability are stronger. If a perturbation is greater than the depth of the potential well, it is likely to bring the system over the barrier of the potential to another stable steady state. The number of steady states is determined by internal feedback mechanisms and the depth shift from deep to shallow (or vice versa) by external control factor λ (figure 3). At the tipping points |λ * | = 1.09, a steady state loses its resilience or its stability and any small perturbation will tip the system to a single stable steady state. The stability of steady states can be examined by a linear stability analysis approach (Nicolis and Prigogine 1977). The perturbed steady state can be written as where x ′ is a perturbation caused by an extreme event.
Assuming |x ′ /x s | ≪ 1, the equation (1) can be linearized by a first-order Taylor expansion, (5) Using the steady state condition dx s /dt = f(x s ) = 0, we obtain the perturbation equation and the perturbation solution is For a stable steady state, the perturbation decays and the system returns to the original state. These stable steady states are called attractors (open circles in figure 3). The absolute value of α(< 0) is called the recovery rate. If |α| is larger, the recovery rate is faster, so the ability of system to absorb perturbations and hence its resilience is stronger (and vice versa). When α > 0, the perturbation increases exponentially and the system runs away from the original state. These unstable steady states are illustrated by filled circles in figure 3. As α becomes less negative, system resilience becomes weaker and recovery speed becomes slower. At the tipping point α = α * (λ * ) = 0, the system loses its resilience and recover time is infinity (i.e. no recovery). From our illustration (3), α = −3x 2 s + 2. Here x s is the set of multiple steady states that are the function of environmental control parameter λ, which can be obtained by solving following algebraic equation Although the above discussion is limited to a small perturbation near a steady state, it is valid in qualitatively understanding of instability behavior of steady states (Nicolis and Prigogine 1977).
Here, we want to clarify concept confusion between equilibrium and steady states used in ecosystem sciences. It is misleading to call the steady state of an ecosystem determined by equation (2) as equilibrium as defined in other fields of scientific study, such as thermodynamics and dynamical systems theory. A biological system or an ecosystem is an open system that is far from equilibrium (Nicolis and Prigogine 1977).
A useful state variable for an ecosystem is entropy. The entropy change dS of an ecosystem during a time interval dt can be written as the sum of two parts, where d e S is the entropy flux due to exchanges of energy and mass with its environment and d i S is the entropy production due to the irreversible processes inside the system. The second law of thermodynamics is expressed as d i S = 0 for equilibrium, at which processes are reversible and hence there are no gradients or fluxes within the system, while d i S > 0 for non-equilibrium, at which processes are irreversible and fluxes are driven by corresponding gradients. The steady state is timeindependent, i.e. dS = 0, thus A life system at high order-level under nonequilibrium conditions is maintained by a sufficient amount of negative entropy flow through energy and mass exchanges. Prigogine (1977) demonstrated that it is impossible for a life system to maintain near equilibrium because entropy production of such nonequilibria is minimal. Therefore, Prigogine (1977) clarified that non-equilibrium may be a source of order in a life system that is far from equilibrium. Equilibrium is one of many steady states for a system, at which there will be no irreversible processes (no mass and heat fluxes or gradients) and the entropy is maximum. The steady state of a life system is far from equilibrium.

Selection of reviewed papers
We conducted our review and analysis in order to compile a list of resilience indices and the published papers in which they were used. We used ScienceDirect to search for papers that empirically quantified the resilience of an ecosystem. The search concluded in August 2018. The search terms were 'resilience index,' 'soil resilience,' 'forest resilience,' and 'ecological resilience'; we then manually examined the first 250 papers of each search by relevance (figure 4). We did not limit our searches by publication year. In some circumstances, we used review papers to locate new index sources. The papers that contained the keywords but did not actively quantify resilience, or were about nonecosystem resilience, were excluded from this review based on abstract and title screening. The papers that examined qualitative definitions or quantified  resilience of a given system through purely theoretical models (i.e. differential equations), screened by examining the content in each 'Material and Methods' section, were also excluded. The papers in the final list were published from 1976 to 2018. The studies were conducted with two principle goals: (a) Identification and characterization of effect modifiers for resilience (b) Quantification of the impact and recovery process on a given area or sample Many studies also analyzed their results in the context of climate change and the changes that the ecosystems in each study would endure in a warmer future.
For each paper, we noted the formulation and the input units (NDVI, BAI, etc) used to set the control, perturbed, and recovered level of the response variable. In addition, we considered the advantages and disadvantages of each formula. After the compiling stage of this review, the final papers were divided into three subgroups based on the system of which resilience was taken: forest resilience, soil microorganism resilience, and hydrological resilience.

State-variable selection
In practice, what measurements can we use to better characterize the state of forest growth? These widely used variables are summarized in the left-box of figure 5. They can be classified into two main groups based on: (a) tree-ring width measurements, and (b) satellite-derived spectral indices of vegetation. The net primary production (NPP) should be the best description of forest growth. A forest NPP can be modeled but not directly measured. The treering width provides a key part for assessing BAI that can be used to estimate forest stocks (Seidl et al 2012). The BAI has been widely used for forest resilience assessments since the most popular approach was developed by Lloret et al (2011). Some forest resilience studies have directly used the measured tree-ring width (e.g. Rust (2015) used the relationship of RWI with the standardized precipitation-evapotranspiration index (SPEI) to be able to identify the threshold SPEI for conifer mortality (Kolb 2015). Yi et al (2018) used the bifurcation characteristics of RWI along with climate to distinguish what climate patterns are healthy or unhealthy for forest growth. No doubt, tree-ring width is a good indicator for forest growth with two distinguished advantages: (a) high resolution of annual variability, and (b) long-term data availability. The limitations of using tree-ring width as a state variable are obviously: (a) the signal of the annual ring of tropical rain forest trees is too weak to use (Lieberman et al 1985, Whitmore 1998 In comparison with the tree-ring-based indices, the satellite-derived vegetation indices as forest-resilience state variables are more powerful in monitoring forest recovery and initial impact of disturbance due to their continuous spatial coverage and higher temporal resolution.

Steady state
How can we measure steady states? Theoretically, steady states are defined by equation (2) and independent of time, which is an idealization that does not exist in reality (Rykiel 1985). In practice, these steady states are approximately equal to average values of a selected variable, i.e. x s ≈x. This is why an average of multiple pre-or post-disturbance event years of state variable data is widely used in forest resilience assessment methods (table 1). The length of the period used by investigators ranges from 1 year to 11 years (table S1), depending on investigators' research interests and focuses. Two or 3 year periods have been frequently used (table S2). This is because most studies concern forest resilience to drought events and drought legacy is significant within about three years (Gao et al 2020).

Disturbance and perturbation
Resilience is the capacity of a system to manage disturbance by resisting and recovering from change (Helfgott 2015). A disturbance is different from a perturbation (Rykiel 1985). Disturbance is a cause of perturbation that results from the interaction between an ecosystem and the disturbance. A disturbance is characterized by type, frequency, and intensity. A perturbation can be measured by a change in the state variable of an ecosystem. Different ecosystems might have different perturbations resulting from the same disturbance. In the face of the same disturbance, the ecosystem with smaller perturbation has stronger resistance (blue curve in figure 5), while the ecosystem with larger perturbation has lower resistance (red curve in figure 5).

Identifying extreme event
A perturbation is a dynamic fluctuation, defined by equation (4), i.e.
which can be negative or positive. We do not study all fluctuations, but only those negative enough to be above a certain level such that they are considered as an extreme event. Lloret et al (2011) used a criterion |x ′ /x| ⩾ 25% to define an extreme event. The variable they selected is averages of the annual BAIs in all trees V = BAI; event value at the time of peak impact is x = V dist and its pre-disturbance value is a 5 year averagē x = V pre (figure 5). Lloret et al (2011) also used summer Palmer drought severity index (PDSI), a disturbance index, to double-check if the event identified by BAI (perturbation) was consistent with occurrence of drought events. The criterion to identify a drought Vpre is the average value of a state variable in years preceding the disturbance, V dist is maximum perturbation, Venv is the value of the disturbed state during events, and Vpost is the average value of post-disturbance years. The averaging period used by different investigators ranged from 1 year to 11 years. The formulae used to calculate resistance, recovery, and resilience are listed in the right boxes. More details can be found in tables 1 and S1. event is usually a 50% or more growth reduction compared to previous years (e.g. Sanchez-Salguero et al 2018, 2013, Rubio-Cuadrado et al 2018b).

Resistance
Resistance is the ability of an ecosystem to persist during a disturbance (Nimmo et al 2015). Forest resistance to a disturbance can be characterized by minimum destruction of biomass or species removal, i.e. by a disturbance severity. Forest disturbance severity is the deviation of a state variable from its steady state, quantified by perturbation (Bhaskar et al 2018) or by their ratio where V pre and V dist are the value of the state variable pre-disturbance and at the time of maximum perturbation, respectively. Although definition (12) makes more sense in theory, most forest resilience scientists use definition (13) (e.g. Lloret et al 2011) to estimate forest resistance since it is dimensionless and therefore can be compared across different types of disturbance and ecosystem. A smaller absolute value of Resistance (13) or (12) indicates that forest resistance is stronger, and vice versa.

Recovery
Recovery is the ability of an ecosystem to return to undisturbed state, which can be defined as a vegetation recovery ratio, where V pre , V dist , and V post are the value of vegetation spectral index (SI) pre-disturbance, at the time of maximum perturbation, and at the time of assessment, respectively. The vegetation recovery index (14) was developed initially by Key and Benson (2006), and then modified by Chompuchan and Lin (2017), mainly for the assessment of recovery from fire disturbance. Many spectral indices exist and can be used as proxies of vegetation variables as listed in the where V post is an average of multiple post-disturbance years and V dist is maximum perturbation. As Lloret et al pointed out, this index is not valid when the system collapses during disturbance (e.g. no growth with zero tree-ring width). The vegetation proxies in recovery index definitions (14) and (15) are all sampled from disturbed locations and areas. However, particularly for fire disturbance, the difference of post-fire weather and pre-burn climate condition will also affect recovery.
To avoid the inaccuracies, Fornacca et al (2018) proposed an adjustable recovery index by comparing disturbed pixels with undisturbed pixels of the same vegetation type where, V ta ref is vegetation SI at selected reference undisturbed pixels of same vegetation type at assessment dates t a , while V ta post is the index at the disturbed location or area. t m is the time of maximum perturbation. Model (16) is more useful for fire disturbance The short-term perturbation pulses are averaged out. Sustained growth increases (%GC > 0, recovery) and declines (%GC < 0, perturbation) periods can be well identified. The perturbation event is not clearly defined, dependent on researchers. There is no clear definition for resistance and resilience. This approach can be used to scale up the biodiversity-resilience relationship with trait databases and remote sensing data.
The TIN can be used as a proxy of annual primary production. The TIN calculation smooths a yearly NDVI curve by removing extreme high/low values that might be false due to contamination by cloud cover or data processing procedures (Pettorelli et al 2005). The smooth or gap-filling procedures might mistreat good data as bad data. Vpost is basal area value at specific year after the hurricane disturbance.
This method is simple and straightforward, particularly for disturbance that is hurricane or forest fire. The indices are site-specific and not normalized to the pre-disturbance state, so it is difficult to compare different study areas or ecosystems, or even the same ecosystem with different state variables.
The following citations are additional papers that used one or more of the formulas in the above row (details in table S1 (available online at stacks.iop.org/ERL/16/053008/mmedia)). and assumes that the burned pixels have similar fluctuation dynamics to the reference unburned pixels.

Resilience
Resilience is the capacity to reach pre-disturbance steady states, and is estimated as the ratio of postdisturbance steady state V post to pre-disturbance steady state V pre , where steady states V pre and V post can be approximated by multiple-year averages of proxies of vegetation condition as discussed in section 4.2. Resilience index (17) converges to 1 for a full recovery. Generally, resistance and resilience of an ecological system have an inverse relationship to each other (Nimmo et al 2015, Matos et al 2019. The system performs resistance during disturbance and resilience following disturbance. For a given disturbance, the same resilience determined by resilience index (17) might correspond to many different recovery processes from the beginning to the end of a disturbance. For a case with high resistance and low resilience, disturbance impact is light but recovery is low, while for a case with low resistance and high resilience, the disturbance impact is heavy but recovery is fast. The justification for one resilience value corresponding to many intermediate recovering processes is that there is no disturbance impact term V dist in the definition (17).
To correct the problem, Lloret et al (2011) proposed a relative resilience index as, Thus, relative resilience can reflect the resilienceresistance inverse relationship, being small for the performance of low resilience with high resistance, and large for the performance of high resilience with low resistance. The linear inverse relationship between resilience and resistance does not always exist but a nonlinear one remains open to explore (Nimmo et al 2015).
To identify forest disturbance history (trend), Nowacki and Abrams (1997) developed a simple but useful BAI-based index named percentage growth change (%GC), defined as, where V pre = preceding 10 year tree-ring-width mean and V post = subsequent 10 year tree-ring-width mean. The 10 year radial-growth averaging technique smooths out short-term forest growth fluctuations related to climate factors. The turning point from negative to positive in (19) occurs at the second year following disturbance. A minimum growth increase ⩾ 50% has been widely used as a criterion to distinguish a disturbance (Wu et al 2014, Xu et al 2017. Major growth release or disturbance is usually identified by %GC ⩾ 100% (Nowacki and Abrams 1997). The index (19) opens an easy way to obtain a long-term record of forest disturbance and postdisturbance growth trend (recovery) using the longterm radial-growth averaging approach. The disadvantage of this approach is that it is difficult to obtain an exact mortality history from tree-ring chronologies (Xu et al 2017). However, it is robust to capture mega-drought disturbance history. To link the index (19) to forest resilience, Xu et al (2017) proposed a concept model to measure forest resilience by comparing growth regeneration and forest mortality against aridity severity.

Soil microorganism resilience
Soil is a nonrenewable life supporting source, containing ∼1500 Gt more carbon than the sum of carbon holdings in vegetation and the atmosphere (Crowther et al 2019). Soil is often disturbed by climate change-induced natural disasters, pollution, pesticides, and agricultural overdevelopment. To protect soil and maintain agricultural sustainability, there is a pressing need to predict soil biota system behavior in the face of such disturbances. Soil microorganisms are the key players that govern the input and output rate of soil organic matter pools, and their functioning and structure is essential for soil stability and sustainability (Cavicchioli et al 2019). Soil resilience scientists face a huge challenge in determining the resistance and resilience of soil biota systems to disturbance because microorganisms living in soil are so abundant and highly diverse, with estimates of up to 10 9 cells and 10 4 species per gram of soil, many of which are unknown (Griffiths and Philippot 2013). Therefore, soil microorganism resilience studies are different from forest resilience studies in state variable selection and approach design, often collecting soil samples in the field and incubating them in the laboratory. Indices used to analyze soil resistance and resilience are tested experimentally, using a control soil sample and a disturbed soil sample, respectively notated as C and P (figure 6). Artificial disturbances include simulated heat, fire, drought, and frost.

Soil state-variable selection
It is a long standing hypothesis that soil stability and function is governed mainly by microorganism species diversity (Botton et al 2006). Soil microorganism species richness has been widely used as a state variable to study soil resilience (Sousa 1980, MacGillivray and Grime 1995, Tilman 1996. The microbial biomasses C, N, and P can also be used as response variables to disturbance instead of counting

Disturbance treatments
Soil resilience studies are lab-based for purposes of mechanistic understanding, and eventually the results are translated into practice. There are many ways to produce different soil disturbance treatments (Napper et al 2009). Lab-based disturbance treatments include heat, fire, desiccation, dry-wet cycling, freeze-thaw cycling, nutrient limitation (C, N, P, S), heavy-metal addition, etc (Griffiths and Philippot 2013). The frequency, duration, and intensity of disturbance treatments depend on the investigators' research purposes (Shade et al 2012).

Resistance
The ability of a soil biota system to withstand a disturbance is tested experimentally by comparing an undisturbed control soil sample and an artificially disturbed soil sample, respectively notated as C and P (figure 6). A simple soil resistance index used by Kaufman (1982) is, where C 0 and P 0 are the values of a state variable for the undisturbed control soil and the disturbed soil, respectively, at the end of the disturbance. Although the soil resistance index (20) seems similar to the forest resistance index (13), the smaller absolute value of resistance (20) means that soil resistance is weaker, opposite to that of index (13). Some soil resilience scientists have used normalized change against control value as a soil resistance index (Sousa 1980, Biggs et al 1999, where D 0 = C 0 − P 0 is the difference between the control (C 0 ) and the disturbed soil (P 0 ) at the end of the disturbance (t 0 ). The disturbance effect usually refers to negative consequences such as reduction in biodiversity or biomass. However, a positive effect is possible for soils such as when a pulse of glucose is added to soil (Orwin and Wardle 2004). Thus, D 0 can be both positive and negative. Given 0 ⩽ P 0 ⩽ 2C 0 , soil resistance (21) is at a maximum (100%) as P 0 = 0 (no disturbance effect) and becomes lower as P 0 increases. To avoid negative values of resistance, Orwin and Wardle (2004) modified the index (21) to be, For 0 ⩽ P 0 ⩽ 2C 0 (i. e.|D 0 ⩽ C 0 |), the index (22) returns a value between 0 and 1. The smallest value 0 corresponds to the strongest resistance.

Resilience
Formulations of soil resilience indices are similar to resistance ones (20)-(22) due to being experimentally tested by similar undisturbed control soil samples and disturbed samples. The simplest soil resilience index is, where P x is the value of the state variable for the disturbed soil at subsequent time points (t x ) (figure 6) (Kaufman 1982). The index (23) indicates how much the disturbed soil has recovered with respect to the initial undisturbed soil (C 0 ). It is straightforward but there is no information as to the difference between undisturbed control soil and treated soil at the initial time (t 0 ). To this end, the following resilience index was used by (Sousa 1980 ). They compared resilience values between plots to deduce the effect of P and N fertilization on resistance and resilience.
Normalizes change against control value.
Can express positive or negative net change through the sign of the expression.
Provides a narrow view of resilience and resilience-does not consider recovery rate or time. (Continued.) where D x = C x − P x is the difference between the control and the disturbed soil at subsequent time points (t x ). The advantage of index (24) is that it does not take into account any changes that are not in response to disturbance. The index can be both positive and negative without bounded values. As D x → 0, the soil biota system approaches a full recovery. The disadvantage of index (24) is that it becomes infinite (→ ∞, singularity) if the system has a perfect resistance (D 0 = 0). To overcome this weakness, Orwin and Wardle (2004) proposed a new resilience index as, The value of this index is bounded between −1(|D x | → ∞) and +1 (|D x | = 0). D x = 0 indicates a full recovery. For a given |D x | ⩽ |D 0 |, the value of the index lies between 0 and 1. The value 1 indicates a full recovery at the measurement time (t x ), while the value 0 indicates that the disturbed soil has not recovered at all (i.e. D x = D 0 ). There is no universally accepted index for experimentally evaluating and quantifying soil resilience. All of the above indices of soil resistance and resilience are highly approximated from nature for soil scientists to use when testing their hypotheses by comparing undisturbed control soil samples and disturbed soil samples. The motivations of these studies are essentially mechanism understanding, which could translate into policy advice and improved land management practices (Griffiths et al 2013). More details concerning the applications of soil resistance and resilience indices are listed in table 2.

Hydrological resilience
Hydrological resilience focuses mainly on responses of catchment ecosystems to droughts or warming events. Its formulation is based on the Budyko framework (Turc 1954, Pike 1964, Budyko 1974, Zhang et al 2001, Yang et al 2008. State variables are not direct measurements of vegetation types and structure, but of a related micrometeorological variable: evapotranspiration (E). The central concept is a dryness index (DI) developed by Budyko (1974), which is a simple but thoughtful combination of water balance and energy balance. The Budyko curve, a function of the DI, is considered as a baseline of steady states to be compared with meteorological and hydrological data to examine the ability of catchment ecosystems to resist and absorb drought and warming related changes.

DI
Annual averages of ground heat fluxes and soil heat storage of an ecological system can be neglected in comparison with annual latent heat flux (LE) and sensible heat flux (H). Thus, annual net radiation (R n ) can be approximated as Potential evapotranspiration (E p ) is a very important parameter in ecohydrology as the maximum of E, and it can be taken as a limit H → 0 where L is a latent coefficient (= 2.5 MJ kg −1 , the enthalpy of vaporization) and R n is annual sum of net radiation (MJ m −2 yr −1 ). The equation (27) defines how the total water flux of an ecosystem is evaporated and transpired into the atmosphere by all available energy R n in a year. A DI is defined by Budyko as where P is annual sum of precipitation (mm yr −1 ). DI indicates a climatological water condition in a straightforward way: DI < 1 is wet, and DI > 1 is dry. This dimensionless index is a thoughtful combination of soil water inputs from the atmosphere by P in denominator and the energy driver of soil water outputs to the atmosphere by E p = R n /L in numerator. The DI contains essential climatological waterenergy support information in ecosystem production (e.g. NPP). Budyko (1974) used a two dimensional space of DI and R n to classify geobotanical zones (figure 7). Due to its simplicity, Budyko's method has been widely used to map vegetation patterns and study their ongoing relation to climate (e.g. Monserud et al 1993, Yanling et al 2008, Cai et al 2014, Constantinidou et al 2016, 2019, Trancoso et al 2016. The DI is applicable on annual or longer time scales but not for sub annual time scales (Yi et al 2012). Budyko (1974) formulated his framework in a two dimensional space spanned by evapotranspiration index EI = E/P and dryness index DI = E p /P. With it, hydrological scientists can test their hypotheses with meteorological and hydrological data based on physical principles, empirical models, and resistance and resilience concepts within certain temporal and spatial conditions (figure 8).

Budyko framework
The theoretical foundation of Budyko framework is the laws of mass and energy conservation. The energy governing equation is where C E is the effective heat capacity of a catchment ecosystem in J m −2 k −1 , T s is surface temperature, Here Rn (y-axis) plays a similar role as temperature. For example, with a given range of DI, the forest type shifts from tropical forest in low latitudes to coniferous forest in high latitudes as Rn decreases. (1) P > Ep, so evapotranspiration is limited by available energy; and (2) P < Ep, so evapotranspiration is limited by available water. Here, EI = E/P and DI = D/P, and all curves are generated by the non-parameter models (Schreiber, Ol'dekop, Budyko, and Pike) and parameter models (Fu, and CY →Choudhury and Yang et al) listed in table 3. and G is ground heat flux. The soil water governing equation is where S w is soil water storage (mm t −1 ), and Q is streamflow (mm t −1 ). The Budyko framework is limited to steady-state conditions, which can be approximated if all terms in equations (29) and (30) are taken as long-term annual averages This is because changes in long-term averages of water and heat storage are negligible in comparison with short time average; for example, ∂T monthly s /∂t ≫ ∂T annual s /∂t and dT monthly s /dt ≫ dT annual s /dt. Thus, the Budyko framework is not valid for temporal scales shorter than annual averages (Youreka et al 2019). In addition, the Budyko framework addresses the water and energy exchanges of a catchment ecosystem with the atmosphere and assumes groundwater inputs and ground heat fluxes are negligible, particularly for ecosystems over larger spatial scales. Thus, the energy balance equation (31) is reduced to (26).
In Budyko space (figure 8), the data of R n and P of four mean annual variables (E p (= R n /L), P, E, Q) are well collected and hence abundant, but E and Q are difficult to measure. The Budyko framework provides a simple way to partition P into E and Q. Two limits in wet (DI < 1) and dry (DI > 1) region are given by two extreme cases with energy and water availability: (a) P > E p , i.e. evapotranspiration E is limited by available energy R n , so EI = DI and runoff occurs; and (b) P < E p , i.e. evapotranspiration E is limited by available water, so EI = 1 and no runoff occurs. Therefore, Table 3. Budyko-curve family.

Equations
References Pike (1964), Turc (1954) Yang et al (2008) EI = E/P, DI = E0/P, m = parameter, n = parameter all long-term water balance data should logically fall below the two limits. What function of EI with DI is to follow based on observational data? We can obtain it by dividing both sides of the energy balance equation (26) by LE, where B r = H/(LE) is a Bowen ratio that is a function of DI (Arora 2002). Thus, we obtain which forms a physical basis of all possible relationships of EI with DI and is called the Budyko hypothesis.
Following the Budyko hypothesis, many empirical equations (table 3) have been developed by the best-fit curves (figure 8) through mean annual climate data. All these empirical fitting curves can be viewed as a Budyko-curve family. The first two empirical equations were developed by Schreiber (1904), and by Ol'dekop (1911) Budyko (1948) found that most hydrological data fell between the predicting curves of equations (35) and (36), so he used the geometric mean of the two equations as a new equation Relationship (37) is a benchmark curve in Budyko space (bold solid black curve in figure 8). In the Budyko curve family, like equations (35)-(37) are called non-parametric models (table 3). The runoff fraction can be inferred as: (a) DI − EI if P > E p ; and (b) 1 − EI if P < E p . The steady-state condition is highly requested by these non-parametric models because there are no adjustable parameters in these models, which become more valid with increasing temporal and spatial scales. As temporal and spatial scales decrease, site-specific characteristics, such as vegetation, topography, soil properties, and even underground water, significantly contribute to a catchment water balance.
To address the site-specific departures from the Budyko curve, many researchers have developed parameter models, of which two are included in table 3 and more of which can be found from the table 1 in Mianabadi et al (2020). Fu (1981) used Buchingham's PI theorem and the assumption that the slope ∂E/∂P is a function of residual evapotranspiration E p − E and precipitation P to derive the well-known Fu's equation, where n is an adjustable parameter, its smaller values correspond to conditions with less soil-water holding capacity and vice versa.

Budyko resilience metrics
Based on the Budyko framework, Creed et al (2014) formulated three indices: static deviation (s), dynamic deviation (d), and elasticity (e), to quantify the resilience of catchment water yield to climate change.
The Budyko curve is considered as a long-term steady-state dimension, which can be obtained from the Budyko-curve family by best-fitting long-term data for a specific catchment. The specific values of parameters such as m or n in the Budyko models (table 3) will determine the position of a catchment system with a specific forest type and set of landscape characteristics in the Budyko curve family (figure 8, Goyal 2017, Sinha et al 2018). The three indices (figure 9) can be used to measure the deviations of EI relative to the Budyko curve so as to assess a catchment's resilience and elasticity by comparing a cool period with a warm period.

Cool and warm periods
Before calculating Budyko resilience metrics, the cool and warm periods are distinguished by calculating  3-or 5-water-year (5-wyr) moving averages of the catchment temperature (Creed et al 2014, Sinha et al 2018). The cool period is defined as the interval with the minimum 5-wyr temperature, while the warm period is defined as a 5-wyr period after the cool period that was warmer than subsequent 5-wyr periods by more than 1 standard deviation ( figure 9(a)). The length of the moving average period depends on researchers' purpose.

Deviation
Deviation refers to a vertical departure of measured EI from the theoretical Budyko curve (figure 9(a)). Static deviation (s) is calculated by from the cool-period observations. EI M,cool is the measured EI during the cool period and EI B,cool is from the Budyko curve. To assess a net change in a catchment's EI relative to Budyko curve during the warm period, Creed et al (2014) assumed static deviation s to be constant with time. Thus, dynamic deviation d is calculated by where EI M,warm and EI B,warm are the EI measured from observations and predicted from the Budyko curve respectively. Data points above the Budyko curve demonstrate smaller-than-expected water yield, while data points below the Budyko curve indicate larger-than-expected water yield.

Elasticity
Creed et al (2014) used the concept of 'Elasticity' to characterize a catchment's ability to maintain water partitioning P into EI and Q, and to test if the catchment stays close to or departs far from the expected Budyko curve. The elasticity (s) is measured by where DI max − DI min is the DI range between the cool and warm period, EI R,max = ER max,warm − EI B,warm is the maximum EI residual value during the warm period, and EI R,min = |EI M,cool − EI B,cool | is the maximum absolute EI residual value during the cool period (figure 9(b)). A highly resilient catchment has a high degree of elasticity (e > 1) and hence its EI data points during the warm period are close to the Budyko curve, while a nonresilient catchment exhibits a low degree of elasticity (e < 1) and its data deviates from the Budyko curve. e = 1 is the threshold value for elastic vs. inelastic catchments. The elasticity (resilience) of a catchment is highly dependent on forest type and warm period selection. The Budyko framework and the elasticity metrics have become useful tools for forest managers to use when assessing the resiliency of catchment water yields to climate warming (

Concluding remarks
The concepts of resilience and tipping point have become more widely-known terms for people to understand socioecological systems' ability to deal with sudden environmental changes. This is evidenced by the fact that I received 185 000 000 results by using Google Search for the word 'Resilience' on 26 September 2020. The theoretical understanding of steady-state stability and bifurcation point is rooted in physics and mathematics and the dynamic models associated with differential equations. The two conceptual systems hold the same aims but both are impractical in the real world. Holling (1973) intended to translate hard-to-understand mathematical concepts of steady-state, fluctuation, stability, and bifurcations into easy-to-follow concepts of resistance, recovery, resilience, and tipping point to pioneer a way for most people to understand sudden change. However, it is still challenging to calculate resilience by a formula or measure it by a universal method.
Hopefully, this review has further clarified the links between the two conceptual systems and summarized the state-of-art methods of measuring resilience for terrestrial ecosystems, which continue to be developed by resilience scientists. At present, there is no single measurable variable that can serve as a state variable for system resilience study. Treering width, NDVI, microbiome mass, and catchment evapotranspiration index can be used as indicators of system state variable, but not as actual state variables. In this case, it is hard to say which is better-it is dependent on investigator research objectives and data availability. A steady state can be approximated by the multiple-year-average of proxy of system state variable. Recovery is a decay process of fluctuation. A highly resistant system would have small fluctuations with a lower recovery and less room for resilience. The Budyko framework and resilience metrics are empirical but practical to study the resilience of catchment yields. The Budyko curve obtained by best-fitting long-term datasets is considered as a long-term steady-state, and the static and dynamic departures are the fluctuations along the steady-state (Budyko curve). Elasticity is also a relative fluctuation. Whether a catchment is resilient or non-resilient can be determined by comparing the departure and relative departure (elasticity) behavior-whether they stay close along the steadystate or deviate far away between cool and warm periods.
Dynamic system theory is a fundamental base of resilience science. Uncertainties in understanding resilience are rooted in the complexities of a nonlinear system, which consists of many positive and negative feedback loops. Positive feedbacks are driving forces for perturbation growth, while negative feedbacks are driving forces for perturbation decay. Resistance, recovery, and resilience are the results of competition and cooperation between positive and negative feedback loops. If policy makers understand the structure of the feedbacks of a nonlinear socialecological system, they are able to manage the feedback loops to reduce the perturbation or speed up the recovery, or to avoid tipping into a new stable steady-state. This review is a starting point of such an endeavor, linking practical resilience indicators to some concepts of stable steady-state in dynamic stability theory but not linking them to feedback loops, regime shifts or tipping points. Further studies are needed: (a) to link practical resilience calculations to feedback loops; (b) to identify which feedback loops govern the others; and (c) to find out how to manage the governing feedback to avoid catastrophic disasters.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Author contribution
C Y conceived and conducted the analysis, and wrote the manuscript. N J conducted literature search and contributed section 3, figure 4 and all tables except table 3 and table S2.

Acknowledgments
We thank Guillaume Wright and editors of Environmental Research Letter for inviting this review. We also thank three anonymous reviewers for providing helpful critique that resulted in substantial improvements.