Quantifying the risk of irreversible degradation for ecosystems: A probabilistic method based on Bayesian inference

Ecological Indicators 107 (2019) 105621

Contents lists available at ScienceDirect

Ecological Indicators journal homepage: www.elsevier.com/locate/ecolind

Quantifying the risk of irreversible degradation for ecosystems: A probabilistic method based on Bayesian inference Sifeng Wua,b, Zhongyao Liangc, Yong Liuc,

T

⁎

a

Complexity Institute, Interdisciplinary Graduate School, Nanyang Technological University, Singapore Complexity Institute, Nanyang Technological University, Singapore c College of Environmental Science and Engineering, State Environmental Protection Key Laboratory of All Materials Flux in Rivers, Peking University, Beijing 100871, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Irreversible degradation Ecosystem resilience Bayesian inference Spatial heterogeneity Disturbance Ecosystem-based management

Ecosystem degradation is usually abrupt and unexpected shifts in ecosystem states that cannot be easily reversed. Some ecosystems might be subject to high risks of irreversible degradation (RID) because of strong undesirable resilience. In this study, we propose a probabilistic method to quantify RID by measuring the probability of the recovering threshold being unattainable under real world scenarios. Bayesian inference was used for parameter estimations and the posteriors were used to calculate the threshold for recovery and thereby the probability of it being unattainable, i.e., RID. We applied this method to lake eutrophication as an example. Our case study supported our hypothesis that ecosystems could be subject to high RID, as shown by the lake having a RID of 72% at the whole lake level. Spatial heterogeneity of RID was significant and certain regions were more susceptible to irreversible degradation, whereas others had higher chances of recovery. This spatial heterogeneity provides opportunities for mitigation because targeting regions with lower RID is more effective. We also found that pulse disturbances and ecosystem-based solutions had positive influences on lowering the RID. Pulse disturbances had the most significant influence on regions with higher RID, while ecosystem-based solutions performed best for regions with moderate RID, reducing RID to almost 0. Our method provides a practical framework to identify sensitive regions for conservation as well as opportunities for mitigation, which is applicable to a wide range of ecosystems. Our findings highlighted the worst scenario of irreversible degradation by providing a quantitative measure of the risk, thus raising further requirements and challenges for sustainability.

1. Introduction Degradation of ecosystems is usually unexpected and abrupt shifts in ecosystem state that cannot be easily reversed (Scheffer et al., 2001; Folke et al., 2004; Ratajczak et al., 2018). The degraded state can be stabilized by the internal reinforcing feedback of the ecosystem and therefore is resistant to mitigation (Walker et al., 2012; Folke et al., 2010). Worse scenarios are where the degradation becomes irreversible if the internal reinforcing feedback is too strong to break (Carpenter et al., 1999). Most studies have focused on the amount of disturbance that will cause the shift of an ecosystem to an undesirable regime and the threshold associated with degradation (Yeung and Richardson, 2016; Walker et al., 2002; Cumming and Peterson, 2017; Andersen et al., 2009), whereas uncertainty and the risk of recovery threshold that is unattainable have received limited attention. Nevertheless, this risk is an extremely important concern as it determines whether the

⁎

degradation is reversible or not(Scheffer et al., 2001; Carpenter et al., 1999). Under the context of global ecosystem degradation (Cumming and von Cramon-Taubadel, 2018; Steffen et al., 2015), this risk will have a profound influence on management perspectives. In theory, when ecosystems lose their key species or functional group, they can lose the ability of reorganization and evolving into the desirable regime; thus, they are subject to the risk of irreversible degradation (Folke et al., 2004; Elmqvist et al., 2003). Our efforts in the restoration of particular ecosystems have shown that most ecosystem degradation, if not irreversible, is extremely difficult to reverse (Nyström et al., 2012; Cumming and von Cramon-Taubadel, 2018; Scheffer and Carpenter, 2003). Previous studies have argued that internal reinforcing feedback and complex cross-scale interactions are the main factors hindering the recovery process (Oliver et al., 2017; Nyström et al., 2012; Steffen et al., 2015; Cumming and von CramonTaubadel, 2018). We further hypothesized an extreme scenario where

Corresponding author. E-mail address: [email protected] (Y. Liu).

https://doi.org/10.1016/j.ecolind.2019.105621 Received 1 April 2019; Received in revised form 25 June 2019; Accepted 4 August 2019 Available online 14 August 2019 1470-160X/ © 2019 Elsevier Ltd. All rights reserved.

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

parts; therefore, certain parts could be more susceptible to irreversible degradation, while other parts could be more responsive to mitigation (Janssen et al., 2017). Identification of the more vulnerable parts could facilitate conservation priority before degradation and leveraging mitigation effects after degradation (Chapin et al., 2010; Folke et al., 2004; Scheffer and Carpenter, 2003). Minimal models, although widely applied in regime shift studies, are often flawed by their homogeneity assumption because the main motivation of these models is to capture the nonlinear dynamics and dominating mechanisms that produce regime shifts (Allen and Hoekstra, 2015; Filatova et al., 2016). This can be overcome by adding interaction terms to the original model functions, which reflect spatial dynamics as well (Ratajczak et al., 2017; Kéfi et al., 2011; Tilman et al., 1994). Therefore, we first modified the original model by adding exchange terms and then explored spatial heterogeneity in undesirable resilience for each part of the system. Irreversible degradation, as defined in the present study, is with respect to long-term external drivers, also known as press (Angeler and Allen, 2016; Ratajczak et al., 2017). However, regime shift could also be triggered by pulse disturbances. Strong and short-term disturbance could lead to huge fluctuations in the system state (Ratajczak et al., 2017; Walker et al., 2012). Given strong enough disturbances, ecosystem states could be driven across the boundary of the alternative stable state (the unstable equilibrium) and then gradually converge to the desirable state (Walker et al., 2012; Ratajczak et al., 2017). Therefore, it is important to incorporate the stochastic environments to better characterize the risk of irreversible degradation (Ratajczak et al., 2018; Carpenter and Brock, 2006). To be consistent that the risk is associated with the ecosystem itself, the stochastic process was associated with the internal process, i.e., sediment release. Both long-term lP reduction as press and random noises as pulse were implemented simultaneously in our simulation. Backward shift from degraded regime to the desirable regime was then compared with the stationary RID calculated by thresholds. Increasing research has argued that ecosystem-based management and nature-based solutions are better methods for conservation and protection (Keesstra et al., 2018; Cohen-Shacham et al., 2016; Sasaki et al., 2015). Other studies have remained skeptical regarding their performance. We testified the usefulness of ecosystem recovery on ecosystems with high RID by testing the effects of changing parameter values on RID. The results provide hints on the performance of ecosystem recovery as it targets ecosystems and attempts to change ecosystem attributes. Moreover, by doing this, we could identify the ecosystem attributes that were most efficient for lowering RID to leverage mitigation effects.

