Kelp Forests: Catastrophes, Resilience, and Management

Resilient kelp forests provide foundation habitat for marine ecosystems and are indicators of the ecosystems’ sustainable natural capital. Loss of resilience and imperfectly reversible catastrophic shifts from kelp forests to urchin barrens, due to pollution or loss of a top predator, are part of an ecological tipping point phenomenon, and involve a loss in sustainable natural capital. Management controls to prevent or reverse these shifts and losses are classified in a number of ways. Systemic controls eliminate the cause of the problem. Symptomatic controls use leverage points for more direct control of the populations affected, urchin harvesting or culling, or kelp enhancement. There is a distinction between ongoing structural (press) controls versus temporary or intermittent perturbation (pulse) controls, and one between shift preventing versus shift reversing or restorative controls. Adaptive management and the options it creates both focus on reductions in uncertainty and control policies with the flexibility to take advantage of those reductions. The various management distinctions are most easily understood by modeling the predator-urchin-kelp marine ecosystem. This paper develops a mathematical model of the ecosystem that has the potential for two different catastrophic shifts between equilibria. Pulse disturbances, originating from exogenous abiotic factors or population dynamics elsewhere in the metacommunity, can activate shifts. A measure of probabilistic resilience is developed and used as part of an assessment of the ecosystem’s sustainable stock of natural capital. With perturbation outcomes clustered around the originating equilibrium, hysteresis is activated, resulting imperfect reversibility of catastrophic shifts, and a loss in natural capital. The difficulty of reversing a shift from kelp forest to urchin barren, with an associated loss in sustainable natural capital, is an example. Management controls are modeled. I find that systemic and symptomatic, and press and pulse, controls can be complementary. Restorative controls tend to be more difficult or costly than preventative ones. Adaptive management, favoring flexible, often preventative, controls, creates option value, lowering control costs and/or losses in sustainable natural capital. Two cases are used to illustrate, Tasmania, Australia and Haida Gwaii, Canada.


INTRODUCTION
Healthy marine ecosystems are rich in natural capital (Mäler and Li, 2010;Bennett et al., 2016;Wernberg et al., 2019). In an ecologically resilient kelp forest community, kelp (order Laminariales) are the foundation species whose physical characteristics provide the habitat for a healthy marine ecosystem (Miller et al., 2018). The density of kelp is an indicator of the richness of the community's stock of natural capital, and the kelp forest's ecological resilience (Holling, 1973) is an indicator of the stock's sustainability (Graham, 2004;Ellison et al., 2005;Small, 2018). A sustained stock has value through the ongoing human life support services provided (e.g., harvestable fish and other food products) (Boyd and Krupnick, 2013;Fenichel and Abbott, 2014;Bond, 2017). However, kelp are consumed by urchins (order Echinoida or Diadematoida), and the kelp forest regime, with cryptic (i.e., low density) urchin herbivores whose density is limited by the availability of safe hiding spaces, typically exists in an ecosystem in which its resilience (Walker et al., 2004) can be compromised. Decreased resilience increases the probability of a catastrophic (abrupt, perhaps over decades) shift from a cryptic urchin-kelp forest regime to a kelp depauperate-urchin barren regime with a much poorer stock of natural capital (Petraitis and Dudgeon, 2004;Ling et al., 2015;Krumhansl et al., 2016;Wernberg et al., 2016;Filbee-Dexter and Wernberg, 2018;Wernberg et al., 2019).
To address the loss of resilience and depletion of kelp forests, this paper develops a mathematical model that integrates what is known about the predator-urchin-kelp ecosystem to analyze catastrophic shifts in urchin and kelp densities. While empirical evidence to support such a model is not definitive, low urchin kelp-forests and kelp depauperate-urchin barrens are prime candidates for alternative regimes (Petraitis and Dudgeon, 2004). The model is a recursively decoupled catastrophe model with two trophic interactions and the potential for two catastrophic shifts, one at the level of the predatorurchin trophic interaction and one at the level of the urchinkelp interaction. It incorporates the concept of resilience at both levels and uses resilience probabilities and hysteresis as part of a definition of the sustainable natural capital stock as expected kelp density. The key concepts used in the model are described in Supplementary Appendix A. Two cases are used to illustrate, Tasmania, Australia and Haida Gwaii, Canada (Lee et al., 2018;Tasmanian Government, 2019). The first is an ocean warming induced invasion of long spined sea urchins (Centrostephanus rodgersii) to coastal Tasmania and subsequent overgrazing of kelp forests (Edgar et al., 2004). The second is the absence of the sea otter predator (Enhydra lutris) in Haida Gwaii and lack of kelp forest recovery. Supplementary Appendix B provides details for both cases. The model is used to assess a range of management controls. See Supplementary Appendix A for details about the controls. I find that there is often complementarity in that combinations of control types can be more successful than one type alone. Moreover, in the face of hysteresis and uncertain restoration costs, preventative controls tend to have an advantage over restorative ones in that they require less investment and are often more amenable to adaptive management with the flexible decision process that creates option value.
The section "Materials and Methods" presents the formal definitions of catastrophe theory, trophic cascades, pulse disturbances, resilience, and hysteresis as they apply to the model of the predator-urchin-kelp ecosystem. It also develops the top-down recursively decoupled catastrophe model (Brummitt et al., 2015). The model is intended to be demonstrative or illustrative. The parameter values used for the catastrophe model and the pulse disturbances are not meant to be precisely correct, but to capture the essence of the various interactions so as to better understand their outcomes and the influence of alternative management policies (Murray, 2002). Finally it characterizes kelp forests and their resilience as sustainable natural capital (Mäler and Li, 2010). The section "Results" uses the catastrophe model, probabilistic resilience, and natural capital to characterize the losses in sustainable natural capital that may occur, or have occurred, in Tasmania and Haida Gwaii. It also examines the use of systemic versus symptomatic controls and press versus pulse controls, and assesses adaptive management control and the option value of preventative versus restorative controls. The section "Discussion" summarizes, considers limitations and extensions, and concludes.

Trophic Cascades
The predator-urchin-kelp model has been characterized as one of best examples of a strong top-down trophic cascade (Estes and Palmisano, 1974;Sala et al., 1998;Hereu Fina, 2004;Borer et al., 2005;Halpern et al., 2006;Sommer, 2008). A top-down trophic cascade describes a set of relationships for which there is a negative relationship between the biomass density of adjacent trophic levels and a positive one between those that are once removed (Carpenter et al., 1985;Sommer, 2008). Predation or grazing mortality, or release from it, is the dominant response (Holling, 1959a,b;Nowak et al., 2008). Due to factors such as larger body size and slower change at higher trophic levels, a trophic cascade exhibits weak numerical responses, and bottom-up effects, focusing on numerical responses (the response of predator or grazer), play a minimal role (Sommer, 2008).
While factors such as loss of top predators, storms or pollution may interfere (Foster and Schiel, 2010;Reed et al., 2011;Dunn and Hovel, 2019), there is statistical analysis supporting the dominance of top-down influence (Halpern et al., 2006), and there is evidence for the absence of numerical response in the behavior of urchins and predators. At the urchin-kelp level, urchin grazers use alternative food resources. With no kelp, urchins in high density barrens consume biofilms and crustose coralline algae (Lyons and Scheibling, 2007;Kraufvelin, 2017). Urchin density in barrens doesn't decline with absence of kelp, instead urchins adapt to less nutritious food sources (Konar and Estes, 2003;Stewart and Konar, 2012). Evidence also suggests that while urchin reproductive potential is higher in kelp forests, the pattern of larval dispersal and settlement is such that settlement in barrens is possible, and is enhanced by chemical cues from coralline algae found in barrens (Rowley, 1989(Rowley, , 1990Baskett and Salomon, 2010;Wilson et al., 2013). In addition, dense kelp forests provide habitat for urchin predators and are able to physically beat back grazing (Andrew, 1993;Konar and Estes, 2003). At the predator-urchin level, both sea otters and lobster are generalist predators, making them less sensitive to changes in urchin density (Edyvane, 2003;Estes et al., 2003;Tinker et al., 2008;Marzloff, 2012;Butler et al., 2016). At high densities, urchins themselves becomes less nutritious and more difficult prey (e.g., longer spines and flight), avoiding encouraging greater predator density (Konar and Estes, 2003;Watson and Estes, 2011). Ling et al. (2019) find the highest lobster predation to be on small urchins within kelp beds.

