Alternative stable equilibria and critical thresholds created by fire regimes and plant responses in a fire‐prone community

55 –––––––––––––––––––––––––––––––––––––––– © 2018 The Authors. This is an Online Open article This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Erica Fleishman Editor-in-Chief: Miguel Araújo Accepted 27 May 2018 42: 55–66, 2019 doi: 10.1111/ecog.03491 doi: 10.1111/ecog.03491 42 55–66


Introduction
Disturbance plays critical roles in shaping plant communities (Levin and Paine 1974, Connell 1978, Collins 2000, Roxburgh et al. 2004, Turner 2010, and aspects of disturbances can act and interact to affect competitive outcomes and plant diversity Alternative stable equilibria and critical thresholds created by fire regimes and plant responses in a fire-prone community (Garrison et al. 2012, Miller et al. 2012a, Seifan et al. 2012, Liao et al. 2016. Global change may alter both disturbance regimes (Dale et al. 2001, Barbero et al. 2015, Abatzoglou andWilliams 2016) and species' responses to disturbance (Zhang et al. 2012, Anderson-Teixeira et al. 2013. In some cases, global change pressures -e.g. non-native invasive species, increased fire ignition rates, or climate changemay affect disturbance-recovery dynamics, producing irreversible shifts in forest cover over extensive areas (D'Antonio and Vitousek 1992, Pausas 2015, Tepley et al. 2016. Theoretical understanding of such dynamics is both of basic ecological interest and needed for informing predictions of climate change effects and potential mitigation strategies. Some environments may support plant communities dominated by either fire-resistant tree species or pyrogenic shrubs or grasses, where self-reinforcing feedbacks could perpetuate either community type as alternative stable equilibria (a.k.a. alternative stable states, D'Antonio and Vitousek 1992, Tepley et al. 2016). Alternative stable equilibria is a condition in dynamical systems where distinct equilibrial states are locally stable, and thus the long-term fate of the system cannot be generally predicted. If forest and shrub states are both stable, then the system can converge on either community, depending on both initial conditions and fire stochasticity.
When a parameter change alters stability of equilibria such that the system transitions from having two alternative stable equilibria to one, the system can have hysteresis. Hysteresis is a phenomenon wherein stability change due to small perturbations of a parameter sends the system to a very different equilibrium, and this change cannot be reversed by returning the parameter to its previous value (Lindenmayer et al. 2011, see also Fig. 1). In ecological communities, hysteresis results in thresholds that, when crossed, cause abrupt shifts in biomass and biological diversity that cannot be readily be reversed (May 1977, Scheffer et al. 2001, Beisner et al. 2003, Scheffer and Carpenter 2003. There are many possible types of thresholds in dynamical systems, and often the system is free to move back to a starting equilibrium after a parameter whose shift drove ecosystem change returns to its original value. In contrast, the term 'critical threshold' is used to specifically convey a point or line where a critical transition (Kuehn 2011) -i.e. a change that will not be reversed even if all parameters return to their original state (Suding andHobbs 2009, Wang et al. 2012). Understanding what types of interactions support bistable dynamics and hysteresis may inform ecosystem management and conservation plans (Knowlton 2004, Scheffer et al. 2015. In the Klamath Mountains of northern California and southwestern Oregon, climate change is expected to alter disturbance-recovery dynamics with potentially substantial consequences for forest cover (Tepley et al. 2017(Tepley et al. , 2018. Ecosystems of the Klamath are affected by a highly variable, mixed-severity fire regime (Halofsky et al. 2011). This roughly four million-hectare region contains steep climatic gradients and lies at the transition zone between two regions. The Cascades and Coast Range Mountains Figure 1. Conceptual dynamics leading to alternative stable equilibria and critical thresholds. Left: general schematic of hysteresis due to climate (adapted from Scheffer et al. 2001). Changes to parameters can change basins of attraction and stability of equilibria, and the system can change states in a non-reversible manner. Right: proposed dynamics of Klamath cover types, in terms of biomass, fire intensity, and forest growth rate. Frequent, high mean intensity fires are hypothesized to result in low mean biomass shrubs with low-amplitude cycles, whereas frequent, low mean intensity fires are hypothesized to lead to high biomass conifer forest with higher amplitude biomass cycles. Thick solid lines indicate stable equilibria; regions are shaded to indicate three zones of stability. Open-headed arrows show typical trajectories, dotted arrows with solid heads show hysteretic, non-reversible change in equilibrium when a parameter changes dynamics from alternative stable equilibria to shrub-or forest-dominated. of Oregon are wetter, fire is rare, and contiguous closedcanopy forest is common. To the south, the region is dryer, frequently burned and chaparral-dominated. The Klamath can support either of two natural vegetation types (i.e. plant functional types; PFTs): shrub and hardwood vegetation (henceforth shrub or shrubland) or mature conifer forest (henceforth conifer) (Whittaker 1960, Odion et al. 2010. Mature conifer forest commonly persists following moderate intensity fires, but the shrub community is both conducive to and perpetuated by high-intensity fire. Mature conifers have adaptations such as thick bark and high crown base heights that allow them to survive low to medium intensity surface fires. Shrub communities are low growing, with fine fuels that enable fires to kill aboveground biomass. They are perpetuated by fire due to adaptations that allow them to rapidly reestablish after fire (Thompson and Spies 2009), including resprouting of hardwood species (Donato et al. 2009) and seedbank persistence (Knapp et al. 2012, Yelenik et al. 2013). This postfire recruitment possibly inhibits conifer establishment (Hobbs et al. 1992) and increases the spread and severity of future fires (Thompson et al. 2007). Through a cycle of frequent severe fires, shrubland may be able to persist for multiple generations (Odion et al. 2010), and would thus be considered stable (Connell and Sousa 1983).
Here, we investigate how the frequency and intensity of fire interact with plant traits such as fire resistance to affect the stability of cover types. We use the Klamath as a case study to construct a theoretical model of plant competition, disturbance, and recovery. Our model is flexible and can be applied to many fire-prone communities.
We use invasion analysis to assess the stability of community types. This method is more commonly applied to evaluate species coexistence (Chesson 1994, Miller et al. 2011), but it can be used to determine whether alternative stable equilibria are present and hysteretic dynamics are possible. Our approach provides an integrated notion of stability in the context of disturbance, using a multidimensional characterization of fire regimes and complex interactions between PFTs. We apply this model to the Klamath Mountains to elucidate how fire regimes can interact with plant life-history traits to produce bistability and critical thresholds. We hypothesize that stochastic cycles of regrowth can occur in three different stability regimes, depending on features of the fire process and plant traits. Two regimes are shrub dominated and conifer dominated, where one community equilibrium is stable and the other is unstable. In this case, all initial conditions lead to the same community type. The third stability regime, alternative stable equilibria, occurs when the equilibria for conifer and shrubland are both stable (Fig. 1). We conduct stability analysis to determine whether alternative stable equilibria occur, and we use the results and insights gained from the analysis to infer how climate change or management might affect the stability of different community types.

Model design
Our model is designed to capture key processes and interactions that affect plant communities in a forest-shrub system influenced by fire (Fig. 2). Although our model can be used to study how single PFT types react to disturbance depending on life history traits, here we use it to develop our theoretical understanding of alternative stable equilibria. It has been suggested that alternative stable equilibria theory can benefit from including 'more realism into simple models' (Scheffer and Carpenter 2003), and we construct our model in this spirit. Compared to models such as the ecosystem demography model (Medvigy et al. 2009, Miller et al. 2016 or LANDIS-II (landscape disturbance and succession II, Thompson et al. 2011a), our model has fewer input parameters and environmental variables. However, it features feedbacks and interactions that are complex relative to previous models of alternative stable equilibria (Scheffer 1989, Shurin et al. 2004, Mumby et al. 2007, Staver et al. 2011. Our model highlights environmental fluctuation, which has not often been included in alternative stable equilibria theory (Scheffer and Carpenter 2003) and we formally characterize stability and alternative stable equilibria under the influence of specific, quantifiable aspects of disturbance regimes.

Klamath features and parameterization
We set parameter values to represent key aspects of the fire and plant ecology of the Klamath (Table 1). Parameter selection was informed by differences in PFT characteristics, and we designed structural features of the model to capture fire and PFT interactions (Supplementary material Appendix 2 Table A1). For example, maturation time of the shrub PFT is set for two years, whereas conifer matures and gains strong fire resistance after 40 yr. The conifer community has high mean germination rate (0.9), and low seedbank survival (95% of non-germinating seeds die after one year). To capture the long-lived seedbank of shrub-hardwood PFTs (and the ability of some to resprout), their mean germination rate is lower (0.4), and seedbank survival is very high (0.01% loss per year). We do not attempt data-driven parameterization, because our goal is to examine features how alternative stable equilibria and critical thresholds arise and change due to environmental or PFT parameters. Our model is designed for theoretical inquiry and although it is conceptually possible to compute estimates for all parameters, most of them have not yet been estimated and published. For example, in our model, germination and fecundity parameters vary randomly to account for environmental variation over time (due, e.g. to precipitation, light, and temperatures), but this variation has not been estimated 4 empirically for most species. Additionally, parameters can be altered to consider different types of PFTs in different communities. For example, mortality could be set to one to represent annual species, or variation in fecundity could be set high to represent species that reproduce sporadically.
Thus, different life history traits can be encoded, and more life stages can be considered with minor modifications. To perform stability analysis for n > 2 PFTs, a stable distribution of size classes for each grouping of n -1 would need to be computed (Chesson 1994).

Temporal structure
Each year, an environmentally dependent random germination rate is set for each PFT. Upon germination, the seedlings compete for free space in which to establish, so each PFT secures space proportional to its representation in the seedling pool. Germinating seeds, and seeds that lose viability are removed from the seedbank. Established seedlings mature at a fixed rate. The density of mature individuals and fire intensity are used to compute the mortality of seedlings and mature trees, which controls free space available for new recruitment from the seedbank. We define intensity as the magnitude of the disturbing force (Shea et al. 2004), which in the case of our model represents the effects of factors such as weather (e.g. temperature, relative humidity, wind speed) or weather-driven fuel moisture.

Invasion analysis
We assess the stability of each PFT via an invasion analysis (sensu Bolker and Pacala 1999), wherein the long-term competitive outcome of disturbance and competition is determined via a long-term, low-density growth rate (Turelli 1978, Chesson 1994, Bolker and Pacala 1999, Adler et al. 2006. We utilize invasion analysis to test whether alternative stable equilibria occur. We model the dynamics of each PFT, considered as a resident. Then we compute the population growth rate of the other PFT as an invader, introduced at low density, that must compete with the resident for resources. The invader's growth rate varies depending on resident and invader densities, environmental variation, and the fire process. If each PFT has mean growth rate greater than one, both PFTs will grow from low density and persist indefinitely. In this case, stable coexistence as the outcome. If one PFT has a growth rate greater than one, then that PFT will exclude the other. Of special interest here is the case in which each PFT has a long-term, low-density growth rate less than one. In this case, the zero-density state is a stable equilibrium for each PFT, and the community types constitute two distinct equilibria with non-global basins of attraction, i.e., alternative stable equilibria. Which stable equilibrium is reached in a specific case of alternative stable equilibria depends on initial conditions and the environmental processes, including stochastic variation in fires, and variation in fecundity and germination. In this case, the system has alternative stable equilibria. If parameter changes move the system from having alternative stable equilibria to one globally stable equilibrium, then hysteresis occurs in the system as a whole, and reversing the parameter change will not generally return to the system to its previous state (Fig. 1A).
We perform an invasion analysis for a wide range of fire regimes. We examine normalized mean intensities in the (0.01, 0.99) interval, ranging from intensities that cause only seedling mortality to those that cause complete mortality of mature forest. For each intensity, we determine PFT stability for fire frequencies F ∈ (0.01, 0.35), so mean return times range from 2.86 to 100 yr. In the Klamath, fire return intervals ranged from 4-87 yr, with a median return interval of 12-19 yr for the time period 1626(Taylor and Skinner 1998. Temporal autocorrelation of disturbance can influence community stability (Garrison et al. 2012), but at present our analysis is restricted to a simple Bernoulli distribution for fire occurrence, as described below.

State variables and dynamics
At time t, N jt is a state variable representing the number of seeds and densities of each size class k.
where S jt , J jt , and M jt are the variables for seeds, juvenile, and mature trees of the jth PFT at time t. The juvenile class includes both seedlings and saplings.
The time evolution of populations is described ∧, (2) We define ∧ such that N j,t+1 =∧ jt N jt . The seedbank evolves by seed removal due to mortality δ and germination G, and seed addition from fecundity f.
indicating recruitment from the seedbank, loss to mortality and size transition rate T j .
Competition follows a lottery model (Chesson and Warner 1981): indicates that recruits of PFT j will fill a proportion of the free space that is equal to its proportion of available seeds. 6 For mature trees, indicating recruitment from juveniles and mortality. Because each size class can only affect neighboring size classes, λ 12 =λ 23 =λ 31 =0. ∧ jt is the time evolution operator that depends on PFT and disturbance parameters, the initial conditions, the disturbance process, and the state of the system: Thus, given an initial condition and a realization of the fire process, the density of each PFT and stage for all future time t is N N jt jt jt A full list of symbols and definitions is in Supplementary Material Appendix 2 Table A2.
The spatial extent of our model is that at which the effect of outside seed immigration is minimal (Loreau andMouquet 1999, Chesson 2000a). Because our model does not explicitly include dispersal or spatial heterogeneity, our results are most realistic when applied to patches that have low levels of spatial heterogeneity but are large enough to not gain or lose many seeds through dispersal. For the purposes of this work, we consider the model to represent a patch of 2-5 km 2 . When patch size exceeds about 25 km 2 , the model results may become unrealistic due to larger amounts of internal heterogeneity that are not considered in the model. When patch size is less than about 0.1 km 2 , model results may not be realistic, due to seed dispersal. If one is applying the model to different systems and PFTs, the spatial extent should be determined on the basis of dispersal distances for the PFT in question. Our model is space-filling, meaning that aboveground densities of all PFT stages sum to 1 at the end of every time step. This common assumption reflects situations where most habitat is quickly colonized.

Fire, mortality, environment
We model fire occurrence as a Bernoulli process with mean F, and we consider intensity as a random variable independent of frequency, which allows us to sample a broad range of frequency and intensity combinations and examine whether similar sites can have different stable outcomes. Fire regimes vary depending on factors such as aspect or elevation. For example, fire mortality on ridgetops in the Klamath tends to be higher than at lower elevations on north-facing slopes (Taylor and Skinner 1998), and these differences can be represented in our model by changing the intensity index.
The intensity process (I t ) is zero when no fire occurs. Otherwise we draw intensity from an independent, identically distributed lognormal process with mean μ I and variance σ I . We use a lognormal distribution for our intensity index because this distribution generally fits empirical distributions of fire size (Hantson et al. 2016) for a variety of systems, and in the conifer forests of the northwestern United States, fire size is correlated with severity (Cansler and McKenzie 2014). We present estimates of sensitivity to variation in intensity in Supplementary material Appendix 1 Fig. A6.
We compute mortality in two passes: one for a primary estimate of the effect of fire, and another that assesses the effects of modifiers that are based on PFT densities. The fire response is computed by a sigmoid function of the intensity index, where m jkl , m jks are the location and steepness parameters of the sigmoid, respectively. Two interactions between mortality and density are possible. First, juvenile and mature size classes can act as ladder fuels that increase the fire-induced mortality of mature conifers, and the contribution to this effect is weighted. The ladder weight for PFT j, class k is W Ljk , and the ladder density where L l and L s are location and steepness parameters for the ladder effect. Thus, mature conifer mortality (δ 23t ) can increase but not decrease due to density effects, and the magnitude of this effect decreases to zero as L d →0.
The second density effect is shading of mature shrubs by conifers. Essentially, when the density of mature shrubs is low and conifer densities are high, the mortality of mature shrubs increases. This modeled effect uses a shade-weight parameter for each PFT class, and the shade density is where S l and S S are the location and steepness parameters of the shading effect.
When no fire occurs, the base mortality represents mortality from all sources other than fire, including disease and drought.

Formal stability
We define as the time-evolution matrix for PFT i in the zero-density limit in all stages (i = SH for shrub PFT, i = CON for conifer PFT). The long-term, low-density growth rate quantifies the ability of a population to grow or decline when perturbed to low density, and is where 'Eigs' is the set of eigenvalues. If λ i > 1 , then the density of that PFT will increase, considering both the disturbance regime and competition from the resident. The magnitude of λ i for each PFT allows us to classify four outcomes. If λ CON > 1 , λ SH < 1 , then conifers will eventually exclude shrubs from the system. If λ CON < 1 , λ SH > 1 , conifers are excluded. If λ CON > 1 and λ SH > 1 , then both PFTs stably coexist. If λ CON < 1 , λ SH < 1 , neither PFT can increase from low density in the presence of the other. Eventually, fluctuations will lead the density of one PFT to decrease to the point that it no longer can increase on average. In this case, the system has alternative stable equilibria, and identity of the excluded PFT depends on initial conditions and stochasticity.

Results
We computed the long-term, low-density growth rate for each PFT across a range of fire processes (Fig. 3A). This generated two surfaces of PFT population growth rates as a function of fire intensity and frequency, and we used these to determine stability for each disturbance regime (Fig. 3B). The base parameterization (Table 1) provided stability results that we used as a baseline for analyzing how thresholds and regions of alternative stable equilibria changed in response to changes in disturbance or PFT parameters. Our base parameterization had alternative stable equilibria in 21% of the fire regimes analyzed (Fig. 3B), and alternative stable equilibria occurred with a mean fire return time as long as 25 yr (e.g. I=0.6, F=0.04), within the range of empirical fire data for the Klamath region.

Factors associated with alternative stable equilibria
Two factors affected which PFT was excluded if the system has alternative stable equilibria (Fig. 4). The first was the set of initial densities for each PFT. When different initial conditions were subjected to the same fire process, different outcomes arose. Relatively small differences in the initial density of PFTs determined which PFT eventually dominates (Fig. 4A).
Second, the same initial conditions led to a different PFT being excluded, depending on the sample of the fire process. For instance, given a parameter set that could allow shrub exclusion, a high-intensity fire early in the simulation could lead to shrub dominance (Fig. 4B). Within the alternative stable equilibria region, points closer to the region of shrub Figure 3. Outcome of dynamics assessed by stability analysis, computed across fire frequency (1 yr -1 ) and intensity index. (A) Long term, low density growth rate for each PFT is represented by a surface, and a surface is included at λ = 1. When exactly one PFT has λ > 1, it excludes the other, regardless of initial conditions or random fluctuation. When both PFTs have λ < 1 alternative stable equilibria result, and both initial condition as well as random fluctuation control the eventual equilibrium. (B) Stability and competitive outcome derived from growth surface, color coded by outcome: yellow is alternative stable equilibria, blue is forest persistence, green is shrub dominance. (C) Same as center, with frequency sampled on a log scale to better visualize outcomes at low fire frequencies.
As mean time to maturity for conifers increased, their likelihood of persisting decreased (Fig. 5A-C). In our base parameterization, conifers reached maturity (and correspondingly increased in fire resistance) after 40 yr (Fig. 5A). The disturbance regime (I = 0.5, F = 0.084) generated alternative stable equilibria, and under many initial conditions (e.g. mature conifer fill 65% of habitat), conifers excluded shrubs. When time to maturity was increased to 67 yr (Fig. 5B), then conifer forest did not persist even if its initial density was high, because shrub dominance became the only stable equilibrium. At even longer times to maturity (Fig. 5C), most fire regimes led shrub dominance. Effects were similar as conifer fecundity decreased (Supplementary material Appendix 1 Fig. A2). Changing variation in fire intensity had little effect on areas of conifer dominance, but very low variation increased the size of the region of alternative stable equilibria (Supplementary material Appendix 1 Fig. A6).
Increases in the non-fire mortality rate of shrubs led to forest persistence over a wider range of fire regimes. Such increases also reduced the range of fire regimes under which shrubland was stable, resulting in an upward shift of the alternative stable equilibria region in (I,F) space, with its area remaining relatively unchanged (Fig. 5E, F). We also determined which parameters had the greatest effects on the existence, size, and location of regions of alternative stable equilibria. Alternative stable equilibria were generated by negative interactions between the PFTs (which stabilized the extinction equilbria) acting against the positive effects of growth and recruitment (which stabilized the dominance equilibria) (Fig. 5G-I). When the weighting of shading effects and ladder fuels was decreased to near zero, the alternative stable equilibria region was essentially eliminated (Fig. 5I). Reducing the magnitude of these interactions also decreased the parameter area in which forest can dominate. This suggests that shrub can outcompete conifer over a broad range of fire regimes without any ladder fuel effect, but the shading effect is necessary for conifer to outcompete shrub under many fire regimes.

Discussion
Using a stochastic, stage-based model of community dynamics and fire disturbance, we demonstrated that portions of the Klamath landscape that are currently forested could alternatively be stable as shrubland, and climate change may alter some areas so that only shrubland is stable. Our methods reveal a route to alternative stable equilibria that depends on stage-specific responses to disturbance, allowing for resistant and resilient plant strategies (Miller and Chesson 2009) to be expressed as part of a PFT's life cycle.
Much debate about altered fire regimes centers around the risk that increased fuel buildup since the suppression of frequent fire could cause extensive forest loss through high-severity fire (Stephens et al. 2012). However, a different pattern may characterize systems having alternative stable equilibria. In these systems, fire-free intervals that are sufficiently long for the development of forest could produce a positive feedback, whereby the higher fire resistance of mature trees or a cool, moist microclimate beneath the dense forest canopy can reduce either the probability of fire occurrence or fire severity (Cochrane et al. 1999, Odion et al. 2004, Ray et al. 2005, Hoffmann et al. 2012, Paritsis et al. 2015, Tepley et al. 2016. If these forests are eventually lost to fire, a different positive feedback could sustain vegetation of lower height that develops after fire if it creates a cycle of burning in which fire-free intervals are too short for forest recovery (Odion et al. 2010, Thompson et al. 2011b, Enright et al. 2015, Kitzberger et al. 2016. Empirical studies have suggested that conifer and shrub communities of the Klamath region may constitute alternative stable equilibria (Odion et al. 2010, Thompson et al. 2011b, and our demonstration provides a theoretical basis for alternative stable equilibria in the system. We also illustrated how factors such as maturation time, fecundity, and aspects of fire can independently destabilize forests, potentially causing a shift in the dominant vegetation type (Fig. 3-5). Our model can be used to examine how the responses of these factors to different climate change scenarios may affect forest persistence.
Sites in the Klamath region have a wide variety of fire regimes and climatic conditions, and our model shows that a currently forested site where both community types are stable can rapidly shift to persistent shrubland after small changes to parameters. For example, growth rates of the dominant conifer (P. menziesii) decrease as the climatic water deficit increases (Restaino et al. 2016), which could slow the rate of forest recovery under anticipated climate change. We also found that stabilizing forest by decreasing fire frequency becomes less effective as conifer time to maturity increases. Thus, if a lengthening of the maturation time coincides with increases in fire frequency or intensity, forests already near critical thresholds are likely to convert to shrubland. Loss of forested area would have many effects, including a large decrease in total carbon stocks and a positive feedback to climate change (Birdsey 1992, Dixon et al. 1994, Campbell et al. 2007. We identified factors that stabilize forests, including increased rates of non-fire shrub mortality (Fig. 5) and increasing conifer fecundity (Supplementary material Appendix 1 Fig. A2). We did not test management scenarios, but some of the effects of management techniques can be incorporated into our model, such as prescribed low intensity burning, shrub removal and conifer planting (Bohlman et al. 2016). The Klamath region generally is predicted to become drier and hotter, and hence increases in fire frequency and intensity are more likely than decreases in many locations. Because Klamath conifer forests store large amounts of carbon and have high conifer diversity, even forest loss at small spatial extents could have considerable ecological consequences.
We focused primarily on alternative stable equilibria, but our model can produce coexistence given large enough environmental variation in fecundity and germination (Supplementary material Appendix 1 Fig. A4, A5). Germination variation has a stabilizing effect on coexistence, which is consistent with previous findings that fluctuation in fecundity and germination can generate distinct mechanisms of coexistence (Miller et al. 2012b). Although this is described by prior theory (Chesson 2000a), our model does not have explicit spatial dynamics, and thus does not include seed dispersal from adjacent stands, which may project an unrealistically low likelihood of stable coexistence. We conjecture that forces stabilizing coexistence are more likely to arise through the spatial interactions at larger extents (Chesson 2000b, Melbourne and Chesson 2006, Staal et al. 2016. Even if many local Klamath sites currently can support either stable conifer forest or persistent shrubland vegetation, regional coexistence still can be stabilized by spatial and temporal environmental heterogeneity. Such interactions will likely reduce the size of regions of alternative stable equilibria as spatial extent increases (Shurin et al. 2004, Staal et al. 2016. The spatial structure of fires will likely affect regional stability of cover types, because spatial autocorrelation of fires differentially affects plants with different dispersal distances (Liao et al. 2016). Dispersal distances also interact with fire size, and these effects are not included in our spatially implicit model. Increase in fire size may have a disproportionate effect on stability of conifer forests, given their short-lived seedbanks and dependence on live seed sources to recruit seedlings (Tepley et al. 2017). We did not directly consider the role of fire duration. Long-burning fires have greater likelihood of being subject to severe fire weather and of forming large fires that deplete conifer seeds. Thus, climate change driven increases in fire duration (Dale et al. 2001) would suppress forest regeneration. The terrain of the Klamath may allow forests to persist in areas with favorable microclimates (Wilkin et al. 2016), promoting regional coexistence. We expect alternative stable equilibria to exist in many small areas, leading to a heterogeneous landscape (Odion et al. 2010). However, our model should not be used to draw inferences regarding regional stability of cover types.
To our knowledge, ours is the first study using invasion analysis to specifically investigate alternative stable equilibria. This analysis of long-term, low-density growth rates offers a rigorous way to estimate effects of plant traits and stochastic disturbance on stability. We believe this model and method can be parameterized to provide reasonable estimates of community stability for a given system. It may also be applied to analyzing whether different life history traits and disturbance regimes interact to change stage distributions within a single PFT category. We see our work as a complement to more quantitative predictions and to a portfolio of methods used to ascertain what features of a system determine stability and alternative stable equilibria of plant communities (Serra-Diaz et al. 2018).
Our results demonstrated how alterations of fire regime and PFT parameters can change the location and shape of critical thresholds, and our analysis demonstrated a novel basis for alternative stable equilibria in fire-prone forestshrub communities, via negative interactions between the two community types and differential responses to fire. Focusing on different disturbance regimes can yield insights into species distributions across climate gradients. For example, many species found in the Klamath occur throughout the western United States. Fire regimes with lower frequency and lower mean intensity could be parameterized to examine the stability of forests in the Pacific Northwest, and regimes with higher frequency and consistently high intensity could be parameterized to study chaparral in California. The PFT traits in our model can be parametrized to represent many other PFTs and communities, such as the fynbos of South Africa or the kwongan of Australia (Miller and Chesson 2009). Similarly, by manipulating ladder fuel and shading parameters, grass and tree PFTs could be modeled in a savanna. We hope that our methods can be applied to more plant communities where PFTs and life stages have different responses to disturbance.