ecosystems could be subject to a high risk of irreversible degradation, making the recovery extremely difficult or even impossible. This worst case scenario which leads to the most severe consequence received much less attention than it should have. This risk is usually characterized indirectly by general resilience indicators, for example, diversity and redundancy (Walker et al., 2012; Oliver et al., 2017; Scheffer et al., 2015). We still lack a concrete indicator to quantify this risk. To verify our hypothesis and quantitatively measure this risk, we proposed a probabilistic indicator, RID, to quantify the risk of irreversible degradation. It is defined as the probability of the threshold for recovery (T2 ) being unattainable under real world scenarios. This could result from strong undesirable resilience where the ecosystems are firmly stabilized by the internal reinforcing feedback and are thereby trapped in the undesirable regime (Holling, 1973; Ibelings et al., 2007; Walker et al., 2012). Traditional resilience indicators, following the original definition of resilience, are usually representations of how much change the system can withstand before flipping into another regime (May, 1977; Walker et al., 2012; Yeung and Richardson, 2016; Scheffer et al., 2015; Carpenter et al., 2001). However, these methods are less applicable to systems that are subject to high risks of irreversible shift because it is illogical to quantify resilience based on an unrealistic threshold that might never be attained. A highlight of this indicator is the extreme scenario where the degradation becomes irreversible, i.e., the worst possible consequence. Moreover, rather than measuring uncertainty associated with external environments that have been commonly observed in previous studies (Rougé et al., 2013; Hazbavi et al., 2018; Carpenter et al., 2012; Ratajczak et al., 2017), RID is measured as an intrinsic property of ecosystems. This is realized by the Bayesian inference that allows for a dynamic view regarding parameterization (Ellison, 2004). Different from the frequentist approach where parameter values are treated as a fixed value, Bayesian inference treats parameters as random variables that vary with certain probabilities (Ellison, 2004; Arhonditsis et al., 2006). Therefore, the uncertainty arises from the parameters that are representations of internal processes and properties of the ecosystem (Wu et al., 2017; Liang et al., 2018; Arhonditsis et al., 2006). Thus, the risk measured in Bayesian approaches also arises from parameters, suggesting that the risks are associated with the intrinsic properties of the ecosystem. Another advantage of the Bayesian method is that we can make full use of the knowledge by incorporating prior information (Ellison, 2004). When given new data, the previous posterior data can then be used as a prior that will be updated based on the new information. This progressive process suits the purpose of adaptive management and could therefore be a promising tool for management practices (Dorazio and Johnson, 2003; McCarthy and Possingham, 2007). We applied our method to the eutrophication of freshwater ecosystems driven by increasing phosphorous (P) inputs. Accumulating empirical evidence has shown that internal sediment release alone contributes much more than external inputs (Wu et al., 2017; GenkaiKato and Carpenter, 2005). This suggests the possibility of irreversible eutrophication, and therefore matches with our study objective. We first attempted to quantify the undesirable resilience as the risk of irreversible degradation and then investigated the influence of several aspects that might provide opportunity for mitigation. We included the following aspects: i) spatial heterogeneity, ii) interactive effects of pulse disturbances and long-term press, represented by random noise and long-term P load (lP ) reduction, respectively (Ratajczak et al., 2017; Hughes et al., 2013); and iii) changes in internal ecosystem attributes resulting from ecosystem recovery, reflected by the overall changes in parameter values. Investigation of these aspects revealed opportunities to reduce RID for ecosystems that are subject to high RID. Spatial heterogeneity and connectivity between different parts of the ecosystems are an important source of ecosystem resilience (Holling, 1973; Oliver et al., 2015). Because of spatial heterogeneity, RID could also vary among different

2. Methodology The indicator was defined by the probability of T2 being unrealistic. We often do not know the real distribution of T2 ; therefore, we can calculate the empirical distribution, the frequency based on large random sampling, as an approximation. The calculation of RID can then be expressed as Table 1

RID = P(T2 < T2c ) =

n (T2 < T2c ) n (T2)

(1)

where, T2c represents the critical threshold for environmental pressure, which will be unattainable in real world scenarios; P(T2 < T2c ) represents the probability for T2 below T2c ; and n (T2 < T2c ) and n (T2) represent the number of samples T2 < T2c and the total sample size, Table 1 Risk of irreversible degradation at the whole lake level and for each subunit.

2

Whole lake

Subunit 1

Subunit 2

Subunit 3

Subunit 4

0.716

0.173

0.864

0.957

0.411

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

resulting in unnecessary complication. In the present study, the lake was divided into four subunits based on proximity polygons, also referred to as nearest neighbor interpolation. All points were divided into the region that contained its nearest point. The number of subunits was decided by the availability of water quality data and the corresponding monitoring site. Results of the spatial division could be found in the Supplementary material (Fig. S.1).

Bayesian inference for the focal system — parameter posterior Large parameter sample drawn from posterior T2 distribution: empirical distribution calculated from parameter sample

2.2. Solving differential equations for threshold

RID: risk of irreversible degradation

Following our rationale on irreversible degradation, to quantify RID, we had to calculate the threshold for recovery (T2 ). T2 is the minimal lP value for the degraded regime. As regimes are defined by a fixed point, they could be obtained by setting Eq. (2) to 0:

Fig. 1. Procedures for calculating risk of irreversible degradation.

lP∗ =

lP∗

dlP∗ d ⎛ rP ∗q ⎞ =s− =0 ∗ dP dP ∗ ⎝ P ∗q + mq ⎠

By solving Eq. (6), we obtained at the threshold, and then lP∗ threshold was obtained by plugging in P ∗ to Eq. (5). For illustration purpose, we usually plot P ∗ against lP∗, indicating that lP is the control variable and P is the response variable. However, to be written as a function, we can only write lP∗ as a function of P ∗ because each lP∗ can map to more than one P ∗ when alternative stable states exist. Similarly, the threshold for each subunit was found by solving Eq. (7):