Urchin-Kelp Interactions
First consider only the lower-level trophic interaction, with urchin as grazers of kelp, and equilibrium urchin density from the predator-urchin subsystem as an exogenous input to the urchin-kelp system. At the urchin-kelp trophic interaction level, kelp forests (high kelp density) and kelp depauperate communities (low kelp density) are examples of alternate states. A forward shift from a kelp forest equilibrium to a kelp depauperate equilibrium can occur, as well as a reverse shift from depauperate equilibrium to kelp forest. Starting from a high kelp density equilibrium, gradual increases in an ecosystem's exogenous parameter (urchin density), which may be result of species invasion encouraged by environmental changes such as ocean warming, or human over-harvesting of urchin predators, causes reductions in kelp density, and some combination of abrupt and gradual decreases in the resilience of the kelp forest equilibrium Filbee-Dexter and Wernberg, 2018). Reduced resilience for the kelp forest equilibrium lessens the propensity of the state variable (kelp density) to return to the original equilibrium after a stochastic perturbation such as a kelp destroying storm event, and increases the probability of a forward tip to a kelp depauperate one (Filbee-Dexter and Scheibling, 2014). Once the kelp forest equilibrium has lost all its resilience, the forward tip becomes unavoidable.
To illustrate, Figure 1 shows, in blue, an isocline of urchinkelp steady state equilibria (the model producing this isocline will be presented in section "Urchin-Kelp Model"). With U as exogenously determined urchin density (in 100 g/m 2 ), and K as the state variable, kelp density (in g/m 2 ), the upper arm, K ≥ 5, contains stable kelp forest equilibria, and the lower arm, K ≤ 0.1, contains stable kelp depauperate equilibria. The upward sloping portion of the curve, from (U, K) = (0.6, 0.1) to (7.7, 5), contains unstable kelp equilibria. For a given U, the unstable kelp equilibrium is a boundary point between the basins of attraction of the stable kelp equilibrium on the upper arm versus one on the lower arm. The inevitable tip points are indicated by vertical solid black arrows.

Kelp Resilience and Catastrophic Shifts
At a given urchin density, a kelp forest equilibrium's resilience is the tendency for a pulse perturbation in kelp density from that equilibrium to be attracted back to it, rather than being induced toward a kelp depauperate equilibrium. See Arrows indicate movement away from equilibria on the unstable arm and toward equilibria on the stable arms. The (U, K) coordinates for the letters are: a− (0.4, 9.9), b− (0.7, 9.7), c− (0.7, 0.17), d− (0.7, 0.06), e− (4, 0.008), f− (7.4, 0.004), and g− (7.7, 0.004).
Supplementary Appendix A for more on catastrophic shifts and resilience. Resilience can be translated into a resilience probability that depends on the size of the equilibrium's basin of attraction and on the probability distribution for the size of the perturbation or, equivalently, the distribution of perturbation displacement locations . If U ≤ 0.6, there is a fully resilient kelp forest equilibrium, and any kelp density perturbation has a probability of one of returning to it. If U ≥ 7.7, the kelp forest equilibrium loses all resilience and disappears. When the kelp forest equilibrium vanishes, the kelp depauperate equilibrium achieves full resilience, and any kelp density perturbation has a probability of one of returning to it. The zone of partial resilience extends over the range 0.6 < U < 7.7 Over that range, there is a change from a fully resilient kelp forest to a fully resilient kelp depauperate regime. At any given urchin density in the partial resilience zone, there are two stable equilibria and two basins of attraction, separated by a boundary point consisting of the unstable equilibrium. Within the zone of partial resilience, the resilience of a stable equilibrium depends on the magnitude of perturbations away from it relative to the size of its basin of attraction . Given a probability density function for perturbation induced displacement locations, the resilience of a kelp forest/depauperate equilibrium can be translated into probabilistic resilience, or the probability that, given a perturbation from a stable kelp forest/depauperate equilibrium, the displacement location will be in the basin of attraction of the stable kelp forest/depauperate equilibrium, and attracted back to it. Otherwise, the perturbation will land in the basin of attraction of the stable depauperate/kelp forest equilibrium, and there will be a shift to the alternative equilibrium.
Designating a stable upper arm equilibrium as a function of U, K h/U = K h (U), a stable lower arm equilibrium as K l/U = K l (U), and an unstable equilibrium boundary point as K b/U = K b (U), the basin of attraction of the uppper arm equilbrium is the vertical distance K h/U − K b/U Standish et al., 2014;Baho et al., 2017). For the lower arm equilibrium, it is K b/U − K l/U . With urchin density fixed at U = 4 in Figure 1, consider a pulse disturbance in kelp density from either of the two stable equilibria that lands at a point somewhere along the vertical line between K l/4 = 0.008 (point e) and K h/4 = 8.46, including the endpoints. The basins of attraction are K h/4 − K b/4 = 6.93 and K b/4 − K l/4 = 1.522. From each stable equilibrium in the partial resilience zone, a perturbation's displacement location may or may not be within that equilibrium's basin of attraction. Specifically, I model the displacement locations of a perturbation from the kelp forest equilibrium in the partial resilience zone as following a combination of two piecewise linear probability density functions. The first, representing the clustering of very small or zero disturbances, gives a displacement location K h/U and occurs with a probability of π h . The second, covering the full range of disturbance locations K h/U to K l/U , is the uniform distribution with a probability density of u h K h/U , K l/U (Schultz et al., 2017). The latter occurs with a probability 1 − π h . Zero disturbances will always end up at K h/U . Each non-zero perturbation outcome will eventually end up at either K h/U or K l/U , depending on the basin of attraction into which it is perturbed. For a given urchin density, negative perturbations from an equilibrium on the high kelp arm toward an equilibrium on the low kelp arm, will have a probability, P K h/h/U , of returning to the upper arm. This probability depends on the probability of a zero perturbation, and the probability, from the uniform distribution, of a perturbation being small enough that its displacement location remains within the basin of attraction of the upper arm equilibrium. The probability of a perturbation from an equilibrium on the high kelp arm ending up on the lower arm is Taking into account all three resilience zones (full, partial, and no resilience), and assuming partial resilience zone perturbations are no greater than K h/U − K l/U , the probability a perturbation from K h/U ending up at K h/U versus K l/U is given as Eq. 1.
Starting from a kelp depauperate equilibrium on the lower arm, K l/U , and assuming the probability, π l , of zero perturbations from K l/U , the probability of ending up at K h/U versus K l/U over three resilience zones is given by Eq. 2.
Over the partial resilience range, 0.6 < U < 7.7, the probability of returning to the original equilibrium will vary between zero and one. The abruptness of the changes in probability in response to changes in U is directly related to the degree of clustering in the perturbation distributions, as is the degree of hysteresis. Hysteresis, or hysteretic memory, is closely related to resilience. If the system has hysteric memory, the probability of a perturbed point landing in one or the other basin of attraction will depend not just on the relative sizes of the basins, but also on which equilibrium it was perturbed from. Supplementary Appendix A contains a diagrammatic example of hysteresis. With displacement locations bounded by K l/4 = 0.008 and K h/4 = 8.46, hysteretic memory is provoked by the perturbation's probability distribution of displacement locations clustering around the perturbation generating equilibrium (Paine et al., 1998). As a result, a perturbation from a kelp forest equilibrium, which lands in the basin of attraction of, and is absorbed by, the kelp depauperate equilibrium, will have an increased probability of returning to the depauperate equilibrium following future perturbations from it (Walker et al., 2004).
The magnitudes of π h and π l in Eqs. 1, 2 determine the degree of clustering around the kelp forest and kelp depauperate equilibrium, respectively. Hysteresis requires that P K h/h/U > P K h/l/U , or equivalently that P K l/l/U > P K l/h/U . Disturbance locations in a kelp forest/kelp depauperate equilibrium's basin of attraction more probably were perturbed from that equilibrium than from the kelp depauperate/kelp forest equilibrium. From Eqs 1, 2, π h + π l > 0 will produce some hysteresis (see the Supplementary Appendix C for proof). I define a moderate amount of clustering and hysteresis as 0 < π h + π l < 2. I use a moderate degree of clustering and hysteresis with π h = π l = 0.5. Exteme clustering, π l = π h = 1, yields strong hystersis, making P K h/h/U = P K l/l/U = 1 and P K h/l/U = P K l/h/U = 0. No clustering, π l = π h = 0, yields no hysteresis, P K h/h/U = P K h/l/U and P K l/l/U = P K h/l/U . The extreme and no clustering cases are presented in more detail in Supplementary Appendix C.