We selected the classic P driven eutrophication model for our study because it is a good representation of the ecosystem and illustrates our method well. The model for the whole lake followed Eq. 2 (Carpenter et al., 1999).

dP = dt

n

∑ i=1

wi

n

eij Pi +

j=1

dPi = lP − dt

∑ j=1

wj

p (θ|X) =

wi

n

∑ i=1

q Pi i

Pi i q + mi i

q Pi i

Pi i q + mi i

n

wi si Pi −

∑ i=1

(3)

q

wi ri

j=1

wj wi

n

∑ j=1

∗q

⎞ d ⎛ ri Pi i ⎞ eji + ⎜ ∗q q ⎟ = 0 ⎟ dPi∗ ⎝ Pi i + mi i ⎠ ⎠

(7)

p (X|θ) p (θ) p (X)

(8)

where, X and θ represent observations and parameters, respectively. p (θ|X) denotes the posterior, the probability distribution of θ given X;p (X|θ) denotes the likelihood function, the probability of X given parameter θ ; and p (θ) denotes the prior assigned to θ before the inference. p (X) denotes the marginalized probability of X , which will be a constant. For our process-based model that involves multiple parameters and noisy observations, the inference was conducted according to Eq. 9 (Hartig et al., 2017; Bishop, 2006):

q

eji Pj + ri

eij +

The parameter posterior was calculated based on Bayes’ theorem (Eq. (8)).

(2)

n

n

∑

2.3. Parameter inference and risk of irreversible degradation

where, P represents the state variable, water P concentration; lP represents the external P load, r represents the potential sediment release; m is the half-saturation coefficient where the actual P release reached half of its potential, and q is the shape parameter that controls the steepness of response between actual sediment P release and P. For the estimation of each subunit within the lake, we added an exchanging term between each unit and obtained Eq. 3:

∑

(6)

P∗

∗ dlPi ⎛ = −si − ⎜ dPi∗ ⎝

dP Pq = lP − sP + r q dt P + mq

(5)

P∗

where, and denote lP and P at equilibrium, respectively. Then, thresholds, as local maximum or minimum, were found by taking the derivative with respect to P ∗:

2.1. Model description

dPi = lPi − si Pi − dt

rP∗q =0 P∗q + mq ∗q rP sP ∗ − P∗q + mq

lP∗ − sP ∗ +

respectively. We directly used the probability of T2 < T2c as an indicator of risk so that the expression was more compact. This implied that the unstated loss function should be L =  { < } ( ) , mapping all T2 < T2c to 1 and the rest to 0. This is one of the typical loss functions. Our rationale for using this binary loss function was that we are only interested in whether the degradation will be irreversible or not. The risk function, as the expectation of the loss function defined above, became the probability written as Eq. 1. For our case of lake eutrophication driven by increasing lP , we set T2c = 0 because it is unrealistic to reduce lP to below 0. To calculate RID, there were two steps: i) Bayesian inference for the posteriors of parameter distribution and ii) calculation of the threshold based on a parameter sample drawn from posteriors of Bayesian inference. We also undertook several dynamic simulations to investigate how factors, including spatial heterogeneity, pulse disturbances, and ecosystem recovery, influenced RID (Fig. 1). In the following part, we elaborate on each step of our method.

p (θ|y , X , σ , α ) ∝ p (y|X , θ , σ ) p (θ|α ) (4)

(9)

N

p (y|X , θ , σ ) =

where, all the subscripts i denote the ith subunit; eij represents the outflow rate from subunit i to its neighboring subunit j and eji the outflow rate from subunit j to i; and wi represents the weight assigned to each subunit based on their water volume (water volume of subunit i over the total water volume of the lake). The second equation acted as a regularization term to ensure mass balance such that the sum of flows in all subunits will not deviate greatly from that at the whole lake level. Note that for the inference at the whole lake level, we should adopt Eq. (2) rather than this regularization term because it involves more parameters than are necessary for the inference at the whole lake level,

∏

N (y|f (X , θ), σ )

n=1

where, y represents the observation and X the input; then p (y|X , θ , σ ) was the likelihood function measuring the likelihood of observations under the given input, model parameters and another parameter σ . This followed a normal distribution with mean f (X , θ) , and variance σ 2 . The rationale behind undertaking this was that the observations always involve noise following a normal distribution. The function f (X , θ) here should be consistent with Eq. (2) or (3) according to the object of the inference. In addition, the p (θ|α ) is the prior set with parameters α , 3

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

2.4.2. RID under ecosystem recovery Another type of mitigation measurement was targeted at the ecosystem itself rather than at reducing the external environmental pressures. We attempted to reveal the effect of ecosystem recovery on RID. Ecosystem-based recoveries lead to changes in ecosystem attributes. Taking the lake ecosystems as an example, a typical ecosystem-based management would be the removal of sediment to decrease sediment release (Zamparas and Zacharias, 2014). This effect will be equivalent to lower values in the sediment release rate r in Eq. (2) (Carpenter et al., 1999; Carpenter, 2005). Therefore, overall changes in parameter values should be a reasonable representation for the effects of ecosystem-based management. Considering that the spatial and temporal variation of ecosystems and thereby the parameter values will not diminish (Beisner et al., 2006), we retained the distribution and variation in the parameter sample but increased or decreased proportionately. RID was calculated for changes in different parameter values. Each time, one of the parameters was set proportional to its original values in the sample while the other parameters remained unchanged. By comparing RID under all these variations, we concluded the usefulness of ecosystem recovery for systems subject to high undesirable resilience, as well as the performance of targeting different attributes of the systems in different parts. In this paper, all simulations were performed with R software (R Core Team, 2017). Bayesian inference was implemented by package BayesianTools (Hartig et al., 2017). Root finding was conducted by package rootSolve (Soetaert, 2009). Spatial division was conducted by package dismo (Hijmans et al., 2017). Data for parameter inference and related simulations (including lake morphology, water quality, and lP ) was based on Lake Erhai, one of the largest lakes in Southeast China. Details of this lake can be found in the Supplementary material (Fig. S.1). The results of parameter inference can be found in the Supplementary material, Fig. S.2 for parameter posteriors and Figs. S.3 and S.4 for model performance at the whole lake level and subunit level, respectively.

which in our case would be the upper and lower bounds of the uniform distributions. We could not obtain a closed-form solution for parameter distributions; therefore, we used a Markov Chain Monte Carlo simulation to obtain an approximation instead (Qian et al., 2003). We had to establish the prior and the Markov Chain for the inference. For the priors, considering that great variation exists among different lakes and no specific information on the mean and variance of parameters for our lake was available, we adopted uniform distributions as priors. The upper and lower bounds for our uniform priors were given such that they covered the range that we found in the literature. The posterior was simulated with the differential evolution Markov Chain MonteCarlo algorithm (DEMCzs), a differential evolution method (Hartig et al., 2017). Convergence was checked using Gelman-Rubin diagnostic (Gelman et al., 2014). We computed two chains with the DEMCzs, each with 500000 iterations, and our model converged on all parameters (Gelman-Rubin’s test of convergence below 1.05 for all parameters). The parameter sample was drawn from the posterior. The first half of the chains was discarded to ensure that we only sampled after the convergence. The size of our parameter sample was 5000. After obtaining the parameter sample, we calculated the values of T2 under each parameter combination. This gave us a set of T2 with size 5000, equal to the size of the parameter sample. Because our resilience indicator was defined as the risk of irreversible degradation, it was easily obtained by calculating the probability of T2 < 0 (Eq. (1)). 2.4. Simulation under dynamic scenarios The calculation of RID was a reflection of the intrinsic characteristic of the lake, indicating the possibility of transition back to the desirable regime after degradation. This theoretical result could be influenced by processes and changes at different scales. Thus, to ensure our results were a better reflection of the real world as well as to identify possible mitigation opportunities, we conducted simulations under several dynamic scenarios in addition to the calculation of the threshold.

3. Results 2.4.1. Interactive effect of pulse and press on RID To simulate the interactive effect of pulse and press on regime shift, in addition to load reduction as a typical long-term press, we added a white-noise term to represent pulses, then the model function was written as Eq. (10).

dP = lPi − si Pi − dt

n

∑ j=1

n

eij Pi +

∑ j=1

wj wi

3.1. RID at different spatial scales The probability distribution of T2 for the whole lake resembled a normal distribution with some skewness (upper plot in Fig. 2). The value of T2 ranged from −0.10 to 0.05 mgP·m−3·month−1 and the probability of T2 < 0 was high. This is better illustrated by the cumulative probability for T2 (lower plot in Fig. 2). RID, according to our definition, was the cumulative probability at T2 = 0 . At the whole lake level, RID was 71.6%. This result showed that Lake Erhai suffered from a high risk of irreversible degradation. Therefore, it is a good study object in the present study. Calculation for each subunit showed differentiated RID among the four subunits. RID for subunits 1 and 4 were substantially lower than the other two subunits (Fig. 2). Subunit 1 had the lowest RID of 16.9%. Subunit 4 had a relatively low risk as well, with a RID of 40.9%, whereas RID for subunits 2 and 3 were significantly higher, 86.8% and 95.3%, respectively (Table 1). These results supported our hypothesis on spatial heterogeneity of RID. Subunits 1 and 4 had a RID that was lower than that of the whole lake, whereas subunits 2 and 3 had a RID that was higher than that of the whole lake (Fig. 2). This indicated higher chances of recovery for subunits 1 and 4, whereas subunits 2 and 3 would be more resistant to mitigation. The reason for differentiated spatial RID, for our case here, could be attributed to the hydrology conditions, which could be seen from the parameter posteriors for inflow/outflow (s) and exchange rate (e). The inflow and outflow rates for subunits 1 and 4 were much higher than those for subunits 2 and 3. Better hydrodynamic conditions for subunits 1 and 4 led to lower RID and less undesirable resilience. Another possible reason, although not directly reflected in our data, could be the

q

eji Pj + ri

q Pi i

Pi i q + σdWσdW + mi i (10)

where, σdWσdW denoted the white noise. σ reflected the strength of the random disturbances. The noise was assumed to obey a normal distribution because noise distributions are well described by normal distributions in a wide range of contexts (Ratajczak et al., 2017; Carpenter and Brock, 2006). Because a random noise was added to the equation, we could not obtain the equilibrium by directly solving the equations. Rather, the stable states were obtained by running the model for long iterations until the final state became steady. The initial P was set to higher than that at T2 (1.2 × PT2 ) to ensure that the system started from a degraded regime. The reduction rate was set proportional to T1 (from 100% to 0), the threshold for degradation. This could not guarantee a recovery but being lower than T1 could make the system become bistable and therefore shift the regime from a eutrophic regime to a clear regime. After long iterations, if P at stable states was lower than that of T2 , we considered the system to be recovered. Calculation under the same parameter sample could give the probability of recovery under the dynamic scenarios accordingly. By comparing this result with previous RID, we obtained further insights on the recovery process for systems with high undesirable resilience. 4

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

count

600

400

200

0

1.00

−0.05

0.00

0.05

−0.05

0.00

0.05

spatial scale wholelake

Cumulative Probability

subunit1 0.75

subunit2 subunit3 subunit4

0.50

0.25

0.00

Threshold for recovery ( T2 ) : mgP / m3 / month Fig. 2. Risk of irreversible degradation measured at different spatial scales. The x axis represents T2 , the threshold for recovery. The upper panel shows the histogram of T2 at the whole lake level; the lower panel shows the cumulative probability of T2 for both whole lake level and different subunits. Cumulative probabilities at T2 = 0 equal the RID as T2c for our case here equals 0.

intensity required to trigger the shift were decided by the difference between the current system state and the alternative stable state boundary. The higher σ values, the higher probability of disturbance intensity that lived up to this requirement, and therefore a higher chance of recovery. Certain discrepancy existed between the recovery rate and that indicated by RID. By definition, the recovery rates of reduction rate equal to 1 and σ equal to 0 were supposed to be 1 − RID ; however, they were somewhat different (for a graphical illustration, please refer to Fig. S.5 in the Supplementary materials). This was caused by a different final P for all subunits. To calculate RID, we assumed that all subunits would have the same P under equilibrium. However, this is not entirely true considering all the processes occurring at different rates in the subunits. Subunits with higher RID will have lower recovery rates, while subunits with lower RID will be more likely to recover. This consistency suggested that RID was a good representation of recovery probability even under dynamic simulation.