Predator-Urchin Interactions
The urchin-kelp system is the lower part of a dynamic trophic cascade. The higher-level trophic interaction in the cascade is the predator-urchin linkage. While predator size relative to prey size has also been shown to be important Dunn et al., 2017), it need not always be, so I simply use predator density as the exogenous driver, and urchin density as the state variable. Depletion of the predator through human harvesting creates greater densities of its prey (urchin) and, reverberating down the food chain, the overgrazing of the lower level food source (kelp) (Estes and Palmisano, 1974;Estes and Duggins, 1995;Steneck et al., 2002;Lafferty, 2004;Estes et al., 2011;Ling et al., 2015;Ripple et al., 2016;Filbee-Dexter and Wernberg, 2018). The overharvesting of sea-otter drove the switch from cryptic urchin-kelp forest to kelp depauperate-urchin barren off Haida Gwaii (Estes and Palmisano, 1974;Sloan and Dick, 2012), and the harvesting of spiny lobster off the Tasmanian coast has reduced the capacity of predation to control urchin density and to maintain kelp bed resilience Ling and Johnson, 2012). Resilience and hysteresis for a cryptic urchin versus urchin barren equilibria can characterized as described above for the kelp forest versus kelp depauperate equilibria. Catastrophic shifts are possible at the predator-urchin level as well as at the urchin-kelp level (Soulé et al., 2003;Brummitt et al., 2015;Ling et al., 2015), and it is the combined catastrophic shifts that matter in managing the overall system to prevent kelp loss.

Modeling Assumptions
Using what is known about the predator-urchin-kelp ecosystem off the Tasmanian Coast (Edgar et al., 2004;Johnson et al., 2011;Marzloff, 2012;Tasmanian Government, 2019) and that in Haida Gwaii (Lee et al., 2016;Lee et al., 2018), I develop a recursively decoupled mathematical model that allows for both predator-urchin and urchin-kelp interactions with the possibility for hysteresis and catastrophic shifts at both trophic interaction levels (Ludwig et al., 1997;Scheffer et al., 2001;Soulé et al., 2003;Pearse, 2006;Flukes et al., 2012;Ling et al., 2015;Johnson et al., 2017;Kenner and Tinker, 2018). Both environmental factors such as ocean warming (Edgar et al., 2004;Gorman et al., 2009;Filbee-Dexter et al., 2016), and predator density, are top level parameters that can drive the system to shift from a kelp forest with cryptic urchin to a kelp depauperate urchin barren.
I follow Scheffer et al. (2001) and specify the minimal model consistent with gradual increases in urchin density triggering a catastrophic decline in kelp density  and with gradual changes in an environmental factor, or decreases in predator density, triggering an explosion in urchin density (Selkoe et al., 2015). Since there are at most very weak numerical responses of urchin density to kelp density and of predator density to urchin density, my model treats them as zero, allowing me to use a downwardly recursively decoupled model; a top down cascade approach with no reciprocal effect (Ludwig et al., 1978;Walker et al., 2012;Brummitt et al., 2015;Pershing et al., 2015;Winnie and Creel, 2017). Changes in the density of urchin prey are driven by environmental factors and predator (lobster or sea otter) density, assumed to be the slower, or exogenous, variables, affecting urchin density but unaffected by it (Estes and Duggins, 1995;Dunn et al., 2017). In the same vein, urchins, as grazers of kelp, are drivers of changes in kelp cover, and urchin density can be treated as a slow, or exogenous, variable, affecting kelp density but unaffected by it, allowing changes at higher trophic levels cascade down the food chain (Paine, 1980;McLaren and Peterson, 1994;Ripple et al., 2016).
For the lower-level trophic interaction, the urchin-kelp model, I assume that kelp exhibit logistic growth in the absence of grazing. The potential for hysteresis and catastrophic shifts comes from a non-linear relationship between urchins and kelp. Specifically, urchin grazing on kelp is modeled as a Holling type-III functional response for which the proportion of kelp consumed per urchin (grazing efficiency or proportion of available kelp consumed per urchin or the average product of kelp density) depends on kelp density, starting off very low at very low kelp densities, increasing at a decreasing rate at low kelp densities, reaching a maximum, and decreasing, at first at an increasing and then at a decreasing absolute rate, approaching zero at very high kelp densities (Holling, 1959a,b;Ludwig et al., 1997;Scheffer et al., 2001). While a Holling type IV function, which at high kelp densities, exhibits declining total (not just the proportion of) kelp consumption per urchin, has been used (Karatayev et al., 2019), a type III is preferred here because it allows both grazer switching and satiation. Support for the type III functional response is provided by the same factors that dampen the numerical response. At the predator-urchin level, I also use a logistic growth function for urchin density, and a Holling type III predation function to model the predator-urchin relationship. As for the urchin-kelp model, the arguments supporting a damped numerical response also support the Holling type III predation function.
The literature contains some ecological modeling, and some empirical work on the possibility of alternative stable equilibria (kelp forest versus urchin barrens) (see the Supplementary Appendix F). I ignore some of the details in the models mentioned in the literature, but nonlinearities, such as the recruitment facilitation posited by Baskett and Salomon (2010), or differing density dependence of predation across kelp beds versus barrens (Dunn and Hovel, 2019), are accommodated by the use of Holling type III grazing and predation functions (Ludwig et al., 1997). The long spined urchin invasion resulting from changed environmental factors (ocean warming) is modeled as a discrete increase in the maximum urchin density growth rate in the predator-urchin model.

Predator-Urchin Model
For the predator-urchin part of the model, with U as the urchin density,Uas its time derivative, and S as exogenous predator density (in g/m 2 ), a minimal dynamic model is given by Eq. 3.
Frontiers in Ecology and Evolution | www.frontiersin.org The first right hand side terms represent the natural growth of urchin density, with α and βas parameters. For the seaotter case the maximum growth rate, α, is set at 0.8. For the lobster case, environmental change of the warm East Australian Current has occasioned the invasion of long spined urchins. Treating all urchins as part of the same population, I follow the approach used for migration in population models (Preston and Wang, 2007), and model this as an exogenous increase in maximum urchin population growth. I give the maximum growth parameter, a, a pre-invasion level of 0.42, increasing to 0.66 and 0.8 as the invasion proceeds. The second term in Eq. 3 represents the loss in urchin biomass from sea otter or lobster predation, with γand θ as parameters. It is a sigmoid function that represents grazing by a given predator population, S, as a function of urchin biomass. Although the predator population is treated as fixed, its predation efficiency, defined as probability per predator of an urchin being consumed, or γU θ + U 2 , changes as the urchin density changes. When U = 0, predation efficiency is zero. When U > 0 but low, predation efficiency will be very low as predators will focus on other types of prey. At higher urchin densities, predation efficiency increases at a positive but declining rate, as predators shift more attention on urchin prey. Predation efficiency reaches a maximum, and then decreases as predators become overwhelmed by the higher prey density, as the nutritional quality of the prey decreases, or (for lobsters) as the spine length increases (Stewart and Konar, 2012;Eurich et al., 2014;Ling et al., 2015;Dunn et al., 2017;Dunn and Hovel, 2019). Watson and Estes (2011) also documented, at high urchin densities, an urchin flight response to evidence of a predation threat (urchin tests discarded by sea otters). To find the steady states I setU = 0 from Eq. 3. Since it is not expected that urchins will be completely extirpated in an area, I assume that, although it may get very small, U > 0. With exogenous S the isocline of steady state equilibria can be characterized as the real roots of the cubic Eq. 4. Figure 2 shows multiple isoclines for Eq. 4, with different values for α, with β = 0.08, θ = 1 and γ = 0.26. For given S and a, there can be either one or three real roots, and one or three equilibria. If there is only one equilibrium, it will be either on the upper or lower arm and locally stable. If there are three, there will be one on each of the three arms, U S,α h = U h (S, α) on the upper stable arm,U S,α l = U l (S, α)on the lower stable arm, and U S,α b = U b (S, α) on the unstable arm (see Supplementary Appendix A for local stability conditions).
The α = 0.8 isocline, the upper heavy weight curve, applies to the sea otter case. The vertical arrows indicate the inevitable forward and reverse shift thresholds for α = 0.8. If the initial equilibrium is a cryptic (low density) urchin equilibrium the resilience probability is described by Eq. 5, and if the initial equilibrium is an urchin barren (high density), by Eq the cryptic urchin (lower) arm of the α = 0.8 isocline of Figure 2. In the full resilience zone (S ≥ 8), perturbations from that equilibrium always return to the cryptic equilibrium. In the zero-resilience zone (S ≤ 5.5) the cryptic equilibrium does not exist. Instead, there is a fully resilient urchin barren. In the partial resilience zone (8 > S > 5.5), perturbations from the cryptic equilibrium are non-negative and bounded by the two stable equilibria. Assuming that perturbation outcomes from a cryptic urchin equilibrium in the partial resilience zone land at the originating equilibrium with a probability of ω l , or follow a uniform distribution, u U S,0.8 l , U S,0.8 h , with a probability 1 − ω l , the probabilistic resilience for ending up on the lower arm versus upper arm equilibrium is given by Eq. 5.
If, instead, the original equilibrium is an urchin barren, U S,0.8 h , on the upper arm of Figure 2, and the probability of a zero disturbance is ω h , the probability of return to U S,0.8 h is given by Eq. 6.
If 0 < ω l + ω h < 2, there is moderate clustering and some hysteresis, giving Q U S,0.8 In the partial resilience zone for both Eqn 5 and 6, 5.5 < S < 8, with hysteretic memory induced by clustering in the probability density function for the perturbationinitiated displacement location, the probability of a perturbation ending up at cryptic urchin/urchin barren equilibrium is greater if it originated from that equilibrium. There is only imperfect reversibility.

Urchin-Kelp Model
Turning to the urchin-kelp subsystem, I now treat the equilibrium urchin density from the predator-urchin subsystem as an exogenous parameter. With K as kelp biomass density,K as its time derivative, and U as urchin biomass, the minimal dynamic model is given by Eq. 7.
The first right hand side term of Eq. 7 represents the natural growth of kelp, with φ and λ as parameters. The second right hand side term represents the loss of kelp from urchin grazing with µ and ψ as parameters. Although some researchers have used a linear grazing function (Marzloff, 2012), others have proposed a sigmoid one . I use a sigmoid, or Holling type III, function because it best represents the urchin-kelp relationship described in most of the literature. With no numerical response, urchin biomass is treated as exogenous. However, the grazing efficiency of urchins (proportion of kelp being consumed per unit of urchin density), (µK) K 2 + ψ , changes as the kelp density changes. When K = 0, grazing efficiency per urchin is zero. When K > 0 but low, grazing efficiency is low as urchins focus on other resources, like detritus (Konar and Estes, 2003). With greater kelp density, grazing efficiency first increases at decreasing rate, reaches a maximum, and then declines toward zero. Increased urchin grazing efficiency at low kelp levels may be due to urchin recruitment facilitation by crustose coralline algae (a kelp competitor) (Baskett and Salomon, 2010) and/or urchins switching from detritus to living plants (Konar and Estes, 2003). At the highest levels of kelp biomass, there is evidence that the motion of kelp inhibits urchin grazing efficiency (Konar and Estes, 2003). From Eq. 7, assuming that although kelp density may get very small kelp will not be completely extirpated (K > 0), the K = 0 isocline contains the kelp-urchin combinations that are real roots of Eq. 8. As with Eq. 4, there are either one or three real roots.
Some research shows that kelp forests are extremely difficult to restore, there is very large amount of hysteresis in the urchin kelp system, and nearly all urchins may have to be eradicated for kelp forests to inevitably recover on barrens Marzloff et al., 2016;Eger et al., 2020).
To mimic these findings, I use most of the same parameter values as above (φ = 0.8, λ = 0.08, µ = 0.26), but with ψ = 0.01. Figure 1 shows the isocline with the arrowed lines indicating inevitable regime shift thresholds. With K h/U , K l/U , and K b/U as the upper, lower, and unstable arm equilibria. Local stability conditions are presented in Supplementary Appendix A for the predator-urchin model. The probability of return to an equilibrium on the high kelp arm is given in Eq. 1, and on the low kelp arm by Eq. 2. Fully resilient high kelp equilibria occur only if U ≤ 0.6. The partial resilience zone is 0.6 < U < 7.7, and in the range 7.7 ≤ U there are only low kelp equilibria, and zero high kelp resilience.
For a given predator density, S, and environmental state, α, over the two tropic relationships, depending on whether each relationship is in its partial resilience zone, there are up to four possible joint stable equilibria. For the predator-urchin system, partial resilience gives two urchin density equilibria, the low-density cryptic urchin equilibrium, U S,α l , and the highdensity urchin barren, U S,α h . Conditional on the cryptic urchin equilibrium, designated C, partial resilience in the urchinkelp system gives two kelp density equilibria, the high-density kelp forest,K S,α h/U l (KC), and the low-density kelp depauperate community, K S,α l/U l (DC). Conditional on an urchin barren, designated B, partial resilience for kelp gives either a low-density kelp depauperate community, K S,α l/U h (DB), or a kelp forest with lots of urchins, K S,α h/U h (KB). KB and DC are lower resilience versions of KC and DB, respectively. The KB can be characterized as a kelp forest with lower resilience because higher urchin density puts it closer to its inevitable tipping point to DB. It is an example of what Tracey et al. (2015) call an incipient barren. The DC is a kelp depauperate community with low resilience because the urchin density is low enough to put it close to its tipping point back to KC. All four joint equilibria exist if both the urchin and kelp equilibria are in their partial resilience zones. For S = 6 and α = 0.8, the equilibrium kelp density for KC is K 6,0,8 h/0.7 = 9.7, for DC is K 6,0,8 l/0.7 = 0.06, for DB is K 6,0.8 l/7.4 = 0.004, and for KB is K 6,0.8 h/7.4 = 5.9. Figure 2 shows the two urchin equilibria, and Figure 1 shows the conditional kelp equilibria.
Outside the zones of partial resilience, not all four joint equilibria exist. At very high or very low predator densities here may be full resilience for both the cryptic urchin and kelp forest (KC), or for both the urchin barren and depauperate kelp equilibria (DB). In each case there is only one joint equilibrium, a KC at very high predator density, or a DB at very low predator density. There is also the possibility of two partially resilient kelp equilibria, conditional on one fully resilient urchin equilibrium, yielding two joint equilibria, either DB and KB, or DC and KC. Finally, there can be two partially resilient urchin equilibria, one or both of which condition a fully resilient kelp equilibrium. Supplementary Appendix E has an example of the former, with three joint stable equilibria, KC, DB and KB.