relatively larger shore area (please refer to the Supplementary materialfor an illustration). Because shore area could play an important role in the reduction of nutrients and recovery of eutrophication, it is possible that subunits 1 and 4 could have a better chance of recovery owing to these having a higher portion of the shore area.

3.2. RID under interactive effect of pulse and press disturbances For all subunits, as the reduction rate increased, the recovery rate increased as well (Fig. 3). Subunit 1, with a RID much lower than the other subunits, displayed a distinctive response pattern. It responded relatively proportional to load reduction in an approximately linear way. The other subunits were much less responsive at low reduction rates, and only became more responsive at high reduction rates. The addition of pulse disturbances (σ values higher than 0) imposed a positive influence on the recovery rate. For all subunits under all reduction rates, pulse disturbances led to higher recovery rates. This positive influence was more significant for subunits with higher RID, as the recovery rates of subunits 2, 3, and 4 differentiated more significantly under different σ values. This could be because subunits with low RID were already mostly reversible even without external disturbances, whereas subunits with high RID were more reliant on external disturbances to provide the opportunities to trigger the backward regime shift. Without external disturbances to drive their state beyond the alternative state boundary, these subunits will have much less chances of recovery. This also explained why this positive influence became more significant with more intense disturbances (increasing σ ). The disturbance

3.3. Influence of ecosystem restoration Varying parameter values could influence RID; however, the effect of different parameters differed in different subunits (Fig. 4). For all subunits, m and r were the most influential parameters for RID. A decrease in r and an increase in m could lower RID. This influence was consistent and obvious in all subunits. Both r and m are relevant to sediment release, r being the potential for sediment release and m being the half saturation coefficient. A decrease in r and an increase in m indicates a lower sediment release rate under the same P. Therefore, 5

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

Subunit 1 0.8

Subunit 3

sigma

sigma

0.3

0.6

0

0.05

Recovery Rate

Recovery Rate

0

0.1 0.2

0.4

0.05 0.1

0.2

0.2

0.1

0.2

0.0

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

Reduction Rate

0.4

0.5

0.6

0.7

0.8

0.9

1

0.8

0.9

1

Reduction Rate

Subunit 2

Subunit 4

0.4

sigma

sigma

0.4

0.3

0

0.05

Recovery Rate

Recovery Rate

0

0.1 0.2

0.2

0.1

0.05

0.3

0.1 0.2

0.2

0.1

0.0

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

Reduction Rate

0.4

0.5

0.6

0.7

Reduction Rate

Fig. 3. Recovery rate for each subunit under different load reduction rate and intensity of external disturbance. σ values represent the disturbance intensities. The reduction rate was set according to the percentage of T1. The recovery rate was calculated as the probability of the final P lower than that at T2 for different reduction rates and σ values.

Subunit 1

Subunit 3 0.98

0.5 0.4

0.96

e

m

0.3

m

RI

RI

e

q 0.2

q

0.94

r

r

s

0.1

s 0.92

0.0 0.5

0.7

0.9

1.1

1.3

1.5

0.5

varying of parameter value

0.7

0.9

1.1

1.3

1.5

varying of parameter value

Subunit 2

Subunit 4

0.95 0.90 0.6 e

e

0.85 m

RI

RI

m q

0.80

q

0.4

r 0.75

r

s

s 0.2

0.70 0.5

0.7

0.9

1.1

1.3

1.5

0.5

varying of parameter value

0.7

0.9

1.1

1.3

1.5

varying of parameter value

Fig. 4. Variation in risk of irreversible degradation (RID) when parameter values varied under ecosystem recovery. Varying parameter values were set proportional to the original parameter sample drawn from their posterior. All parameter values ranged from 50% to 150% of its original values. RID were calculated under the variation of each parameter when other parameters remained unchanged. 6

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

Brock, 2006). Thus, the system is liable to secondary degradation. Ecosystem recovery could then be a way to prevent secondary degradation, because it could skew the dominant attractor from the undesirable regime to the desirable regime, making the desirable regime more resilient (Ratajczak et al., 2018; Carpenter et al., 2012). Nevertheless, reducing external environmental pressure is still important because it is the preliminary requirement for the existence of a desirable regime (Scheffer and Carpenter, 2003; Ratajczak et al., 2018). Additionally, without the reduction of external environmental pressure, the system will be held in a vulnerable position (Walker et al., 2012). Just like after rebuilding the bridge, we still need to carefully maintain the pressure within the threshold.

high sediment release could be the main reason for irreversible degradation, and decrease in sediment release could make the degraded regime less resilient. Parameter s, i.e., the loss rate through outflow and sedimentation, could lower RID with increasing parameter values (Fig. 4). This result suggested that a higher loss rate through outflow and sedimentation could make the system less resilient in a degraded regime. This had much more influence on subunits 2, 3, and 4, and was less influential for subunit 1. This was caused by the relatively high loss rate in subunit 1. Therefore, s cast stronger limits on the reversal of subunits 2, 3, and 4. The effects of parameter e, the exchange rate, on the four subunits were controversial. An increase of e could lead to lower RID in subunit 1, but higher RID in subunits 2 and 3. For subunit 4, e did not make an obvious difference (Fig. 4). This was due to the difference in the net contribution of exchange. For subunit 1, exchange between neighboring subunits acted as a net sink, whereas for subunits 2 and 3 it was a net source. Therefore, an increase in e could decrease RID in subunit 1 and increase RID in subunits 2 and 3. For all subunits, q displayed little influence on RID (Fig. 4). Parameter q was the shape parameter that represented the nonlinear response of sediment release to P. This result suggested that all subunits are highly nonlinear and changes in q in this range (−50% to +50%) were not sufficient to induce a significant influence on the reversibility of system degradation.