Sustainable Natural Capital
Kelp Density, Resilience, and Hysteresis Kelp density, resilience, and hysteresis are all important to the measurement of sustainable natural capital. Although I use equilibrium kelp density as the resource input to the indicator of the sustainable stock of natural capital, it is conditional on urchin density, which is in turn conditional on predator density. The resilience probabilities and hysteresis for both the kelp and urchin equilibria, which describe their ability to persist, are the sustainability inputs.
Resilience probabilities are tied to pulse perturbations in both the urchin and kelp systems. If any joint equilibrium has full resilience in both the urchin and kelp systems, all perturbations return to it and the kelp at that equilibrium is sustained. If the initial joint equilibrium is a fully resilient KC, its kelp density or natural capital will be sustained at a high level. If that KC does not exhibit full resilience in both systems, the kelp density at the KC equilibrium is not sustainable, and a perturbation can cause a shift to another stable joint equilibrium with a lesser amount of natural capital. For a measure of sustainable natural capital, I consider multiple rounds of perturbations, giving the ability to look at multiple generations of shifts from an originating equilibrium.
For any given perturbation round, there is both an urchin density perturbation and a kelp density perturbation. I assume that the urchin density perturbation occurs first, and that urchin density has settled at its post-perturbation equilibrium before kelp density is perturbed. Kelp density adjusts to shifts in urchin density along the upper (kelp forest) or lower (kelp depauperate) arm of Figure 1 but then faces separate perturbations. A round of pertubations allows for transitions between the joint equilibrium states, with the resiience probabilites determining the transition probabilities. If a perturbation round starts at U S,α l , there is a probability of Q U S,α l/l of returning to that equilibrium following the round, and a probabilty of Assuming S = 6 and α = 0.8, other parameters as in Figures 1,  2, and with the resilience probabilities from Eqs 1, 2, 5, 6 with π h = π l = ω h = ω l = 0.5, Eq. 9 becomes Eq. 10.
Let ( KC0   DC0  DB0 KB0 ) be the initial vector of probabilities for the four possible joint equilibria, with KC0 + DC0 + DB0 + KB0 = 1. Starting from a KC equilibrium gives an initial probabilty of KC0 = 1. After one perturbation round, using the transition matrix in Eq. 10, the probabilties of each equilibrium are given by the right most vector in Eq. 11. With the equilbrium kelp densities for KC, DC, DB, and KB, expected kelp density after one round of perturbations from a KC equilbrium is Eq. 12. Expected kelp has declined from its pre-perturbation KC level of K 6,0.8 h/0.7 = 9.7 because partial resilience allows there to be non-zero probabilities of displacement to the basins of attraction of other equilibria with lower kelp densities.
If perturbation rounds continue, hysteresis will limit the ability to return to the KC equilibrium. Starting from KC, after two rounds the expected kelp density will be Eq. 13.
At some point there will be a stationary probability vector, which remains the same from one round to the next and depends only on the transition probabilities, not on the originating equilibrium. The equilibrium probabilities for this case are given by Eq. 14.
Using the equilbrium probabilities, there will be an equilibrium expected kelp density, Eq. 15.
The factors potentially affecting the expected kelp measures are: (i) kelp densities at the four stable equilibrium, (ii) the probabilty vector for the initial equilibrium, and (iii) the transition probability matrix, whose elements are affected by both the resilience of the four stable equilibria and the degree of hysteresis present. A more highly resilient KC equilbirum would mean a larger a KCKC in the northwest corner of the A matrix in Eq. 10 and smaller entries in the rest of the row. This would increase K 6,0.8 EKC1 , K 6,0.8 EKC2 , through to K 6,0.8 E * . Hystersis also influences the A matrix. Strong hystersis (π h = π l = ω h = ω l = 1) would produce an identity matrix for A, causing the probabilty vector for the initial equilibrium to be repeated after each round of perturbations, and expected kelp density to remain at the initial level. There would be no convergence to a common equilibrium. If the intitial equilibrium was a KC,K 6,0.8 h/0.7 = K 6,0.8 EKC1 = K 6,0.8 EKC2 = K 6,0.8 EKC * = 9.7. If it was a DB, K 6,0.8 l/7.4 = K 6,0.8 EDB1 = K 6,0.8 EDB2 = K 6,0.8 EDB * = 0.004. With more moderate hysteresis, there is convergence, but, depending on hysteresis strength, it may take a large number of perturbations rounds to achieve. The case in Eqs 14 and 15 takes around ten perturbation rounds. With no hysteresis, the probability of ending up at any equilibrium will be independent of the perturbing equilibrium, and the A matrix will have identical rows. The result is immediate convergence to K 6,0.8 E * = 5.1. Supplemental Appendix C shows the transition matrices and outcomes for the strong and no hysteresis cases. Whether the better measure of expected kelp density involves few or many rounds of perturbations, depends on the context in which it is being used.

Losses in Sustainable Natural Capital
Natural capital, as measured by expected kelp density, can be lost through a decrease in the predator density parameter, S, or an increase in the urchin growth rate, α. For the Haida Gwaii case, a decrease in S (or loss of sea-otters) caused the loss in expected kelp. The urchin growth parameter, α, remained fixed at 0.8, while the fur trade caused sea otter extirpation is modeled as reducing S from S ≥ 8 to S ≤ 5.5. For the Tasmanian case predator density remains fixed at S = 5.5, but α, reflecting the urchin invasion, increases from 0.42 to 0.8.
First, consider the Haida Gwaii case. The situation prior to the post-contact sea otter fur trade would have been a low urchin density equilibrium in which both the sea-otter population and kelp beds were largely intact and which exhibited a high degree of resilience (Sloan and Bartier, 2000;Sloan and Dick, 2012). In Figure 2, this suggests S ≥ 8. With S ≥ 8 and α = 0.8 there was a fully resilient low urchin equilibrium at U ≥8,0.8 l ≤ 0.4 (Figure 2), and with U ≥8,0.8 l ≤ 0.4 there was a fully resilient high kelp equilibrium at K ≥8,0.8 h/≤0.4 ≥ 9.8 (Figure 1). Since both the low urchin and high kelp equilibria were in their full resilience zones, there was only a KC equilibrium, with A = a KFKF = Q U >8,0.8 l/l = P K ≥8,0.8 h/h/≤0.4 = 1. This ensured sustainable expected kelp, K ≥8,0.8 h/≤0.4 = K ≥8,0.8 EKC1 = K ≥8,0.8 EKC2 = K ≥8,0.8 EKC * ≥ 9.8. The fur trade reduced S below the reverse shift threshold of S = 8, and, as that happened, the low urchin equilibrium lost some if its resilience to perturbations. In the range 8 > S > 5.5, large enough perturbations from the cryptic urchin equilibrium would have had the potential to land in basin of attraction of the urchin barren equilibrium and cause a move to that equilibrium.
For example, if S was reduced to S = 6, the initial equilibrium would have been a KC with a slightly lower kelp density, K 6,0.8 h/0.7 l = 9.7. But, due to a lack of resilience and hysteresis, π l = π h = ω l = ω h = 0.5, there would have been a much lower sustainable expected kelp density, K 6,0.8 EKC1 = 7.32 > K 6,0.8 EKC2 = 5.73 > K 6,0.8 EKC * ≥ 3.29. Further decreases in sea otter density past the forward shift threshold of S = 5.5 would have caused low urchin resilience to be completely lost, the low urchin equilibrium to disappear, and would have made the shift to an equilibrium on the high urchin arm inevitable. In the urchin-kelp system, the urchin shift would have caused an inevitable shift to the lower kelp arm. Once S had been reduced to S ≤ 5.5, there was only fully resilient DB equilibrium, ≤ 0.004. Further decreases in S resulted in further gradual increases in urchin density along the high urchin arm and gradual decreases in kelp density along the lower kelp arm. There the system has remained (Sloan et al., 2001).
For the Tasmanian lobster-urchin case, the pre-invasion locus assumes α = 0.42, the lowest (dashed) isocline in Figure 2. This isocline does not bend back on itself, leaving no possibility for an abrupt regime shift or hysteresis. With S = 5.5, U 5.5,0.42 / = 0.4, and K 5.5,0.42 h/0.4 = 9.9 (Figures 1, 2), with all equilibria stable and fully resilient. The long spined urchin invasion transformed the isocline in Figure 2 from the dashed α = 0.42 locus, to the normal thickness α = 0.66 locus and is causing it to approach the heavier weight α = 0.8 locus. With S constant at S = 5.5, the result of the shift from α = 0.42 to α = 0.66 was a move from a fully resilient KC equilibria at U 5.5,0.42 / = 0.4 , K 5.5,0.42 h/0.4 = 9.9 to another fully resilient equilibrium just inside the boundaries of the full resilience zones for both urchin and kelp at U 5.5,0.66 / = 0.6, K 5.5,0.66 h/0.6 = 9.8. Once α > 0.66, the equilibria moved into the zones of partial resilience, and there was the possibility for a perturbation to kick the urchin equilibrium from the lower arm to the upper arm, and/or the kelp equilibrium from the upper arm to the lower arm. As shown in the Supplemental Appendix D, α = 0.7, arrived at via increasing a, gave a partially resilient KC equilibrium at K 5.5,0.7 h/0.64 = 9.75. With perturbations clustered according to π l = π h = ω l = ω h = 0.5, hysteresis results in lower expected kelp, K 5.5,0.7 EKC1 = 8.37,K 5.5,0.7 EKC2 = 7.7, and K 5.5,0.7 EKC * = 6.5. When α passes 0.8, there will be only a high urchin equilibrium,U 5.5,≥0.8 h ≥ 7.7, and a low kelp equilibrium, K 5.5,≥0.8 l/≥7.7 ≤ 0.004. Both will be fully resilient, giving a DB equilibrium with A = a DBDB = 1 and sustainable natural capital of K 5.5,≥0.8 EDB * ≤ 0.004.

Management Controls
Systemic, Symptomatic, Press, Pulse, Preventative, and Restorative Controls I consider management controls in the context of the two cases, Tasmania and Haida Gwaii. Management controls are necessary to prevent or reverse losses in sustainable natural capital. These controls may be systemic, symptomatic, press, pulse, preventative, or restorative (Scheffer et al., 2001;Lessard et al., 2005;Layton et al., 2020). Systemic controls aim to remove the ultimate causes of the losses. They are typically ongoing press controls (like permanent reductions in the harvesting of predators). In contrast, symptomatic control implies using leverage points for more direct control of affected populations, urchin harvesting or culling, or kelp enhancement (seeding or transplantation) (Lessard et al., 2005). Symptomatic controls may be either pulse (like one time or intermittent urchin culling in response to perturbations in urchin density) or press. Preventative controls aim to prevent a loss of kelp forest, and restorative controls to restore it after the loss has occurred. See Supplementary Appendix A for background on these controls.
In the Tasmanian case, the ulimate cause or threat is ocean warming, which for scientific and governance reasons, is not easily removable ). Harvesting of the lobster predator is an exacerbating factor , but there are questions as to whether restrictions on lobster harvest, be they systemic or symptomatic press controls, can be sufficent to prevent significant kelp forest losses (Tracey et al., 2015). Starting from a KC equilibrium, if the long spined urchin invasion proceeds to α = 0.8 in Figure 2, imposing lobster harvesting restrictions to ensure S > 5.5, will keep the equilibrium in the partial resilience zone. But, in that zone, natural perturbations in urchin and/or kelp density have the potential to initiate a shift away from the KC equilibrium.
In the Haida Gwaii case, the cause of the decrease in sustainable natural capital was overharvesting and extirpation of the sea otter predator (Sloan and Dick, 2012). The postextirpation equilibrium is effectively a fully resilient DB. Systemic press control would have to produce a significant increase in effective sea otter density (supplemented by press uchin harvesting or culling to mimic sea otter predation) to move to the partial resilience zone, and a likely infeasible increase to get to a fully resilient KC equilibrium. While the sea otter has been reintroduced in the Pacific Northwest, only recently have there been indications that it will return to Haida Gwaii (Gwaii Haanas Archipelago Management Board, 2020). Press human harvesting of urchins is difficult or costly, and the pending return of the sea otter is not necessarily encouraged by human harvesters of other sea otter prey like abalone (Lee et al., 2016). Systemic control seems unlikely to be successful on its own.
While systemic or symtomatic press controls may be necessary to move into, and remain in, the zone of partial resilience (Layton et al., 2020), in that zone pulse controls may provide a cost effective way to cushion or enhance perturbations, disuading or encouraging a kick to an alternative equilibrium. Since both the Tasmanian and Haida Gwaii probrems are difficult to solve with press controls alone, pulse controls are useful complements that are designed to produce one time or periodic changes in the transition probability matrix. For simplicity of exposition, I illustrate with changes that increase crucial transition probabilities to one.
In the Tasmanian case, suppose press lobster harvesting restrictions have increased S to S = 6 while the long spined urchin invasion has produced α = 0.8. Assuming the KC equilibrium still prevails, KC0 = 1. Even if press lobster harvest controls can maintain S = 6 in the partial resilience zone, perturbations still have the potential to kick the system to other equilibria, including a DB equilibrium. Over many perturbation rounds sustainable natural capital would be reduced to K 6,0.8 EKC * = 3.29. However, supplementing press lobster harvest controls with others that are symptomatic, pulse and preventative may be able to sustain natural capital at K 6,0.8 EKC * = 9.7. These would be periodic uchin culling, perhaps supplemented by periodic kelp enhancement (Flukes et al., 2012;Tracey et al., 2015). Pulse urchin culling could be used to counteract positive urchin perturbations from the low urchin equilibrium in both fthe first and subsequent rounds. ideally increasing Q U 6,0.8 l/l from 0.6 to 1. Pulse kelp enhancement could counteract and negative perturbations from the hign kelp equilibrium in multiple rounds, ideally raising P K 6,0.8 h/h/0.7 from 0.99 to 1, and producing an increase in overall transition probability, a KCKC = Q U 6,0.8 l/l P K 6,0.8 h/h/0.7 , from 0.594 to 1. This changes the transition probabilty matrix from Eq. 10 to A in Eq. 16. With culling and kelp enhancement used to offset successfully all urchin and kelp perturbations, from KC, K 6,0.8 h/U l = K 6,0.8 EKC1 = K 6,0.8 ELC2 = K 6.0;8 EKC * = 9.7, and natural capital would be preserved.
Turing to the Haida Gwaii case, if effective sea otter density can be increased to S = 6, the intitial equilbrium is a DB in the partial resilience zone. With no further action, other than maintaining S = 6, K 6,0.8 l/7.4 = 0.004 < K 6,0.8 EDB1 = 1.33 < K 6,0.8 EDB2 = 2.17 < K 6.0;8 EDB * = 3.29. The restoration of a kelp forest would required pulse urchin culling to create or enhance negative urchin density perturbations, possibly paired with pulse kelp enhancement. Starting at a DB equilibrium ( DB0 = 1), urchin culling could ideally increase Q U 6,0.8 l/h from 0.1 to 1, and pulse kelp enhancment would ideally create positive perturbation in kelp density to increase P K 6,0.8 h/l/0.7 from 0.49 to 1. Overall the increase in a KCDB = Q U 6,0.8 l/h P K 6,0.8 h/l/0.7 would be from 0.049 to 1. The transition probabilty matrix becomes A " in Eq. 17.
The restoration involves a shift from a DB to a KC. Since S remains in the partial resilience zone, there would need to be pulse preventative culling and enhancement to maintain the KC equilibrium. Maintaining a kelp forest equilbrrium is not as hard as restoring it, requiring only an sustained increase in a KCKC from 0.549 to 1 as in Eq. 16. Using restorative control followed by ongoing, but intermittent, preventative control, expected kelp density will increase from K 6,0.8 l/7.4 = 0.004 < K 6,0.8 EDB1 = 1.33 < K 6,0.8 EDB2 = 2.17 < K 6.0;8 EDB * = 3.29 to K 6,0.8 h/7.4 = K 6,0.8 EKC1 = K 6,0.8 ELC2 = K 6.0;8 EKC * = 9.7.