4.2. Opportunities and challenges for future generalization of this method Based on the successful application of our proposed method, we would like to further argue its generalization to other ecosystems. Although we demonstrated our method in a freshwater ecosystem, the method was not tailored to this specific ecosystem but more to a general framework that is applicable to a wide range of ecosystems. Increasing model studies on all types of ecosystems and the application of Bayesian inference for model parameterization ensured the basis for wider application of our method (Bagnara et al., 2019; Wu et al., 2017). Quantifying the risk of irreversible degradation could have a profound influence on our perspectives and perceptions of ecosystem degradation. Much emphasis has been placed on avoiding and reversing ecosystem degradation for sustainable development (Rossberg et al., 2017). To achieve this, the degradation threshold and its associated uncertainty has been widely studied as it determines the sustainable development boundary (Scheffer et al., 2001; Rockström et al., 2009). Safe operating space was proposed as a practical guide to this problem. This requires anthropogenic intervention to be maintained under the threshold for degradation, thus avoiding the existence of an alternative stable state (Rockström et al., 2009; Steffen et al., 2015). In this way, the risk of degradation was controlled at the highest level. For systems with high RID, the safe operating space might not be allowed because the chance that this system will have only one desirable regime without an alternative stable state is rare (Fig. 1). This raised further challenges for safe operating space and imposed a higher requirement for sustainable development. The spatial aspect of our method is especially important. Spatial patterns associated with responses and risks are an important factor for conservation. Identification of such spatial patterns could provide information on conservation priority and mitigation opportunities for leveraging effects (Gaynor et al., 2019; Wearn et al., 2012). We testified the spatial heterogeneity of RID in a freshwater ecosystem. This spatial heterogeneity could be even more evident in other ecosystems because aquatic ecosystems are more free-flowing than other rigid ecosystems that have a lower rate of internal exchanges. Therefore, spatial heterogeneity due to different morphology and other features intrinsic to each subunit are less averaged by these internal exchange processes. This provides a valid reason for the wider application of our method on other types of ecosystems. Promising as our method is, the main challenge remains with its uncertainty. In a sense, RID quantified a certain fraction of the uncertainty associated with T2 , the probability of it falling below 0. Because RID itself is derived from uncertainty, it becomes complicated if we wish to further analyze its uncertainty. Nevertheless, with such significant results, both for RID at the whole lake level and spatial differentiated RID (Fig. 2), it should be convincing enough to draw our attention to this risk of irreversible degradation.