Preventative Control as Adaptive Management
Although sometimes restorative controls are necessary, when there is a choice preventative controls are usually preferable to restorative ones (Young, 2000). I show this in the context of pulse controls, although similar aguments can be made for press controls. There are three arguments to be made. First, due to hysterisis, getting back to a KC from an equilibrium like a DB, involves a larger transition probability increase than preventing the loss of the KC. Second, the cost function for transition probability increases (via urchin culling and kelp enhancement) may be strictly convex, with the increasing marginal cost for the probabilty increases required to restore a KC, making resoration of badly degraded barrens prohibitively costly (Layton et al., 2020). Third, a control which prevents a shift away from a KC can be used as adaptive management, providing more flexibilty for subsequent controls than doing nothing and allowing the shift.
In the presence of current uncertainty about whether innovations will bring down restoration costs (Sunnset et al., 2010;Fredriksen et al., 2020), research to reduce the uncertainty, combined with the flexibility to resort to the option of restortative control only if it can be done cheaply and effectively, can be worthwhile. The extra benefit provided by using preventative control as adaptive management is its option value.
Using the Tasmanian case with S = 6 and α = 0.8, I first compare a pure preventative control plan (three round prevention controls) versus restorative control plan (after a first round of prevention there is a round of no control followed by a round of restoration control). Starting at a KC equilibrium, KC0 = 1, both plans put us back at a a KC after three rounds, KC3 = 1. This allows us to make comparisons on the basis of the first three rounds alone. With the initial probability of a KC equilibrium at KC0 = 1 and Eq. 10 as the no control transition matrix, the pure preventative plan requires increasing the probability of a return to the KC from 0.594 to 1 in round one. This changes the transition matrix from Eqs 10-16. Continuing preventative control sustains the KC equilibrium.
The restoration plan also starts at a KC equilibrium, KC0 = 1. With preventative control in round one, Eq. 16 is the transition matrix and the outcome is KC1 = 1. With no control on the second round, the transition matrix reverts to A in Eq. 10, and the outcome vector is ( KC2 DC2 DB2 KB2 ) = (0.594 0.006 0.136 0.264). Expected kelp after a round of perturbations with no control drops toK 6,0.8 EKC1 = 7.32. Restoration occurs in round three, with the goal that KC3 = 1, with K 6,0.8 h/0.7 = 9.7. To generate this outcome, restoration must change the transition matrix from A to A in Eq. 18, with a left-hand column of ones.
Now consider the adaptive plan that is useful when there is uncertainty about an important decision factor like restoration costs. It is structured to allow for adaptive management and the creation of option value. It uses research to reduce uncertainty and a decision structure that has the flexibility to take advantage of the reduction in uncertainty. The first-round choice is crucial. To be useful research must be undertaken early, and the firstround choice should be one that is not difficult to change in the second round. For example, doing nothing in the first round may make it impossible to use prevention in the second round. Successful kelp forest restoration would have to be undertaken first. However, prevention in the first round can easily be changed to doing nothing in the second round. The adaptive plan combines first round preventative control with research which, by the beginning of the second round, reveals whether restoration costs are strictly convex (high) or linear (low). Once the restoration costs are known, the second-round decision, to continue prevention in rounds two and three (like the pure prevention plan) or to do nothing in round two and pursue restoration in round three (like the restoration plan), both followed from then on by prevention, can be taken.
Like the other two plans, the adaptive plan starts with KC0 = 1. In round one, prevention changes Ain Eq. 10 to A in Eq. 16, producing KC1 = 1. With research into restoration costs yielding cost certainty by the beginning of round two, at this point there is a choice to be made between continued prevention versus a sequence of no control followed by restoration. The former maintains the transition probabilities as A and yields KC2 = KC3 = 1. The latter goes back toA for the second round yielding ( KC2   DC2  DB2 KB2 ) = (0.594 0.006 0.136 0.264), and a switch from A to A in (??) is required for restoration to KC3 = 1 in the third round. To compare the costs of the three plans, it is necessary to model cost function uncertainty. I assume the cost of increasing a transition probability, a KCKC , a KCDC , a KCDB , or a KCKB , to certainty, through a combination of pulse urchin culling and kelp enhancement, has a 0.5 probability of being given by a strictly convex cost function and a 0.5 probability of being given by a linear one. With KC as the original equilibrium, the convex and linear functions are C c = 6.06 (1 − a KCKC ) 3 and C n = (1 − a KCKC ), respectively. For the smaller transition probability increase required for preventative control, an increase in a KCKC from 0.594 in A in Eq. 10 to 1 in A in Eq. 16, uncertainty with respect to the cost function doesn't matter. Figure 3 shows 1 − a (without the subscript) on the horizontal axis versus the cost, C. As shown in Eq. 19 and Figure 3, with 1 − a = 0.406 the prevention cost is C P = 0.406 (at point b) regardless of the shape of the cost function. The per round costs of prevention are the same regardless of whether the cost function is strictly convex or linear.
However, with restoration from A in Eq. 10 to A in Eq. 18 it does matter whether the cost function is strictly convex or linear. Other equilibria, such as a DB, have positive probability ( KC2 DC2 DB2 KB2 ) = (0.549 0.006 0.136 0.264), and all three of them require a larger probability increase to restore a KC than is required to maintain it. For the strictly convex cost outcome, the expected one-time restoration cost is Eq. 20.
With a 0.5 probability of each of the two cost function outcomes, one-time expected restoration costs are E C R = 0.5 C R c + C R n = 1.003. Figure 3 shows C P = 0.406, C R n = 0.559, C R c = 1.446, and E C R = 1.003 as b c, d, and e, respectively. Now consider the cost of the three alternative plans for the Tasmanian case, starting from a KC and considering three potential rounds of perturbations. I use no discounting, as it would not affect the nature of the results, only their magnitude. The goal is to minimize the cost of returning to KC by the end of the third round. If there is no control during a round, the cost includes losses of services from the temporary loss of natural capital.
The pure preventative policy is a continuing preventative policy that always brings the system back to a KC equilibrium. Over three perturbation rounds, the control cost is C PPP = 3C P = 1.218, with natural capital sustained at K 6,0.8 EKC * = 9.7. The restorative policy is to revert to no control for the second round of perturbations, letting natural capital temporarily drop from K 6,0.8 h/0.7 = 9.7 to a level of K 6,0.8 EKC1 = 7.32, for a lost value of services in round two of L. In round three, restoration to the kelp forest state is undertaken at an expected cost of E C R = 1.003 and there is a return to a sustained K 6,0.8 h/0.7 = K 6,0.8 EKC * = 9.7. The three-round expected cost of the restorative policy is Eq. 22.