4. Discussion 4.1. Mitigation for systems with high undesirable resilience Our results showed that ecosystem recovery could be a more efficient way for mitigation of systems with high RID. As reflected in our model results, varying model parameter values significantly changed the undesirable resilience (Fig. 4). Ecosystem recovery led to improvements in ecosystem attributes and, therefore, T2 gradually tended toward the positive side. We could have better intuition on this when relating to the fundamental idea of catastrophe bifurcation, which is the theoretical basis for regime shift (May, 1977; Scheffer et al., 2001). Catastrophe bifurcation was originally widely observed in the field of physics and engineering. It describes the situation where the system first changes smoothly with external drivers until it reaches a threshold where it jumps discontinuously to another state. The recovery does not follow the same trajectories; however, the driver needs to be reverted far back from the previous threshold (Strogatz, 2018). This phenomenon is known as hysteresis and has been more recently observed in ecological and social systems (Ratajczak et al., 2018; Scheffer et al., 2001; van Nes et al., 2016). In physics and engineering, it is quite common that the critical transition in a catastrophe bifurcation will be irreversible. For example, the collapse of a bridge after the pressure exceeds the threshold. (Strogatz, 2018). There is no other way for the recovery of the bridge other than by building a new bridge. For ecosystems, the rebuilding process might be obscured by their complexity and on-going dynamics. The essence of a regime shift lies in the change of ecosystem structure and function (Folke et al., 2004; Oliver et al., 2015). Therefore, after the collapse of the ecosystem, rebuilding the system is the most straightforward method. In particular, when confronted with high risks of irreversible degradation because reducing environmental pressures might never achieve this goal. In addition, ecosystem recovery is a good way to build resilience into the system so that it will be less vulnerable to change. Although our results only revealed the positive effect of disturbances, because our system was subject to high RID, the effect of disturbances can act both ways. It could bring the system from a desirable regime to a degraded regime, and vice versa. For systems with high RID, the degraded state could be a much stronger attractor (Ibelings et al., 2007; Carpenter and

5. Conclusions We proposed a novel indicator RID that measured the risk of an ecosystem degradation being irreversible, in terms of the uncertainty 7

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

associated with T2 . By adopting Bayesian inference, the whole framework for implementing this indicator was straightforward and flexible, making it convenient for management practice. By applying RID to the example of lake eutrophication of Lake Erhai, we testified that ecosystem could suffer from high RID, as shown by the RID of 72% at the whole lake level in our case study. This finding suggested that there should be more awareness regarding ecosystem degradation because the consequences could be more severe, and therefore raise further challenges for sustainability. On the other hand, by applying RID, we also identified the following opportunities for mitigation.

transition. Ecol. Lett. 9, 311–318. Carpenter, S.R., Ludwig, D., Brock, W.A., 1999. Management of eutrophication for lakes subject to potentially irreversible change. Ecol. Appl. 9, 751–771. Carpenter, S.R., Walker, B., Anderies, J.M., Abel, N., 2001. From metaphor to measurement: resilience of what to what? Ecosystems 4, 765–781. Carpenter, S.R., Arrow, K.J., Barrett, S., Biggs, R., Brock, W.A., Crépin, A.-S., Engström, G., Folke, C., Hughes, T.P., Kautsky, N., et al., 2012. General resilience to cope with extreme events. Sustainability 4, 3248–3259. Chapin, F.S., Carpenter, S.R., Kofinas, G.P., Folke, C., Abel, N., Clark, W.C., Olsson, P., Smith, D.M.S., Walker, B., Young, O.R., et al., 2010. Ecosystem stewardship: sustainability strategies for a rapidly changing planet. Trends Ecol. Evol. 25, 241–249. Cohen-Shacham, E., Walters, G., Janzen, C., Maginnis, S., 2016. Nature-based Solutions to Address Global Societal Challenges. IUCN, Gland, Switzerland, pp. 97. Cumming, G.S., Peterson, G.D., 2017. Unifying research on social-ecological resilience and collapse. Trends Ecol. Evol. 32, 695–713. Cumming, G.S., von Cramon-Taubadel, S., 2018. Linking economic growth pathways and environmental sustainability by understanding development as alternate social-ecological regimes. Proc. Natl. Acad. Sci. 115, 9533–9538. Dorazio, R.M., Johnson, F.A., 2003. Bayesian inference and decision theory-a framework for decision making in natural resource management. Ecol. Appl. 13, 556–563. Ellison, A.M., 2004. Bayesian inference in ecology. Ecol. Lett. 7, 509–520. Elmqvist, T., Folke, C., Nyström, M., Peterson, G., Bengtsson, J., Walker, B., Norberg, J., 2003. Response diversity, ecosystem change, and resilience. Front. Ecol. Environ. 1, 488–494. Filatova, T., Polhill, J.G., van Ewijk, S., 2016. Regime shifts in coupled socio-environmental systems: review of modelling challenges and approaches. Environ. Model. Softw. 75, 333–347. Folke, C., Carpenter, S., Walker, B., Scheffer, M., Elmqvist, T., Gunderson, L., Holling, C.S., 2004. Regime shifts, resilience, and biodiversity in ecosystem management. Annu. Rev. Ecol. Evol. Syst. 35, 557–581. Folke, C., Carpenter, S.R., Walker, B., Scheffer, M., Chapin, T., Rockström, J., 2010. Resilience thinking: integrating resilience, adaptability and transformability. Ecol. Soc. 15. Gaynor, K.M., Brown, J.S., Middleton, A.D., Power, M.E., Brashares, J.S., 2019. Landscapes of fear: Spatial patterns of risk perception and response. Trends Ecol. Evol. Gelman, A., Hwang, J., Vehtari, A., 2014. Understanding predictive information criteria for bayesian models. Stat. Comput. 24, 997–1016. Genkai-Kato, M., Carpenter, S.R., 2005. Eutrophication due to phosphorus recycling in relation to lake morphometry, temperature, and macrophytes. Ecology 86, 210–219. Hartig, F., Minunno, F., Paul, S., 2017. BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics. URL:https://CRAN.R-project.org/ package=BayesianTools r package version 0.1.4. Hazbavi, Z., Baartman, J.E., Nunes, J.P., Keesstra, S.D., Sadeghi, S.H., 2018. Changeability of reliability, resilience and vulnerability indicators with respect to drought patterns. Ecol. Ind. 87, 196–208. Hijmans, R.J., Phillips, S., Leathwick, J., Elith, J., 2017. dismo: Species Distribution Modeling. URL:https://CRAN.R-project.org/package=dismo r package version 1. 1-4. Holling, C.S., 1973. Resilience and stability of ecological systems. Annu. Rev. Ecol. Syst. 4, 1–23. Hughes, T.P., Carpenter, S., Rockström, J., Scheffer, M., Walker, B., 2013. Multiscale regime shifts and planetary boundaries. Trends Ecol. Evol. 28, 389–395. Ibelings, B.W., Portielje, R., Lammens, E.H., Noordhuis, R., van den Berg, M.S., Joosse, W., Meijer, M.L., 2007. Resilience of alternative stable states during the recovery of shallow lakes from eutrophication: lake veluwe as a case study. Ecosystems 10, 4–16. Janssen, A.B., De Jager, V.C., Janse, J.H., Kong, X., Liu, S., Ye, Q., Mooij, W.M., 2017. Spatial identification of critical nutrient loads of large shallow lakes: implications for lake Taihu (China). Water Res. 119, 276–287. Keesstra, S., Nunes, J., Novara, A., Finger, D., Avelar, D., Kalantari, Z., Cerdà, A., 2018. The superior effect of nature based solutions in land management for enhancing ecosystem services. Sci. Total Environ. 610, 997–1009. Kéfi, S., Rietkerk, M., Roy, M., Franc, A., De Ruiter, P.C., Pascual, M., 2011. Robust scaling in ecosystems and the meltdown of patch size distributions before extinction. Ecol. Lett. 14, 29–35. Liang, Z., Wu, S., Chen, H., Yu, Y., Liu, Y., 2018. A probabilistic method to enhance understanding of nutrient limitation dynamics of phytoplankton. Ecol. Model. 368, 404–410. May, R.M., 1977. Thresholds and breakpoints in ecosystems with a multiplicity of stable states. Nature 269, 471. McCarthy, M.A., Possingham, H.P., 2007. Active adaptive management for conservation. Conserv. Biol. 21, 956–963. Nyström, M., Norström, A.V., Blenckner, T., de la Torre-Castro, M., Eklöf, J.S., Folke, C., Österblom, H., Steneck, R.S., Thyresson, M., Troell, M., 2012. Confronting feedbacks of degraded marine ecosystems. Ecosystems 15, 695–710. Oliver, T.H., Heard, M.S., Isaac, N.J., Roy, D.B., Procter, D., Eigenbrod, F., Freckleton, R., Hector, A., Orme, C.D.L., Petchey, O.L., et al., 2015. Biodiversity and resilience of ecosystem functions. Trends Ecol. Evol. 30, 673–684. Oliver, S.K., Collins, S.M., Soranno, P.A., Wagner, T., Stanley, E.H., Jones, J.R., Stow, C.A., Lottig, N.R., 2017. Unexpected stasis in a changing world: Lake nutrient and chlorophyll trends since 1990. Glob. Change Biol. 23, 5455–5467. Qian, S.S., Stow, C.A., Borsuk, M.E., 2003. On monte carlo methods for bayesian inference. Ecol. Model. 159, 269–277. Ratajczak, Z., D’odorico, P., Collins, S.L., Bestelmeyer, B.T., Isbell, F.I., Nippert, J.B., 2017. The interactive effects of press/pulse intensity and duration on regime shifts at multiple scales. Ecol. Monogr. 87, 198–218.

i) We found spatial heterogeneity could provide opportunity windows as certain regions could have higher chances of recovery, and we also need to pay more attention to regions that are more susceptible to irreversible degradation. ii) External pulse disturbances could have an overall positive influence on ecosystem recovery for ecosystems with high RID, because it could drive the ecosystem beyond the boundary for alternative stable states and, therefore, trigger the shift towards desirable regime without further reducing environmental pressures. iii) Our results also suggested that ecosystem recovery resulting from ecosystem-based management could be more efficient for such systems. Changes in certain ecosystem attributes could lower RID significantly. Based on the findings of our study, we concluded that RID could be a meaningful and useful indicator. It also has high potential for future application to a wide range of ecosystems. Acknowledgments We acknowledge the help of the editor and anonymous reviewers. Their comments improved the quality of this paper significantly. We appreciate the insightful comments given by Prof. Zhang Shiqiu during the thesis defense. Her comments partially motivated the risk perspective we employed in this paper. We also acknowledge our colleagues Dr. Zhang Xiaoling and Dr. Zou Rui for providing the monitoring data of Lake Erhai, Dr. Wu Zhen and Dr. Sheng Hu for their help in external load estimation, and Dr. Chen Huili for helpful discussion on spatial division. The work presented in this paper was supported by National Basic Research Program (973) of China (No. 2015CB458900) and the National Natural Science Foundation of China (No. 51779002). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, athttps://doi.org/10.1016/j.ecolind.2019.105621. References Allen, T.F., Hoekstra, T.W., 2015. Toward a Unified Ecology. Columbia University Press. Andersen, T., Carstensen, J., Hernandez-Garcia, E., Duarte, C.M., 2009. Ecological thresholds and regime shifts: approaches to identification. Trends Ecol. Evol. 24, 49–57. Angeler, D.G., Allen, C.R., 2016. Quantifying resilience. J. Appl. Ecol. 53, 617–624. Arhonditsis, G., Stow, C., Steinberg, L., Kenney, M., Lathrop, R., McBride, S., Reckhow, K., 2006. Exploring ecological patterns with structural equation modeling and bayesian analysis. Ecol. Model. 192, 385–409. Bagnara, M., Gonzalez, R.S., Reifenberg, S., Steinkamp, J., Hickler, T., Werner, C., Dormann, C.F., Hartig, F., 2019. An r package facilitating sensitivity analysis, calibration and forward simulations with the lpj-guess dynamic vegetation model. Environ. Model. Softw. 111, 55–60. Beisner, B.E., Peres-Neto, P.R., Lindström, E.S., Barnett, A., Longhi, M.L., 2006. The role of environmental and spatial processes in structuring lake communities from bacteria to fish. Ecology 87, 2985–2991. Bishop, C.M., 2006. Pattern Recognition and Machine Learning. Springer. Carpenter, S.R., 2005. Eutrophication of aquatic ecosystems: bistability and soil phosphorus. Proc. Natl. Acad. Sci. 102, 10002–10005. Carpenter, S.R., Brock, W.A., 2006. Rising variance: a leading indicator of ecological