C PLR
= C P + L + E C R = 0.406 + L + 1.003 = L + 1.409 (22) Comparing the costs of the purely preventative and restoration plans gives C PLR − C PPP = L + E C R − 2C p = L + 0.191. In Figure 3, L = 0, but E C R = 1.003 (point e) is high enough to exceed 2C P = 0.812 (double the height of point b). The pure prevention plan is clearly more cost effective. The crucial difference is that E C R > C P by a substantial amount. The two factors accounting for this difference are the greater probability increases required for restoration, and the increasing marginal cost in the event that C R c is the cost function outcome. A large loss in natural capital services, L, would mitigate even more strongly in favor of the preventative control policy. Now consider the adaptive plan. It uses research to reduce uncertainty and a decision structure that has the flexibility to take advantage of the reduction in uncertainty. The first-round choice is crucial. The adaptive plan combines first round preventative control with research which, by the beginning of the second round, reveals whether restoration costs are strictly convex (high) or linear (low). Once the restoration costs are known, the second-round decision can be taken, to continue prevention in rounds two and three (like the pure prevention plan) or to do nothing in round two and pursue restoration in round three, (like the restoration plan). Including research costs, W, the firstround cost is C P = 0.406 + W. In the second round, the choice involves using the flexibility built into the preventative control to minimize costs for the second and third rounds based on newly acquired certainty about restoration costs. Over three rounds the costs for the adaptive plan are given by Eq. 23.
The difference, V = C PPP − C PAA = 1.218 − C PAA , is the option value generated by the combination of knowledge about restoration costs and flexibility to change control options. Whether V > 0 depends upon whether the cheap restoration option, including the lost natural capital services and the cost of research, is less than the cost of continuing prevention. In Eq. 23 this requires 0.5L < 0.1265 − W. If 0.5L < 0.1265 − W, knowing after the first round that there will be cheap restoration costs, makes it worthwhile to do nothing in round two and then undertake restoration, giving C PAA = 1.0915 + W + 0.5L and an option value, V = C ppp − C PAA = 0.1265 − W − 0.5L. In Figure 3, L = 0 and V = 0.1265 − W. The possibility for a positive V comes from being able to choose one later round of cheap restoration (point c in Figure 3 with a restoration cost of 0.559) if it ends up being available (0.5 probability), rather than two more rounds of prevention costs for a total of 0.812 (double the height of point b in Figure 3) if cheap restoration is not available (0.5 probability). As long as L and the cost of research, W, are not too great, V > 0 and research to ascertain restoration costs is worthwhile. The adaptation plan provides an extra argument, based on flexibility, for current use of preventative controls. The overall conclusion is that initial prevention is always preferred. If it is paired with research that is not too costly, and it turns out that restoration costs, including lost natural capital services, are small, prevention can be abandoned in favor of later restoration. Otherwise prevention should continue to be used as needed.

Adaptive Manegement as Part of Restorative Control
For the Haida Gwaii case, the system is better characterized as already in a fully resilient DB state, with restorative control the only alternative to shift back to a KC. In a fully resilient DB state, all perturbations return to the DB, and pulse controls alone will not be helpful in producing the shift. Until the sea otter predator returns in sufficient numbers, there must also be ongoing press urchin culling control, simulating sea otter predation and pushing the system back at least into the partial resilience zones, where pulse urchin culling, with or without kelp enhancement, can kick the system back to a KC. Whether a given amount of press urchin culling will move the effective S from a fully resilient urchin barren into the partial resilience zone depends on the rate of sea otter recovery, and the thresholds at which the urchin density shifts from high to low and kelp density shifts from low to high, all of which are uncertain (Goldman, 2020). Here I simplify by characterizing the uncertainty as doubt about how much resilience the urchin barren equilibrium will lose for a given amount of press culling. Hypothetically, assume there is a 0.5 probability that otter recovery, combined with some fixed amount of press urchin culling, will increase S to S = 6, with a DB equilibrium and a transition probability a KCDB = 0.049 (from Eq. 10) of ending up at KC following negative urchin and positive kelp pulse perturbations. Or there is a probability of 0.5 that press culling alone will increase S to S = 9, with a KCDB = 1, a move to a fully resilient kelp forest equilibrium. Without strong monitoring of either sea-otter numbers or kelp and urchin abundance (Estes et al., 2010;Williams et al., 2018), we will not know whether we are approaching S = 6, requiring a pulse of urchin culling and kelp enhancement to get to a KC, or S = 9, requiring none. Instead, we know only that we have fifty-fifty chance of approaching either S = 6or S = 9 for an expected predator density of E(S) = 7.5. Adaptive management can be beneficial here as well. Knowing whether S = 6 or S = 9 enables the policy maker to tailor investment in pulse culling and enhancement control to the sea-otter recovery level, with the possibility of making that investment more cost effective and generating positive option value. Details are presented in the Supplementary Appendix D.

DISCUSSION
Addressing the loss of resilience and depletion of kelp forests as a foundation species in an important marine ecosystem, this paper developed a mathematical model to analyze catastrophic shifts in urchin and kelp densities. The recursively decoupled catastrophe model has two trophic interactions, with the potential for two catastrophic shifts, one at the level of the predator-urchin trophic interaction and one at the level of the urchin-kelp trophic interaction. It incorporates the concept of resilience at both levels and uses resilience probabilities and hysteresis to describe the probabilities of shift patterns between alternative urchin and kelp equilibria and as part of a definition of a sustainable natural capital stock as expected kelp density. To prevent or reverse losses in the sustainable natural capital stock, I consider systemic versus symptomatic, press versus pulse, and preventative versus restorative controls. There is often complementarity in that combinations of controls can be more successful than one type alone. Pulse controls alone cannot be successful in converting a DB to a KC unless the originating equilibrium is in the partial resilience zone. But they may provide a cost-effective complement to a systemic or symptomatic press control that moves the originating equilibrium into, and keeps it in, that zone. With hysteretic memory activated by clustered perturbations, and convex cost function for controls like urchin culling and kelp enhancement, it is likely to be more difficult and costly to restore a kelp forest, than to prevent its loss. I also suggest a role for adaptive management, in that improving knowledge about such uncertain factors as restoration cost or the rate of sea otter recovery can, if control policies are flexible, create an option value that reduces the cost of maintaining or restoring a kelp forest state.
The mathematical model has some limitations. The pure topdown predator-urchin-kelp cascade is often used in the literature and makes the modeling more straight forward, but also ignores some ecological complexity (Wallington et al., 2005). This may limit application to other systems, particularly those where the downward cascade is less prominent. Exogenous factors, such as climatic changes or ongoing pollution, can be incorporated by parameter manipulation, as was done for the invasive long spined sea urchin, but complexity may be compromised here as well (Mooney and Cleland, 2001;Karatayev and Baskett, 2020). Because it generates multiple stable equilibria, the Hollingtype III predation function plays a critical role in the model. While findings of alternative stable states with different shift thresholds (hysteresis) support type III (Ling et al., 2019), direct empirical estimates of functional form do not always do so (Dunn and Hovel, 2019;Karatayev et al., 2019). Another crucial assumption is the treatment of pulse perturbations and hysteresis. While pulse perturbations clearly have a role in initiating shifts between equilibria, and in creating hysteresis, there has been little or no incorporation of them into mathematical models. The piecewise approach is an initial attempt. For both predation functions, and pulse perturbations, future model adjustments are possible depending on what scientific evidence suggests. Other case studies, such as the California marine heat wave that combines climate perturbations with predator decline (Rogers- Bennett and Catton, 2019;McPherson et al., 2021), could be usefully explored. Finally, while sustainable natural capital is a clear enough concept, both physical and value measures are mostly missing (Loomis, 2006). Overall, improved empirical evidence, both ecologic and economic, is needed to give more substance to the final verdict. Nevertheless, despite its limitations, the model is based on available scientific evidence and does generate important implications for control policies. Pairing with the economic concepts of sustainable natural capital and the cost functions for controls, enriches the policy conclusion.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.