8

Ecological Indicators 107 (2019) 105621

S. Wu, et al.

analysis of ordinary differential equations. R package 1.6. Steffen, W., Richardson, K., Rockström, J., Cornell, S.E., Fetzer, I., Bennett, E.M., Biggs, R., Carpenter, S.R., De Vries, W., De Wit, C.A., et al., 2015. Planetary boundaries: guiding human development on a changing planet. Science 347, 1259855. Strogatz, S.H., 2018. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. CRC Press. Tilman, D., May, R.M., Lehman, C.L., Nowak, M.A., 1994. Habitat destruction and the extinction debt. Nature 371, 65. van Nes, E.H., Arani, B.M., Staal, A., van der Bolt, B., Flores, B.M., Bathiany, S., Scheffer, M., 2016. What do you mean’, tipping point’? Trends Ecol. Evol. 31, 902–904. Walker, B., Carpenter, S., Anderies, J., Abel, N., Cumming, G., Janssen, M., Lebel, L., Norberg, J., Peterson, G.D., Pritchard, R., 2002. Resilience management in socialecological systems: a working hypothesis for a participatory approach. Conserv. Ecol. 6. Walker, B.H., Carpenter, S.R., Rockstrom, J., Crépin, A.-S., Peterson, G.D., 2012. Drivers, slow variables, fast variables, shocks, and resilience. Ecol. Soc. 17. Wearn, O.R., Reuman, D.C., Ewers, R.M., 2012. Extinction debt and windows of conservation opportunity in the brazilian amazon. Science 337, 228–232. Wu, Z., Liu, Y., Liang, Z., Wu, S., Guo, H., 2017. Internal cycling, not external loading, decides the nutrient limitation in eutrophic lake: a dynamic model with temporal bayesian hierarchical inference. Water Res. 116, 231–240. Yeung, A.C., Richardson, J.S., et al., 2016. Some conceptual and operational considerations when measuring ’resilience’: a response to hodgson. Trends Ecol. Evol. 31, 2–3. Zamparas, M., Zacharias, I., 2014. Restoration of eutrophic freshwater by managing internal nutrient loads. A review. Sci. Total Environ. 496, 551–562.

Ratajczak, Z., Carpenter, S.R., Ives, A.R., Kucharik, C.J., Ramiadantsoa, T., Stegner, M.A., Williams, J.W., Zhang, J., Turner, M.G., 2018. Abrupt change in ecological systems: inference and diagnosis. Trends Ecol. Evol. R Core Team, 2017. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing Vienna, Austria. URL:https://www.R-project. org/. Rockström, J., Steffen, W., Noone, K., Persson, Å., Chapin III, F.S., Lambin, E.F., Lenton, T.M., Scheffer, M., Folke, C., Schellnhuber, H.J., et al., 2009. A safe operating space for humanity. Nature 461, 472. Rossberg, A.G., Uusitalo, L., Berg, T., Zaiko, A., Chenuil, A., Uyarra, M.C., Borja, A., Lynam, C.P., 2017. Quantitative criteria for choosing targets and indicators for sustainable use of ecosystems. Ecol. Ind. 72, 215–224. Rougé, C., Mathias, J.-D., Deffuant, G., 2013. Extending the viability theory framework of resilience to uncertain dynamics, and application to lake eutrophication. Ecol. Ind. 29, 420–433. Sasaki, T., Furukawa, T., Iwasaki, Y., Seto, M., Mori, A.S., 2015. Perspectives for ecosystem management based on ecosystem resilience and ecological thresholds against multiple and stochastic disturbances. Ecol. Ind. 57, 395–408. Scheffer, M., Carpenter, S.R., 2003. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol. 18, 648–656. Scheffer, M., Carpenter, S., Foley, J.A., Folke, C., Walker, B., 2001. Catastrophic shifts in ecosystems. Nature 413, 591. Scheffer, M., Carpenter, S.R., Dakos, V., van Nes, E.H., 2015. Generic indicators of ecological resilience: inferring the chance of a critical transition. Annu. Rev. Ecol. Evol. Syst. 46, 145–167. Soetaert, K., 2009. rootSolve: nonlinear root finding, equilibrium and steady-state

9