Spatial modeling of forest-savanna bistability: impacts of fire dynamics and timescale separation

TL;DR

This paper explores how fire dynamics and vegetation regrowth influence the coexistence of forests and savannas in tropical ecosystems.

Contribution

The study introduces a spatial stochastic model that explicitly integrates fire dynamics and vegetation regrowth to explain forest-savanna bistability.

Findings

Bistability arises from periodic fires maintaining low forest cover and dense forests inhibiting fire spread.

Deterministic mean-field models predict bistability but differ in transient dynamics from spatial models.

Timescale separation between fire and vegetation processes is critical for ecosystem resilience.

Abstract

Forest-savanna bistability – the hypothesis that forests and savannas exist as alternative stable states in the tropics – and its implications are key challenges for mathematical modelers and ecologists in the context of ongoing climate change. To generate new insights into this problem, we present a spatial Markov jump process model of savanna forest fires that integrates key ecological processes, including seed dispersal, fire spread, and non-linear vegetation flammability. In contrast to many models of forest-savanna bistability, we explicitly model both fire dynamics and vegetation regrowth in a mathematically tractable framework. This approach bridges the gap between slow-timescale vegetation models and highly resolved fire dynamics, shedding light on the influence of short-term and transient processes on vegetation cover. In our spatial stochastic model, bistability arises from periodic fires that maintain low forest cover, whereas dense forest areas inhibit fire spread and preserve high tree density. The deterministic mean-field approximation of the model similarly predicts bistability, but deviates quantitatively from the fully spatial model, especially in terms of its transient dynamics. These results also underscore the critical role of timescale separation between fire and vegetation processes in shaping ecosystem structure and resilience.

Figures

Simulations of the spatial-stochastic GBA model. A. Example of the system switching between large fires and mean-field-like behavior (i.e. fluctuations about an ash-dominated state). The y-axis is the proportion of space covered by each state and time *t* in years. B. Snapshot of the spatial system in the stochastic analogue of the ash-dominated mean-field GBA steady state at time $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \be

Two different possible outcomes for forest (dark green line) in the spatial-stochastic FGBA model; the proportion of space covered by each state is on the y-axis with time in years. Both simulations were run with the same parameter values and initial conditions. Forest was randomly distributed and occupied $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$40\%$$\end{document}$ of sites. The remaining sites were all

A simulation of the spatial stochastic FGBA model demonstrating the ability of forest cover to limit the extent of fire spreading and specifically to stop smaller fires from growing into larger fires on the scale of the entire spatial domain (compare to Fig. (9)CD). The color legend is shown in Fig. (9) with time *t* in years. Parameter values: $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50000$$\end

A plot of the computed branch points governing the stable to unstable transition for the GBA steady state for system (4) alongside the analytically predicted branch points in Eq. (13). Branch points are plotted as a function of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document}$ , the rate of fire spread over forest, and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym}

One-parameter bifurcation diagrams for the mean-field FGBA system (3) with partially timescale-separated parameter values with the forest proportion of the solution of the y-axis. A. Varying the fire spread rate over grass ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document}$ ). B. Varying the fire spread rate over forest ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasys

Equilibrium land state proportions (i.e., steady-state values) of the mean-field FGBA system (3) reached starting from a low forest state ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0.01$$\end{document}$ initial condition) for partially timescale separated parameter values varying one system parameter at a time. A. Varying the fire spread rate over grass ( $\documentclass[12pt]{minimal} \usepackage{amsmath

One-parameter bifurcation diagrams for the mean-field FGBA system (3) for fully timescale-separated parameter values with the forest proportion, *F*, on the y-axis and the parameter being varied on the x-axis. A. Varying the fire spreading rate in forest ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document}$ ). B. Varying the grass regrowth rate from ash ( $\documentclass[12pt]{minimal} \usepa

State transition diagram of the spatial FGBA model. Transition arrows are labeled with the relevant parameters and/or flammability functions. Forest and grass/fire timescale transitions are shown in green and orange, respectively

FGBA model state transition diagram in the mean-field limit. Transition arrows are labeled by the transition rates. Forest and fire timescale transitions are shown in green and orange, respectively

Plots of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G$$\end{document}$ (denoted by $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepa

A plot of the computed branch points governing the stable to unstable transition for the GBA steady state for system (4) alongside the analytically predicted branch points in Eq. (13). Branch points are plotted as a function of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document}$ , the rate of fire spread over forest, and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym}

One-parameter bifurcation diagrams for the mean-field FGBA system (3) with partially timescale-separated parameter values with the forest proportion of the solution of the y-axis. A. Varying the fire spread rate over grass ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document}$ ). B. Varying the fire spread rate over forest ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasys

Equilibrium land state proportions (i.e., steady-state values) of the mean-field FGBA system (3) reached starting from a low forest state ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0.01$$\end{document}$ initial condition) for partially timescale separated parameter values varying one system parameter at a time. A. Varying the fire spread rate over grass ( $\documentclass[12pt]{minimal} \usepackage{amsmath

One-parameter bifurcation diagrams for the mean-field FGBA system (3) for fully timescale-separated parameter values with the forest proportion, *F*, on the y-axis and the parameter being varied on the x-axis. A. Varying the fire spreading rate in forest ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document}$ ). B. Varying the grass regrowth rate from ash ( $\documentclass[12pt]{minimal} \usepa

Simulations of the spatial-stochastic GBA model. A. System dynamics on the fire timescale (with time *t* in years), with the proportions of space covered by each state on the y-axis. Note that a slight increase in burning area (fire ignition) is followed by a large increase in ash and a large decrease in grass. The grass cover then gradually recovers from the ash state. B. System dynamics on a grass timescale showing repeated large fire events followed by grass recovery. C. Montage of fire ignition and spreading in grassland on a fire timescale as shown in A. D. Montage of grass regrowth follo

Full text

Introduction

Satellite data of tree cover in Sub-Saharan Africa combined with mean annual rainfall (MAR) records have shown that tree cover exhibits a distinctly bimodal distribution in regions with intermediate MAR (1000 to 2000 mm/yr). Savannas (regions with sparse open-canopy tree cover less than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 40 \%$$\end{document} per unit area) and forests (regions with high closed-canopy tree cover above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 70 \%$$\end{document} per unit area) are frequently observed in the tropics, while intermediate tree cover (40 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$70 \%$$\end{document} ) is rarely observed (Van Nes et al. 2018). The notion of forests and savannas existing as alternative stable states has attracted a significant amount of attention and concern since it suggests that perturbations caused by climate change, drought, or human activity can cause large-scale conversion of forests to savanna, which may be extremely difficult to reverse due to hysteresis effects (Cochrane et al. 1999).

The preservation of forests, in turn, is a crucial concern in mitigating climate change, as forests have significantly higher carbon storage capacities than savannas. Forests also support local human populations by providing building materials. In addition to forests, the maintenance of savannas is also critical to biodiversity and conservation efforts due to the high species richness of ecosystems located at the savanna-forest ecotone (Behling et al. 2007; Myers et al. 2000; Overbeck et al. 2007). Therefore, a model that accurately predicts the responses of forest-savanna mosaics to various environmental disturbances would usefully inform various climate- and biodiversity-related interventions, such as selecting and maintaining protected areas and regulations of human ecological disturbances. As such, a maximally valuable model must encompass short-term disturbances, such as fires, and long-term effects, such as forest encroachment or loss. The model ideally must also be spatially explicit since the fundamental driving processes of the system (e.g., seed dispersal and fire spread) have a highly spatially inhomogeneous structure.

A widely accepted explanation for the observed bimodal distribution of forests and savannas is that regular disturbance by fires can maintain open-canopy savannas. In contrast, denser closed-canopy forests sufficiently suppress the flammable grassy layer, which limits fire spread and lowers the probability of large fires, stabilizing the forest state. Mechanistically, fires typically spread readily within savannas fueled by highly flammable grasses but do not significantly penetrate forest boundaries (Biddulph and Kellman 1998; Hennenberg et al. 2008; Kellman and Meave 1997; Bond 2008; Hoffmann et al. 2009; Swaine et al. 1992). Thus, regions with sufficiently high grass cover may burn frequently enough to maintain savanna conditions despite climatological conditions supporting forest growth (Brando et al. 2012, 2014; Cochrane and Laurance 2002). The original Staver-Levin model (Staver et al. 2011; Staver and Levin 2012) and related models (Accatino et al. 2010; Tilman 1994), which incorporate fire-vegetation interactions in a spatially implicit mean-field approximation, have been widely used to explain and study forest-savanna bistability (Staal et al. 2018; Nes et al. 2014; Yatat et al. 2018; Wuyts et al. 2017). We will use an extension of the Staver-Levin modeling framework as the basis for our study, but our conclusions apply more broadly across this class of spatially implicit forest-savanna models.

Many forest fire models, such as the Staver-Levin model, use mean-field approximations to remain interpretable and analytically tractable, while only indirectly incorporating spatial processes. For example, Klimasara and Tyran-Kamińska (2023); Magnani et al. (2023); Hoyer-Leitzel and Iams (2021) are recently proposed models that use a mean-field (i.e., spatially implicit) approach but add discontinuous, stochastic changes in vegetation proportions to implicitly account for fire spreading events. However, explicitly spatial models have some advantages over non-spatial models, which implicitly incorporate spatial effects (DeAngelis and Yurek 2017). For instance, spatial models have greater practical applicability and more accurately replicate the behavior of systems in which local mechanisms amplified by positive feedback, such as fire, can cause large-scale effects. Thus, to construct an accurate and practically applicable model, it is necessary to include sufficient detail to model both local conditions and the dynamics that they generate. Spatial structure can also profoundly affect the behavior of a system. For instance, Durrett and Ma (2018) found that in a spatial model of grass-forest competition, both states can co-exist only when transition rates are spatially heterogeneous.

Various other spatial models of fires in savanna-forest systems have already been proposed and studied, but have shortcomings. For example, some forest fire models are parametrized based on empirical data (Lasslop et al. 2016; Archibald et al. 2012; Aleman and Staver 2018; Rammer and Seidl 2019), but such models cannot transparently demonstrate the system’s underlying dynamics nor be applied to novel situations. Other studies based on field sampling of species diversity, fire frequency, and/or soil quality have similar drawbacks (Beckett et al. 2022; Sagang et al. 2024). Some mean-field models have been extended to incorporate spatial processes (seed dispersal, fire) by adding diffusion terms to the nonspatial dynamics (Zelnik et al. 2024; Goel et al. 2020), but these models do not explicitly model fire dynamics. Explicit fire dynamics are essential for examining the behavior of grass-forest systems on shorter timescales and for studying the effects of seasonally dependent flammability on forest dynamics (Scholes and Archer 1997). Nonetheless, these studies have already produced interesting results, such as the dependence of forest front travel speed on the curvature of the front and the need for a critical patch size to ensure growth, adding a hysteresis-like quality to the system’s behavior (Durrett and Ma 2018; Goel et al. 2020).

Most spatially explicit forest fire models employ a square grid and use nearest neighbor-dependent transition probabilities to model spatial processes such as fire spread and seed dispersal (Wuyts and Sieber 2023; Staal et al. 2018; Van Nes et al. 2018; Hochberg et al. 1994; Accatino et al. 2016; Li et al. 2019; Abades et al. 2014; Fair et al. 2020). These models have already led to notable findings. For example, Wuyts and Sieber (2023) explicitly demonstrated that fire-vegetation feedback can maintain bistability in tree cover states, Fair et al. (2020) found that forest-savanna mosaics occur under limited conditions, suggesting that such mosaics are vulnerable to loss under climate change, and Abades et al. (2014) found a second order phase transition in spanning cluster formation probability. However, these models have some intrinsic limitations. For instance, they cannot incorporate inhomogeneities in vegetation density and long-distance or anisotropic spreading of fire and seeds. Additionally, the ability of these models to explain observations on short time horizons may be limited by the relatively short timescale of fire dynamics compared to the timestep increment (typically 1 year).

The paper is organized as follows: Section 2 introduces the spatial stochastic modeling framework that forms the basis for this study, including a discussion on plausible parameter regimes and rigorous links to continuum models. Section 3 provides a bifurcation analysis of a nonspatial version of the model for later comparison with the spatial version and other similar nonspatial forest-savanna models in the extant literature. Section 4 studies the dynamics of the full spatial stochastic model and contrasts them with the dynamics of the nonspatial model from the previous section, emphasizing the key role of transient dynamics on short timescales. Finally, Section 5 concludes with a synthesis of the main findings and highlights some future directions and open problems in the mathematical modeling of forest-savanna bistability.

A Spatial Stochastic Forest Fire Model

Ecologically, tropical forests differ markedly from coniferous forests in several key ways. First, savanna trees are not readily killed by fire. Typically, they are only top-killed (i.e., only aerial biomass is burned) and can readily resprout (Hoffmann et al. 2003, 2009; Balch et al. 2011) or have thickened bark to prevent stem death (Hoffmann et al. 2012). Second, fire does not propagate readily through tropical forests (Hennenberg et al. 2006; Archibald et al. 2009; Pueyo et al. 2010) but does propagate rapidly through savanna grassland (Wragg et al. 2018). Mechanistically, this is due to forest understory shade excluding flammable C4 grasses (Hoffmann et al. 2012) in addition to reduced wind speeds and increased moisture in the forest microclimate (Hoffmann et al. 2012; Hennenberg et al. 2008). For brevity, we will henceforth refer to tropical forests as “forest” and savanna regions with lower tree density as “grassland”.

Mathematical Description of the Model

We use a spatially extended Markov jump process to model the dynamics of fires in forest-grass systems using the mathematical framework introduced in Patterson et al. (2020). The set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega \subset \mathbb {R}^2$$\end{document} denotes the spatial domain and the random variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{r_i \in \Omega :\, i = 1,\dots , N \}$$\end{document} are sites (locations) where vegetation is present or has the potential to grow. Each site i, with respective location \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_i$$\end{document} , undergoes a Markov jump process with possible states F (forest), G (grass), B (burning), and A (ash). We refer to models with this state space as the “FGBA model” for brevity. The “ash” state refers to areas of land that have recently burned and thus do not contain flammable cover or vegetation that can immediately burn (Hébert-Dufresne et al. 2018; Wuyts et al. 2017). We can thus consider ash to be a transitory or refractory state, and it is often omitted from models that only consider vegetation on longer timescales (years, decades). Still, it can play an important role in ecosystem dynamics on shorter timescales, which is the focus of the present work.

The transition rates between the FGBA model states depend on the spatial locations and states of all the other sites at locations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{r_j \in \Omega \,: \, j = 1,\dots ,N, j \ne i\}$$\end{document} . This allows us to model both spatial spreading processes (e.g. forest and fire spreading) as well as non-spatial spontaneous transitions (e.g. fires ignited due to lightning strikes or human activity (Archibald et al. 2012)). The spatial locations of the vegetation can be drawn from any sufficiently regular distribution on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} , but, in this paper, we draw from a uniform distribution on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} to reflect a homogeneous domain scenario. In the limit as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \rightarrow \infty $$\end{document} , we completely fill the spatial domain with possible sites for vegetation. We can then obtain an approximate continuum model that can be studied as either a spatial mean-field model (via a system of integro-differential equations) or as a non-spatial mean-field model that presents as a system of ordinary differential equations (see Section 2.2 below for details).

We make the following assumptions about the possible transitions in our model:

1. Forest trees can expand into nearby grass and ash due to local spreading of seeds,
2. burning sites can expand into nearby forest and grass due to local fire spread,
3. Burning sites are spontaneously quenched into ash at a fixed rate,
4. Grass can regrow spontaneously from ash at a fixed rate due to homogeneous dispersal of grass seeds,
5. Forest can spontaneously transition to grass at a fixed rate due to non-fire related mortality. Mathematically, we represent these assumptions by allowing each site to transition between the states F, G, B, A at exponentially distributed times with rates given by the sum of spatial spreading and spontaneous processes (see Table 1). The parameters in the transition rates are defined as follows:

• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _G, \varphi _A$$\end{document} are constants controlling the rate of forest seeding into grass and into ash, respectively.
• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F, \beta _G$$\end{document} are constants controlling the rate of fire spread within forest and within grass, respectively.
• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_F, W_B: \Omega \times \Omega \rightarrow \mathbb {R}^+$$\end{document} are forest spread and burning spread kernels, respectively, which control the extent of spatial spreading as a function of distance between the interacting sites. We assume that the kernels are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {C}^{\infty }$$\end{document} . Table 1. Transition rates between states in the spatial-stochastic FGBA model and their corresponding ecological processesTransitionTransition rate at site iEcological process \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G \rightarrow F$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \frac{ \varphi _G}{N} \sum _{j=1}^N W_F(r_i,\,r_j) 1\!\!1_{\{X^j(t)=F\}}$$\end{document} forest spreading into grass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A \rightarrow F$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\varphi _A}{N} \sum _{j=1}^NW_F(r_i,\,r_j) 1\!\!1_{\{X^j(t)=F\}}$$\end{document} forest spreading into ash \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F \rightarrow B$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\beta _F}{N}\sum _{j=1}^N W_B(r_i,\,r_j)1\!\!1_{\{X^j(t)=B\}}$$\end{document} fire spreading into forest \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G \rightarrow B$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\beta _G}{N}\sum _{j=1}^NW_B(r_i,\,r_j)1\!\!1_{\{X^j(t)=B\}}$$\end{document} fire spreading into grass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F \rightarrow G$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} non-fire forest mortality \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B \rightarrow A$$\end{document} qfire quenching \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A \rightarrow G$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} grass regrowth from ash

For the neighbor (local) transition rate parameters for fire and forest invasion ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _G$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _A$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} ), the subscript is the first English letter of the state before the relevant transition and the main letter is the Greek letter of the state after the relevant transition. We use Greek letters without subscripts for spontaneous transition rates ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu , q, \gamma $$\end{document} ). We normalize the neighbor spread rates by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\text {area}(\Omega )}{N}$$\end{document} to ensure that the spreading rates remain bounded in the limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \rightarrow \infty $$\end{document} .

In addition to the diffusive mode of fire spread described by the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_B$$\end{document} kernel and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G, \beta _F$$\end{document} parameters, our modeling framework also includes a “cascade” mode of vegetation burning motivated by percolation theory. A key finding of detailed percolation-based fire models is that the probability of ignition for grass sites increases from a low baseline level to near one once the fraction of grass sites in the landscape exceeds a critical percolation threshold (with the exact value depending on the specific model) (Abades et al. 2014; Schertzer et al. 2015). This mechanism is also well supported by empirical data for forest-savanna ecosystems (Archibald et al. 2009; Pueyo et al. 2010; Van Nes et al. 2018). We incorporate this effect of a critical threshold for fire spread as a function of flammable cover by adding the terms

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Phi _G\bigg (\frac{1}{N}\sum _{i=1}^N W_G(r_i,\,r_j)1\!\!1_{\{X^j(t)=G\}}\bigg ) \quad \text { and } \quad \Phi _F\bigg (\frac{1}{N}\sum _{i=1}^N W_G(r_i,\,r_j)1\!\!1_{\{X^j(t)=G\}}\bigg ) \end{aligned}$$\end{document}

to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G \rightarrow B$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F \rightarrow B$$\end{document} transition rates of site i at time t. Here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_G: \Omega \times \Omega \rightarrow \mathbb {R}^+$$\end{document} is another smooth spreading kernel, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G, \Phi _F:\mathbb {R}^+ \rightarrow \mathbb {R}^+$$\end{document} are smooth sigmoidal functions which output flammability as a function of the local grass cover and the kernel \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_G$$\end{document} . More explicitly, we assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G(\cdot )$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _F(\cdot )$$\end{document} will take the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Phi _G(x) = g_0 + \frac{g_1 - g_0}{1 + e^{-(x-\theta _G)/s_G}} \quad \text { and } \quad \Phi _F(x) = f_0 + \frac{f_1 - f_0}{1 + e^{-(x-\theta _F)/s_F}} \end{aligned}$$\end{document}

where the parameters are defined as follows:

• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0, g_0$$\end{document} are the baseline spontaneous flammability of forest or grass sites when no other grass sites are present, e.g. ignitions due to lightning strikes or human activity (Archibald et al. 2012).
• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1, g_1$$\end{document} are the total flammability of a grass or forest site, respectively in the limit when the land patch is purely grassland.
• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{F}, \theta _G$$\end{document} are the percolation thresholds in forest and grass, respectively. We will use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta = \theta _F = \theta _G \approx 0.6$$\end{document} as used for percolation in a square lattice (Gebele 1984).
• \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_F, s_G$$\end{document} are non-negative constants controlling the width of forest and grass sigmoids, respectively. The sigmoids are assumed to be smooth approximations of increasing step functions, so we set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_F = s_G = 0.05$$\end{document} . The FGBA model outlined above spans multiple distinct timescales. Fire dynamics occur on a time scale of hours, grass regrowth on a time scale of months, and forest dynamics occur on a time scale of decades. In particular, we will estimate the parameter values from the expected time between events. Since we assume the transitions follow an exponential distribution, the rates are the inverse of the expected time between events. The rate parameters all have the same units of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} , and we assume that the vegetation sites are separated by an average distance of 1 m.

The vegetation sites are chosen randomly within a compact square domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[0, L] \times [0, L] \subset \mathbb {R}^2$$\end{document} for some fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L >0$$\end{document} . We use periodic boundary conditions to reduce boundary effects and to model an infinite domain. The distance between two points \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r, s \in [0, L] \times [0,L]$$\end{document} is then computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} |r - s|^2 = \frac{L}{2\pi } \Big ( \Big [ \text {arg}\Big (e^{i \frac{2 \pi (r_x - s_x)}{L}}\Big ) \Big ]^2 + \Big [ \text {arg}\Big (e^{i \frac{2 \pi (r_y - s_y)}{L}}\Big ) \Big ]^2 \Big ) \end{aligned}$$\end{document}

where the argument function has range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[0, 2\pi )$$\end{document} .

We use Gaussian functions for the spreading kernels \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_G$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_F$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_B$$\end{document} , although in principle, different functions could be used to model alternative spreading mechanisms if desired (such as very long-range seed dispersal modeled by heavy-tailed kernels). A vegetation site located at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_i \in \Omega \subset \mathbb {R}^2$$\end{document} has spreading kernels

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} W_{\mathcal {X}}(r_j, r_i)&= \frac{1}{2 \pi \sigma _{\mathcal {X}}^2}\text {exp}\Big (-\frac{|r_j - r_i|^2}{2\sigma _{\mathcal {X}}^2}\Big ) \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {X} \in \{G,\, F,\,B\}$$\end{document} . To preserve the physical meaning of the three parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\mathcal {X}}$$\end{document} under changes in L and/or N, we choose the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\mathcal {X}}$$\end{document} to be in units of the average spacing between the vegetation sites in the domain which we define as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x \equiv \frac{L}{\sqrt{N}}$$\end{document} . The normalization constant of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_{\mathcal {X}}$$\end{document} is thus chosen such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lim _{N \rightarrow \infty } \frac{1}{N}\sum _{i=1}^NW_{\mathcal {X}}(r_j, r_i)&= \int _{\Omega } W_{\mathcal {X}}(r_j, r_i) \, dr_i = 1, \end{aligned}$$\end{document}

We refer to the spatial Markov jump-process model outlined above as the spatial-stochastic FGBA model to avoid confusion with the mean-field versions (both spatial and non-spatial) that will also feature in our analysis and discussion. The model described in this section is similar to previous models proposed by Hébert-Dufresne et al. (2018) and Wuyts and Sieber (2023), particularly the choice of the state space, but has several distinguishing features. Firstly, our model allows the sites to be distributed arbitrarily in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} , while the related models use a square lattice of sites. Secondly, Hébert-Dufresne et al. (2018) and Wuyts and Sieber (2023) include additional spontaneous transitions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G, A, \rightarrow F$$\end{document} to model the homogeneous long-distance dispersal of forest seeds. Our model more realistically incorporates long-distance forest seed dispersal by allowing an appropriately heavy-tailed kernel for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_F$$\end{document} . Lastly, the related models assume that spreading processes occur strictly via nearest-neighbor spreading across adjacent lattice sites, while our model allows for more general spreading mechanisms via heavy-tailed kernels often used for seed dispersal (Nathan et al. 2012). Moreover, as we show below, our spatial stochastic model has immediate, rigorous connections with both spatial and non-spatial mean-field versions, which greatly assist with analysis and interpretability (see Section 2.2).

Parameter estimates

In this section, we set reasonable ranges and approximate values for all parameters in the FGBA model. We assume that L represents the physical side length of the area of land under study. Then, the average spacing between vegetation sites, assuming a uniform probability distribution, is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x \equiv \frac{L}{\sqrt{N}}$$\end{document} . We now estimate the parameter values by assuming \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x = 1$$\end{document} m to fix the spatial scale of the model and obtain plausible order-of-magnitude parameter estimates relative to this scale. Our results are not intended to be quantitative in nature, but rather we will illustrate dynamics possible over broad ranges of the plausible parameter space. When adjusted for average grid size differences, our parameter estimates are in line with those of related spatial stochastic forest-savanna models that used empirically validated parameters (Hébert-Dufresne et al. 2018; Wuyts and Sieber 2023). Our baseline parameter estimates are as follows:

• A square patch of grassland of side length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x$$\end{document} surrounded by fire will burn after several minutes. Using Table 1 and Eq. (1), the burning rate can be approximated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G \approx 10^5$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} .
• A square patch of forest of side length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x$$\end{document} , which is surrounded by fire, will burn after about an hour. The burning rate can be approximated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F \approx 10^4$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} .
• A burning site is expected to burn for several hours before turning into ash, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \approx 10^3$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document}
• An ash site is expected to regrow grass in several months, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \approx 10^1$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} .
• A square patch of grass or ash of side length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x$$\end{document} surrounded by forest will grow a forest tree after several years to a decade, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _A \approx \varphi _G \approx 10^{0}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-1}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} .
• A tree at a square forest site of side length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x$$\end{document} will spontaneously die from non-fire related causes after about 100 years, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \approx 10^{-2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} .
• A square grass or forest vegetation site of side length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x$$\end{document} completely surrounded by forest will spontaneously catch fire once every 100 years, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0, g_0 \approx 10^{-2}$$\end{document} .
• A grass or forest vegetation site will have increased flammability when completely surrounded by grass compared to a forest. We set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = 10^{-1}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_0 = 10^0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} .
• The standard deviation of a fire site spreading throughout grass or forest vegetation is approximately 5 m, so set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _F \approx 5$$\end{document} m.
• The standard deviation of a forest site spreading throughout grass or ash is approximately 5 m, so set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _B \approx 5$$\end{document} m.
• The flammability of a vegetation site is dependent on grass proportions in a 5 m radius up to one standard deviation, so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _G \approx 5$$\end{document} m. We summarize the parameter values, and a representative fixed value is given below in Table 2. In this table, G(t) represents the fraction of the landscape occupied by grass sites.Table 2. Estimates and representative values of the parameters for the FGBA model in units of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} ParameterValueUnitsEcological Interpretation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^5$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} magnitude of fire spread rate over grass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^4$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} magnitude of fire spread rate over forestq \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^4$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} rate at which fire is quenched \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^2$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} rate at which grass regrows from ash \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _A$$\end{document} 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} magnitude of forest spread rate over ash \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _G$$\end{document} 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} magnitude of forest spread rate over grass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} tree mortality rate due to non-fire-related causes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} spontaneous fire rate for each forest site when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(t) \approx 0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-1}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} spontaneous fire rate for each forest site when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(t) \approx 1$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} spontaneous fire rate for each grass site when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(t) \approx 0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-1}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {yr}^{-1}$$\end{document} spontaneous fire rate for each grass site when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(t) \approx 1$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _F$$\end{document} 5mwidth of forest spreading kernel \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _B$$\end{document} 5mwidth of fire spreading kernel \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _G$$\end{document} 5mwidth of spontaneous fire kernel

There are three levels of time separation in the rate parameters given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G,\beta _F, q > \gamma \gg \varphi _A, \varphi _G, \mu , f_0, g_0$$\end{document} . This concludes our mathematical description of the spatial FGBA model. We summarize the model in the state transition diagram in Fig. (1), showing the relevant parameters and functions governing transitions between states at each vegetation site.Fig. 1. State transition diagram of the spatial FGBA model. Transition arrows are labeled with the relevant parameters and/or flammability functions. Forest and grass/fire timescale transitions are shown in green and orange, respectively

Mean-Field Approximations of the Spatial FGBA Model

Theorem 2.2.1 in Patterson et al. (2020) shows that the spatial-stochastic FGBA model with uniformly randomly distributed sites on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} converges to a spatial McKean-Vlasov process in the limit as N (the number of sites) tends to infinity. We may think of this limiting processes as a type of spatial mean-field limit since it retains information about spatial structure, but the vegetation at each location only interacts with other locations through weighted spatial averages of all other locations. Let X(r, t) indicate the probability that the limiting McKean-Vlasov process at location r and time t is in state X. Then the forward equation governing the dynamics of the probability densities of the limiting process are given by the following system of integro-differential equations (IDEs):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{d}{dt} F(r,t)&= (\varphi _G G(r,t) + \varphi _A A(r,t))\int _{\Omega }W_F(r,r')F(r',t)dr' - \mu F(r,t)\nonumber \\&\quad -\beta _F F(r,t) \int _{\Omega }W_B(r,r')B(r',t)dr' -\Phi _F\bigg [\int _{\Omega }W_G(r,r')G(r',t)dr'\bigg ] F(r,t),\nonumber \\ \frac{d}{dt}G(r,t)&= \gamma A(r,t) - \varphi _G G(r,t) \int _{\Omega }W_F(r,r')F(r',t)dr' +\mu F(r,t) \nonumber \\&\quad -\Phi _G\bigg [\int _{\Omega }W_G(r,r')G(r',t)dr'\bigg ] G(r,t) - \beta _G G(r,t) \int _{\Omega }W_B(r,r')B(r',t)dr',\nonumber \\ \frac{d}{dt}B(r,t)&= ( \Phi _F F(r,t) + \Phi _G G(r,t) )\bigg [\int _{\Omega }W_G(r,r')G(r',t)dr'\bigg ] - qB(r,t) \nonumber \\&\quad + (\beta _G G(r,t) + \beta _F F(r,t)) \int _{\Omega }W_B(r,r')B(r',t)dr',\nonumber \\ \frac{d}{dt} A(r,t)&= q B(r,t) - \gamma A(r,t) - \varphi _A A(r,t) \int _{\Omega }W_F(r,r')F(r',t) dr', \end{aligned}$$\end{document}

for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r \in \Omega $$\end{document} . Moreover, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(r,t)+G(r,t)+B(r,t)+A(r,t) = 1$$\end{document} for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(r,t)\in \Omega \times \mathbb {R}_+$$\end{document} .

To gain a rough, intuitive understanding of the stochastic FGBA model’s dynamics, we will analyze a more analytically tractable nonspatial mean-field model derived from (2) that takes the form of a system of ordinary differential equations. The general structure of the solution space of the spatial model (2) is not analytically tractable. One can investigate the linear stability of spatially homogeneous steady-state solutions, but we expect similar stability criteria to the nonspatial mean-field case for most reasonable choices of the kernel functions, as has been shown for related forest-savanna models posed in this framework (Patterson et al. 2024). Using the nonspatial mean-field approximation results as a baseline will also allow us to isolate the impact of spatial structure in the stochastic FGBA model. In addition to spatial structure, the stochastic FGBA model allows us to retain stochastic, discrete and finite-size effects, which we expect to be especially important for short-time and transient dynamics (Hastings et al. 2018), which we can contrast against the nonspatial mean-field model.

To obtain a nonspatial (spatially implicit) mean-field model from (2), we assume that the types of states are well-mixed throughout the spatial domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} so that interactions depend only on the fractions of land occupied by each state. One way to achieve this is to take the kernels functions to be uniform on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} and then define the land cover proportions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(t):= \int _\Omega F(r',t)dr'/\text {area}(\Omega )$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(t):= \int _\Omega G(r',t)dr'/\text {area}(\Omega )$$\end{document} , and so on. The system (2) then simplifies to the system of ordinary differential equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{d}{dt}F(t) = (\varphi _G G + \varphi _A A) F -\beta _F B F -\Phi _F(G)F- \mu F \\ \frac{d}{dt} G(t) = \gamma A +\mu F- \varphi _GFG -\Phi _G(G) G - \beta _GBG \\ \frac{d}{dt} B(t) =\Phi _G(G) G + \Phi _F(G)F + (\beta _G G + \beta _F F)B -q B\\ \frac{d}{dt} A(t) = q B - \gamma A - \varphi _AF A. \end{array}\right. } \end{aligned}$$\end{document}

The state transition diagram for the spatial FGBA model in the mean field limit is illustrated in Fig. (2). First notice that due to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G \xrightarrow {\Phi (G)G} B \xrightarrow {qB} A \xrightarrow {\gamma A} G$$\end{document} cycle then if any one of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G}, \overline{B}$$\end{document} , or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{A}$$\end{document} is nonzero then all three must be nonzero. Next, since there is an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F \xrightarrow {\mu F}G$$\end{document} transition then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{F} > 0$$\end{document} implies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G} > 0$$\end{document} . Thus, we can classify all steady states as either GBA (where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{F} = 0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G}, \overline{B}, \overline{A} \ne 0$$\end{document} ) or FGBA (where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{F}, \overline{G}, \overline{B}, \overline{A} \ne 0$$\end{document} ).Fig. 2FGBA model state transition diagram in the mean-field limit. Transition arrows are labeled by the transition rates. Forest and fire timescale transitions are shown in green and orange, respectively

Analysis of the nonspatial FGBA model

Before studying the full spatial stochastic FGBA model introduced above, we perform some qualitative analysis of solutions to the nonspatial mean-field system (3).

GBA Steady States

As it is the simplest subsystem of the model, we first consider the case when forest is absent from the landscape. Set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0$$\end{document} and use the relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 1 - G - B$$\end{document} to eliminate A from (3) to obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{d}{dt}G(t) & = \gamma (1-G-B) - \Phi _G(G) G- \beta _GBG \\ \frac{d}{dt} B(t) & = \Phi _G(G)G + \beta _G GB - qB. \end{array}\right. } \end{aligned}$$\end{document}

Establishing the forward invariance of the GBA subspace is straightforward, and so the proof of the following proposition is deferred to Appendix A.

Proposition 1

For any ecologically relevant initial value (i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(0) \ge 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B(0) \ge 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(0) + B(0) \le 1$$\end{document} ), the solution (G(t), B(t)) for equations (4) remains ecologically relevant for all times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \in \mathbb {R}_+$$\end{document} .

Next, we solve for the steady-state solutions of the model by setting both equations in the system (4) equal to 0. Rearranging the first equation in (4), shows that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{B} = \frac{\gamma (1-\overline{G})}{\gamma + \beta _G \overline{G}} - \frac{\Phi _G(\overline{G})\overline{G}}{\gamma + \beta _G \overline{G}}. $$\end{document}

Insert this relation into the second equation in (4) and multiply across by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma + \beta _G \overline{G}$$\end{document} to obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ 0 = \Phi _G(\overline{G})\overline{G} \left( \gamma + \beta _G \overline{G} \right) +\beta _G \gamma (1-\overline{G})\overline{G} - \beta _G \Phi _G(\overline{G})\overline{G}^2 -q \gamma (1-\overline{G}) + q \Phi _G(\overline{G})\overline{G}. $$\end{document}

Rearranging the equation above shows that the steady-state solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G}>0$$\end{document} obeys:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Phi _G(\overline{G})&= (1-\overline{G})\frac{\gamma q}{\gamma + q} \bigg (\frac{1}{\overline{G}}-\frac{\beta _G}{q}\bigg ). \end{aligned}$$\end{document}

Unfortunately, closed-form expressions for the equilibria cannot be obtained except in special cases due to the sigmoidal properties of the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G$$\end{document} . However, we can readily deduce the uniqueness of the steady state. This analysis naturally involves studying properties of the function

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {F}(G) := (1-G)\frac{\gamma q}{\gamma + q} \bigg (\frac{1}{G}-\frac{\beta _G}{q}\bigg ), \quad G \in (0,1), \end{aligned}$$\end{document}

which is the right-hand side of equation (5).

Proposition 2

For any positive parameter values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma , q,$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} there exists a unique GBA steady state.

Proof of Proposition 2

By definition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G}$$\end{document} is a steady state if and only if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}(\overline{G}) = \Phi _G(\overline{G})$$\end{document} . Let us first consider the case where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q < \beta _G$$\end{document} . Then in the interval [0, 1], \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}$$\end{document} has exactly two roots located at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{q}{\beta _G}$$\end{document} and 1. Furthermore, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F} > 0$$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, q/\beta _G)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F} < 0$$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(q/\beta _G,1)$$\end{document} . Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G > 0$$\end{document} on (0, 1] and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 \le \Phi _G(0) \le \infty $$\end{document} then any roots of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F} - \Phi _G$$\end{document} in [0, 1] must occur in the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, q/\beta _G)$$\end{document} . Next note that the derivative of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}$$\end{document} is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {F}'(G) = \frac{\gamma q}{\gamma + q}\bigg (\frac{\beta _G}{q} - \frac{1}{G^2} \bigg ) \end{aligned}$$\end{document}

and has a single root in [0, 1] located at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{q/\beta _G}$$\end{document} . Furthermore, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}' < 0$$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, \sqrt{q/\beta _G})$$\end{document} . Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, q/\beta _G) \subset (0, \sqrt{q/\beta _G})$$\end{document} it follows that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}$$\end{document} is strictly monotonically decreasing on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, q/\beta _G)$$\end{document} . Next since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G$$\end{document} is monotonically increasing on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, q/\beta _G)$$\end{document} then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F} - \Phi _G$$\end{document} is strictly monotonically decreasing on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, q/\beta _G)$$\end{document} . Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F} - \Phi _G$$\end{document} is continuous on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {R}^+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lim _{G \rightarrow \infty }(\mathcal {F} - \Phi _G)(G) = +\infty $$\end{document} while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\mathcal {F}- \Phi _G)(q/\beta _G) = -\Phi _G(q/\beta _G) < 0$$\end{document} then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F} - \Phi _G$$\end{document} must have exactly one root in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0, q/\beta _G)$$\end{document} . It follows that there is a unique steady state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G} \in [0,1]$$\end{document} .

When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \ge \beta _G$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}$$\end{document} has exactly one root at 1 in the interval [0, 1] and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F} > 0$$\end{document} on (0, 1). By the same argument as before, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}' < 0$$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0,\sqrt{q/\beta _G})$$\end{document} , so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}$$\end{document} is strictly monotonically decreasing on (0, 1). The rest of the argument is the same as the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q < \beta _G$$\end{document} case with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q/\beta _G$$\end{document} replaced by 1. Thus, a unique GBA steady state always exists for any choice of system parameters. See Fig. (3) below for a graphical illustration.Fig. 3. Plots of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G$$\end{document} (denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} ) in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q < \beta _G$$\end{document} case (left) and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \ge \beta _G$$\end{document} case (right)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

Graphically, from Fig. (3), we note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G}$$\end{document} increases as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\gamma q}{\gamma + q}$$\end{document} increases. Noting that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial _{\gamma }(\frac{\gamma q}{\gamma + q}) > 0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial _{q}(\frac{\gamma q}{\gamma + q}) > 0$$\end{document} when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma , q > 0$$\end{document} it follows that increasing either \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} or q will increase the proportion of grass in the steady state. This is logical since increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} allows faster regrowth of grass from ash, and increasing q allows faster fire quenching into ash from which grass can regrow. Similarly, increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{q}{\beta _G}$$\end{document} also increases \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G}$$\end{document} . This also follows intuition: the case of increasing q has already been discussed, while decreasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} reduces the rate at which grass burns due to fire spread. We also note that increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_0$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1$$\end{document} reduces \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G}$$\end{document} as expected, since doing so increases the flammability of grass. Lastly, we note that the presence of the GBA steady state depends on the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G$$\end{document} is non-decreasing. Thus, if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G$$\end{document} were constant, the result of the system containing a single unique GBA steady state would still hold.

Stability of the GBA Steady State

In this section, we demonstrate two results regarding the stability of the GBA steady state. The first concerns stability within the GBA subsystem, and the second gives an interpretable condition characterizing when we can expect the GBA subsystem to be stable to forest invasion.

Theorem 3

The GBA steady state is stable to perturbations within the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0$$\end{document} subspace for any parameter values.

Proof

We compute the linearization matrix at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0$$\end{document} of the system (4), which yields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf {J}}(F,\,G,\,B)\vert _{F=0} = \begin{bmatrix} -\gamma -\Phi _G'(G)G-\Phi _G(G)-\beta _G B,~\,\,\,\, & -\gamma - \beta _G G \\ \Phi _G'(G) G + \Phi _G(G) + \beta _G B, & -q + \beta _G G \end{bmatrix}. \end{aligned}$$\end{document}

The fixed point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\overline{G}, \overline{B})$$\end{document} is stable to perturbations within the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0$$\end{document} boundary exactly when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_{F=0}(\overline{G}, \overline{B})$$\end{document} has eigenvalues with strictly negative real parts. By the Routh-Hurwitz stability criteria, this occurs if and only if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\,\textrm{Tr}\,}}({\textbf {J}}_{F=0}(\overline{G}, \overline{B})) < 0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det ({\textbf {J}}_{F=0}(\overline{G}, \overline{B})) >0$$\end{document}  (Edelstein-Keshet 2005). Note that since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G} < \frac{q}{\beta _G}$$\end{document} then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-q + \beta _G \overline{G} < 0$$\end{document} so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\,\textrm{Tr}\,}}(J_{F=0}(\overline{G}, \overline{B})) < 0$$\end{document} is clearly satisfied. The condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det ({\textbf {J}}_{F=0}(\overline{G}, \overline{B})) >0$$\end{document} is equivalent to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Phi '(\overline{G}) > \frac{\gamma q}{(\gamma + q)^2}\bigg (\frac{\beta }{q}-\frac{1}{\overline{G}^2} \bigg ) = \mathcal {F}'(\overline{G}), \end{aligned}$$\end{document}

which holds from graphical inspection or more rigorously by recalling from our earlier argument that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G'(\overline{G}) \ge 0$$\end{document} while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}'(\overline{G}) < 0$$\end{document} . \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

Proposition 4

The GBA steady state is stable to invasions by F exactly when

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta _F \overline{B} + \Phi _F(\overline{G}) + \mu > \varphi _G \overline{G} + \varphi _A \overline{A}. \end{aligned}$$\end{document}

Proof

We compute the linearization matrix for the full FGBA system and evaluate it at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0$$\end{document} . We use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 1 - F - G - B$$\end{document} to obtain a dynamical system in the three variables F, G, and B which gives the matrix

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&{\textbf {J}}(F,\, G,\, B)\vert _{F = 0} \\&= \begin{bmatrix} \varphi _G G + \varphi _A A - \beta _F B -\Phi _F(G) - \mu , & 0, & 0 \\ -\gamma -\varphi _G G +\mu , & -\gamma -\Phi _G'(G)G-\Phi _G(G)-\beta _G B,~ & -\gamma - \beta _G G \\ \Phi _F(G) + \beta _F B, & \Phi _G'(G) G + \Phi _G(G) + \beta _G B, & -q + \beta _G G \end{bmatrix}. \end{aligned}$$\end{document}

The GBA steady state is stable to invasion by F if all the eigenvalues of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {J}}(\overline{F}, \overline{G},\overline{B})\vert _{\overline{F} = 0}$$\end{document} have negative real parts. We have already shown that the eigenvalues of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 \times 2$$\end{document} submatrix in the lower right have strictly negative real parts so the GBA steady state is stable if and only if the entry in the upper left is negative. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

An intuitive interpretation of equation (7) is that the GBA equilibria is resistant to invasion by forest when the rate of tree mortality by fire spread through forest, spontaneous forest fires, and natural mortality (left-hand side) exceeds the rate of forest seeding into the steady state grass and ash land (right-hand side). Since a unique GBA steady state exists for any value of the system parameters, but its stability can vary based on the parameter values, we expect to find transcritical bifurcations within the parameter space when equality occurs in equation  (7). In other words, a stable GBA steady state bifurcates into an unstable no-forest state and a stable forest state.

Some example bifurcation diagrams in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = \varphi _G = \varphi _A$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} using partially timescale-separated parameter values are given in Fig. (6). As expected, when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} is increased past a critical value or when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} are decreased past a critical value, the proportion of forest sites within \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} becomes non-zero. However, the parameter values produce a mostly ash grassland steady state which is not very ecologically realistic. Nonetheless, it will later be shown that the spatial FGBA model has grassy steady states distinct from the ash-dominated, spatially homogeneous GBA steady state.

Estimates of the GBA Steady State

We will now show that using our estimates of the parameter values, the steady state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\overline{G}, \overline{B}, \overline{A})$$\end{document} can be easily estimated. Notice that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}(G)$$\end{document} becomes asymptotic to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\gamma q}{\gamma + q} (\frac{1}{G} - \frac{\beta _G}{q})$$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G \rightarrow 0$$\end{document} and the coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\gamma q}{\gamma + q}$$\end{document} is of order 2. Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G(G)$$\end{document} is of order 0 or lower for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G \in [0,1]$$\end{document} then we can approximate the condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _G(G) = \mathcal {F}(G)$$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 = \mathcal {F}(G)$$\end{document} which gives \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{G} \approx \frac{q}{\beta _G}$$\end{document} . It follows that to increase the portion of land covered in grass, we must either increase q, i.e. increasing the fire quenching rate, or decrease \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} , i.e. decreasing the rate of fire spread through grass, both of which make sense intuitively. By solving System 4 and using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G+B+A=1$$\end{document} we can compute an estimate of the GBA steady state as well as the ratio \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{B}/\overline{A}$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} (\overline{G}, \overline{B}, \overline{A}) \approx \bigg (\frac{q}{\beta _G}, \frac{\gamma (1-\frac{q}{\beta _G})}{q+\gamma }, \frac{q(1 - \frac{q}{\beta _G})}{q+\gamma }\bigg ) \Longleftrightarrow \frac{\overline{B}}{\overline{A}} \approx \frac{\gamma }{q}. \end{aligned}$$\end{document}

Based on the estimate, the ratio between the portions of land in burning and ash states is determined by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} , the rate of grass regrowth from ash, and q, the rate of fire quenching in a manner that agrees with intuition.

Then Eq. (8) can be used to give Eq. (7) purely in terms of the model parameters:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\beta _F \gamma (1 - \frac{q}{\beta _G})}{q+\gamma } + \Phi _F\Big (\frac{q}{\beta _G}\Big ) + \mu > \frac{q\varphi _G}{\beta _G} + \frac{q\varphi _A(1-\frac{q}{\beta _G})}{q+\gamma }. \end{aligned}$$\end{document}

In particular, if we vary each of the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F, \gamma , f_0, \mu ,$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \varphi = \varphi _A = \varphi _G$$\end{document} one at a time and fix all other parameters at the estimates given in Table 2, we can compute the transcritical bifurcation points for each parameter. To simplify the calculations, we will treat \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _F$$\end{document} as a Heaviside step function, that is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Phi _F(x) = f_0 + (f_1-f_0)H(x - \theta ), \quad \text { where } \quad H(x) = {\left\{ \begin{array}{ll} 1 & x \ge 0, \\ 0 & x < 0, \end{array}\right. } \end{aligned}$$\end{document}

which can be viewed as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _F$$\end{document} in the limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_F \rightarrow 0$$\end{document} .

FGBA Model Steady States

We next solve for the FGBA steady states. We first use the relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = 1 - F -G - B$$\end{document} to eliminate the ash land cover fraction from (3), which gives the following system of equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{d}{dt}F(t)&= (1-F-B) \varphi F -\beta _F B F -\Phi _F(G)F- \mu F \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{d}{dt}G(t)&= \gamma (1-F-G-B) +\mu F- \varphi FG -\Phi _G(G) G - \beta _GBG \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{d}{dt} B(t)&=\Phi _G(G) G + \Phi _F(G)F + (\beta _G G + \beta _F F)B -q B . \end{aligned}$$\end{document}

From this point forward, we will set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = \varphi _A = \varphi _G$$\end{document} to simplify calculations. This is an ecologically reasonable assumption since the rate of forest spread into ash should not differ substantially from the rate of forest spread into grassland. Note that a similar assumption was made in Wuyts and Sieber (2023). Then after setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{G}= \dot{B}= \dot{A} = 0$$\end{document} we can solve for roots \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\overline{F}, \overline{G}, \overline{B}, \overline{A})$$\end{document} . As for the GBA subspace, the full model retains forward invariance of the ecologically relevant region of the phase-space (proof deferred to Appendix A).

Proposition 5

For any ecologically relevant initial value (i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(0) \ge 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(0) \ge 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B(0) \ge 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(0) + G(0) + B(0) \le 1$$\end{document} ), the solution (F(t), G(t), B(t)) remains ecologically relevant for all times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \in \mathbb {R}$$\end{document} .

We can now investigate the transcritical bifurcation anticipated to occur at the transition from the GBA steady state to states with non-zero forest. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{B}(F)$$\end{document} be a function that outputs the equilibrium burning proportion(s) B for any given input F. The function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{B}(F)$$\end{document} can be obtained from Eqs. (9) to (11) by solving for B in terms of F after setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{G}= \dot{B}= 0$$\end{document} and enforcing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F + G+ B+ A = 1$$\end{document} . In general, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{B}(\cdot )$$\end{document} is a complicated function we will analyze later. The vector field for forest cover can be expressed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} Q(F)&= (1 - F - \overline{B}(F))F \varphi - \beta _F \overline{B}(F)F - (f_1 + \mu )F \\&= -\overline{B}(F) F (\varphi +\beta _F) - \varphi F^2 + F(\varphi -\tilde{\mu }) \end{aligned}$$\end{document}

where we introduce the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\mu }:= \mu + f_1$$\end{document} to reduce notational clutter. The normal form for a transcritical bifurcation is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{x}=a_1x+x^2$$\end{document} . To verify that the genericity conditions for a transcritical bifurcation hold, we first note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q(F = 0) = 0$$\end{document} as expected since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0$$\end{document} is an equilibrium. Furthermore,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} Q'(F = 0)&= -\overline{B}(0)(\varphi + \beta _F) + \varphi - \tilde{\mu } \end{aligned}$$\end{document}

which indicates the stability of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0$$\end{document} equilibrium (i.e. the GBA equilibrium). For instance, note that the inequality \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q'(F) < 0$$\end{document} is exactly the same as Eq. (7) after making the simplifications \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = \varphi _G = \varphi _A$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _F(\overline{G}) = f_1$$\end{document} ; the latter assumption is motivated by assuming we are in the high-fire incidence regime with grass the dominant cover type (see Table 2 for the definition of the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1$$\end{document} ). Lastly, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} Q''(0)&= - 2\overline{B}'(0)(\varphi +\beta _F) -2\varphi \end{aligned}$$\end{document}

and it follows that as long as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q''(0) \ne 0$$\end{document} a non-degenerate transcritical bifurcation occurs at parameter values where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q'(F=0)$$\end{document} i.e. when

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \overline{B}(0)&= \frac{\varphi - \tilde{\mu }}{\varphi + \beta _F} \end{aligned}$$\end{document}

Next, we know that after passing the transcritical bifurcation, the system displays two distinct equilibria, one stable and the other unstable. Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0$$\end{document} is always an equilibrium, it remains to be determined whether or not the second equilibrium is ecologically plausible, i.e., contained in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {F}_{FGBA}$$\end{document} . To do this, we first note that the nonzero forest equilibrium can be calculated from Eq. (9) by dividing out the factor of F and setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{F} = 0$$\end{document} . This gives an implicit equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \overline{F}(\overline{B})&= \frac{(\varphi -\tilde{\mu }) - (\varphi + \beta _F)\overline{B}(\overline{F})}{\varphi }. \end{aligned}$$\end{document}

Next we can expand about the equilibria by writing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F = 0 \rightarrow \delta F$$\end{document} while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi \rightarrow \varphi + \delta \varphi $$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F \rightarrow \beta _F + \delta \beta _F$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\mu } \rightarrow \tilde{\mu }+ \delta \tilde{\mu }$$\end{document} . Then, after keeping only first-order terms, we find that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \bigg [\frac{\varphi }{\varphi + \beta _F} + \overline{B}'(0)\bigg ] \delta F = \bigg [ \frac{\tilde{\mu } + \beta _F}{(\varphi + \beta _F)^2}\bigg ] \delta \varphi = \bigg [ \frac{\tilde{\mu } - \varphi }{(\varphi + \beta _F)^2}\bigg ] \delta \beta _F = \bigg [\frac{-1}{(\varphi + \beta _F)^2}\bigg ] \delta \tilde{\mu }. \end{aligned}$$\end{document}

Then, the new equilibria emerging at the transcritical point after perturbation of a parameter value is ecologically plausible exactly when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta F > 0$$\end{document} . Thus, whether or not the non-zero forest equilibria are ecologically plausible can be determined based on the parameter values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\mu }, \beta _F, \varphi $$\end{document} as well as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{B}'(0)$$\end{document} .

Next, after some calculations (see Appendix B for details), we find that the parameter values along which the transcritical bifurcation occurs can be determined to various corrections in orders of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\varphi - \tilde{\mu }}{\varphi + \beta _F} = \frac{\gamma (\beta _G - q)}{\beta _G(q+\gamma )} + \frac{q g_1}{\beta _G(\beta _G - q)} + \frac{q(q+\gamma )g_1^2}{(q-\beta _G)^3 \gamma } + \mathcal {O}(g_1^3) \end{aligned}$$\end{document}

A plot of numerically computed branch points and Eq. (13) plotted to various orders in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0$$\end{document} is given in Fig. (4).Fig. 4A plot of the computed branch points governing the stable to unstable transition for the GBA steady state for system (4) alongside the analytically predicted branch points in Eq. (13). Branch points are plotted as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} , the rate of fire spread over forest, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} , the rate of forest spreading over grass and ash (assuming \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = \varphi _A = \varphi _G$$\end{document} ). Parameter values: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 30$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = 0.1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 1$$\end{document}

Bifurcation diagrams (partially timescale separated)

As the previous section shows, analyzing the FGBA system with a non-zero forest is complex. As such, a numerical study of the system is helpful for a more general understanding of the system’s behavior. Numerical continuation computations with timescale-separated parameter values are extremely slow and often challenging to perform. Thus, we first carried out our analysis with a set of “partially timescale-separated parameter values”. By partially timescale-separated, we mean that the values used likely underestimate the speed at which fire dynamics operate relative to vegetation dynamics and related processes (e.g., seed dispersal). This analysis can still provide valuable insight into overall system dynamics and serve as a baseline for further analysis of scenarios with more realistic levels of timescale separation (see Section 4).

Fig. (5) below illustrates the presence or absence of forest from stable steady-states of the FGBA model through one-parameter bifurcation diagrams. In each panel, the forest proportion is plotted on the y-axis against a different system parameter on the x-axis, with solid red lines indicating stable equilibria and solid black lines denoting unstable equilibria. Fig. (5)ABCDE all contain parameter regimes with multiple stable equilibria (bistability), with the only parameter not showing a bistable range being Fig. (5)F, which varies the tree mortality rate due to non-fire-related causes ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} ). These bistable ranges are often marked by saddle node bifurcations (denoted LP in the diagrams) at which new equilibria emerge or vanish, such as in Fig. (5)ABC. As expected, the zero forest cover state is stable at low values of the rate of forest spread into grass and ash ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = \varphi _A = \varphi _G$$\end{document} ) and high values of the rate of fire spread into forest ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} ). Another generic feature of these one-parameter bifurcation diagrams is the change of stability of the GBA subsystem via transcritical bifurcations (denoted BP in the diagrams) at which forest either enters or is excluded from the system. For example, in Fig. (5)D there is a stable forest-dominant steady-state for all values of the rate at which grass regrows from ash ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} ), and bistability emerges when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} reaches a critical value ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 22$$\end{document} ) at which the GBA steady state becomes stable. Similarly, a transcritical bifurcation is present at the critical parameter values of both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} where the GBA steady state changes stability (Fig. (5)BE). Additional two-parameter bifurcation analysis illustrating the size of the bistable regions as multiple parameters are varied together can be found in Appendix C.Fig. 5. One-parameter bifurcation diagrams for the mean-field FGBA system (3) with partially timescale-separated parameter values with the forest proportion of the solution of the y-axis. A. Varying the fire spread rate over grass ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} ). B. Varying the fire spread rate over forest ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} ). C. Varying the rate at which fire is quenched. D. Varying the rate at which grass regrows from ash ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} ). E. Varying the rate of forest spread over ash and grass ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} ). F. Varying the rate of tree mortality due to non-fire-related causes ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} ). Red lines indicate stable equilibria, while black lines indicate unstable equilibria. LP indicates a saddle-node bifurcation point (limit point), BP indicates a transcritical bifurcation point (branch point), and N indicates a neutral saddle equilibrium with no stability change. Parameter values (fixed if not the bifurcation parameter): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 30$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = \varphi _A = \varphi _G = 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.5$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 1$$\end{document}

Fig. (6) gives more detail on the structure of the steady-states shown in Fig. (5) by once more varying model parameters systematically but now plotting all four cover type proportions on the y-axis. In each case, the system is started from a very low forest cover proportion ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0.01$$\end{document} ) and run until it reaches a steady state. This serves to highlight the balance between grass, ash and burning states, and also the abrupt nature of the transitions possible when parameters are varied. For example, in Fig. (6)C, a small change in the rate of forest spread over ash and grass (around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi \approx 2.8$$\end{document} ) takes us from stable forest exclusion to a steady state in which forest occupies nearly the entire landscape.Fig. 6. Equilibrium land state proportions (i.e., steady-state values) of the mean-field FGBA system (3) reached starting from a low forest state ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0.01$$\end{document} initial condition) for partially timescale separated parameter values varying one system parameter at a time. A. Varying the fire spread rate over grass ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G$$\end{document} ). B. Varying the fire spread rate over forest ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} ). C. Varying the rate of forest spread over ash and grass ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} ). D. Varying the rate of tree mortality due to non-fire-related causes ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} ). Parameter values (fixed if not the bifurcation parameter): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 30$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.5$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 1$$\end{document}

Bifurcation diagrams (timescale separated)

Several one-parameter bifurcation diagrams were plotted using fully timescale-separated parameter values using MatCont (Dhooge et al. 2008) to illustrate the role of timescale separation and contrast with the results of the previous section. Due to prohibitively slow computation time, two-parameter bifurcation diagrams were not computed for fully timescale-separated parameter values.Fig. 7. One-parameter bifurcation diagrams for the mean-field FGBA system (3) for fully timescale-separated parameter values with the forest proportion, F, on the y-axis and the parameter being varied on the x-axis. A. Varying the fire spreading rate in forest ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} ). B. Varying the grass regrowth rate from ash ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} ). C. Varying the forest spreading rate into ash and grass ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} ). Red lines indicate stable equilibria while black lines indicate unstable equilibria. LP indicates a limit point, BP indicates a branch point, and N indicates a neutral saddle equilibrium. Parameter values (fixed if not the bifurcation parameter): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 20000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 500$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.05$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 0.1$$\end{document} Fig. 8A. One parameter bifurcation diagrams for the mean-field FGBA system (3) with the forest tree mortality rate ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} ) as the bifurcation parameter. Red lines indicate stable equilibria while black lines indicate unstable equilibria. BP indicates a branch point. B. Same as panel A with the linear regime near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0$$\end{document} overlaid as a solid blue line. C. Same as panel A with the linear regime near \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=1$$\end{document} overlaid as a solid blue line. Parameter values (fixed if not the bifurcation parameter): \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 25$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 20000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 500$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi = 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.05$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 0.1$$\end{document}

Fig. (7) shows one-parameter bifurcation diagrams for the forest steady-state as different system parameters are varied in the mean-field FGBA model (3) for the fully timescale separated parameter set (full parameter set included in the figure caption). We immediately note the stark differences between Fig. (7)A and Fig. (5)B, which both vary the fire spread rate over forest ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} ), but with very different levels of timescale separation between fire and vegetation dynamics. For the more realistic parameter set used in Fig. (7)A, we observe a qualitatively different picture of the forest dynamics as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} varies; we no longer have saddle-node bifurcations and a narrow bistable parameter regime, but instead we observe bistability between the GBA steady state and a forest steady state for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F$$\end{document} sufficiently large. We note that such a significant qualitative change is not entirely unexpected since we have moved numerous parameters in shifting to the fully timescale separated parameter set and thus we do not expect Fig. (7)A to be a simple smooth deformation of Fig. (5)B. However, Fig. (7)B and Fig. (7)C do qualitatively match the corresponding partially timescale separated bifurcation diagrams (Fig. (5)D and Fig. (5)E), albeit with the locations of the saddle node and branch points significantly shifted in parameter space. Overall, the results of Fig. (7) illustrate that the degree of timescale separation can significantly alter the underlying bifurcation structure of the model and increased timescale separation tends to enlarge the bistable regions of parameter space (at least from the perspective of varying any one parameter in isolation).

Fig. (8)A shows the one-parameter bifurcation diagram for the forest steady-state as the forest tree (non-fire) mortality rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} varies for the full timescale separated parameter set (cf. Fig. (5)F). However, the branch point occurs very close to zero so we set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 25$$\end{document} to ensure that the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} branch point occurs at a reasonable positive value and allows us to observe the trends in F. The steady states curves for F notably has two highly linear regimes in this case. In this parameter regime, we can deduce asymptotic formulae for these curves and these asymptotic approximations are shown as solid blue lines in Fig. (8)BC. The associated calculations and formulae for these approximating curves can be found in Appendix E.

Dynamics of the spatial stochastic FGBA model

Simulation algorithms and codes

The spatial FGBA model was simulated in MATLAB using the Gillespie algorithm (Gillespie 1977). All code used in this project can be accessed at https://github.com/patterd2/vegetation_fire_models. Testing and validation of the simulation codes are outlined in Appendix D.

Grassland without forest

Unless otherwise noted, all simulations in this section and the next section were run with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L = 1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 2000$$\end{document} . Parameter values were identical to those used in the timescale-separated bifurcation diagram analysis, with the widths of all spreading kernels set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _F=\sigma _B=\sigma _G=0.05$$\end{document} . The mean-field analysis predicts that the GBA model will exhibit a single ash-dominated stable fixed point. On longer timescales, we typically observe fluctuations about the mean-field dynamics in most parameter regimes, as expected, but the stochastic spatial model displays more complex dynamics on shorter timescales, and hence we focus on this case. Transient dynamics play a crucial role in ecological applications by capturing short- to medium-term ecosystem responses that differ from long-term equilibria, especially following disturbances or environmental change (Hastings 2001; Hastings et al. 2018). In the context of forest–savanna ecosystems, these dynamics may help to explain shifts in vegetation structure, species composition, and fire regimes observed in remotely sensed data on shorter timescales.

Using timescale-separated parameter values enables the spatial-stochastic GBA model (forest excluded) to exhibit various behaviors that can only be observed appropriately on vastly differing timescales. For example, the dynamics of the cover proportions of the different states at a fire scale vs a grass time scale are shown in Fig. (9)AB. Fig. (9)A shows the profile of a single ignition event over a day, while Fig. (9)B shows multiple small fires and grass regrowth phases taking place over several months. Fig. (9)C shows a montage of the spatial progression of a typical fire ignition and spreading event occurring over a very short time (several hours).Fig. 9. Simulations of the spatial-stochastic GBA model. A. System dynamics on the fire timescale (with time t in years), with the proportions of space covered by each state on the y-axis. Note that a slight increase in burning area (fire ignition) is followed by a large increase in ash and a large decrease in grass. The grass cover then gradually recovers from the ash state. B. System dynamics on a grass timescale showing repeated large fire events followed by grass recovery. C. Montage of fire ignition and spreading in grassland on a fire timescale as shown in A. D. Montage of grass regrowth following a fire ignition event in grassland on the grass timescale as shown in A. Parameter values: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 20000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 500$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _A=\varphi _G = 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.05$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 0.1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=2000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _B=0.035$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _F=\sigma _G=0.05$$\end{document}

On the fire timescale, the spatial structure of a fire front spreading through space and leaving behind a region of ash is readily apparent. On the grass timescale, the fire dynamics are obscured, but the regrowth of grass following fire spreading events is easily observed. Compared to the fire spreading, grass regrowth occurs relatively homogeneously spatially as shown in Fig. (9)D.Fig. 10. Simulations of the spatial-stochastic GBA model. A. Example of the system switching between large fires and mean-field-like behavior (i.e. fluctuations about an ash-dominated state). The y-axis is the proportion of space covered by each state and time t in years. B. Snapshot of the spatial system in the stochastic analogue of the ash-dominated mean-field GBA steady state at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.02$$\end{document} years and a snapshot of the system in the grass-dominated regime at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.042$$\end{document} years. The color legend used is shown in Fig. (9). Parameter values: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 15000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 500$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _A=\varphi _G = 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.05$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 0.1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=2000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _F=\sigma _G=\sigma _B=0.05$$\end{document}

In Fig. (9)B, several ignition events occur well separated in time on the grass timescale. This suggests that, in this parameter regime, the spatial GBA model spends a lot of time in a grass-dominated regime, distinct from the (unique) mean-field steady state, consisting of mostly ash. This is an especially interesting example of long-lived transient dynamics maintaining the system away from the mean-field prediction for extended periods. The spatial model can also more closely reflect the mean-field grass steady state under certain conditions. For example, the spatial model spends significant time near the mean-field GBA steady state for certain parameter values, which allows a uniform mixing of the states across the land area. For example, setting the fire quenching rate q to 20000 instead of 15000 produces the behavior shown in Fig. (10)A, with the system experiencing stochastic fluctuations that cause it to switch between regimes of big fire spreading events with spatial structure and mean-field behavior. The analogue of the mean-field GBA steady state is shown in Fig. (10)B for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.02$$\end{document} years, while Fig. (10)B with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.042$$\end{document} years shows a grass-dominated regime produced by transient and spatial dynamics not present in the mean-field model. A detailed investigation of the nature of these transient phenomena is beyond the scope of the present work, but it has been noted that high-dimensionality and spatial extent can promote long transients (Hastings et al. 2018), as observed here. In particular, the grass-dominated regime in Fig. (10)B takes time to build up a sufficiently dense and continuous layer of flammable cover to facilitate the large fire required to initiate the transition back to the ash-dominated state predicted by the mean-field. We hypothesize that this mechanism plays a key role in determining the duration of these extended transients.

Grassland with Forest

When studying grassland containing forest trees, we find additional new behavior at the forest timescale. In general, ignition events in grass cause large fires that propagate to an extent determined by the surrounding distribution of forest and destroy trees near the perimeter of the forest regions. During periods between fire events, the trees steadily regrow.

In general, there are two distinct outcomes over longer time periods. In the first scenario, the forest trees grow faster than the transitions that destroy them and eventually become so dense that fires can no longer propagate; thus, the high forest state is stably maintained. Alternatively, the grass fires are sufficiently destructive to reduce the forest proportion to a vanishing or near-vanishing proportion, where the forest can no longer spread noticeably. These two cases are illustrated in Fig. (11) and confirm that our model produces the desired forest-savanna bistability in its complete spatial stochastic form. As forest cover decreases, fires become larger and more frequent.Fig. 11. Two different possible outcomes for forest (dark green line) in the spatial-stochastic FGBA model; the proportion of space covered by each state is on the y-axis with time in years. Both simulations were run with the same parameter values and initial conditions. Forest was randomly distributed and occupied \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$40\%$$\end{document} of sites. The remaining sites were all grass. The color legend is shown in Fig. (9). Parameter values: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 20000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 500$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _A=\varphi _G = 1.5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.05$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 0.1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=2000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _B=0.03$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _F=\sigma _G=0.05$$\end{document}

Due to the stochastic nature of the model and the underlying bistable structure, we cannot determine the long-term cover proportions at the outset; we may even observe noise-induced switching (Weinan et al. 2021) between grass and forest-dominated steady states, which cannot be understood in terms of the mean-field approximation. In simulations starting at an initial condition of high forest cover, fire size also shows a noticeable increase as the forest cover decreases. A montage demonstrating the ability of forest cover to limit the extent of fire spreading is given in Fig. (12), corresponding to the fire occurring at the start of the simulation.Fig. 12A simulation of the spatial stochastic FGBA model demonstrating the ability of forest cover to limit the extent of fire spreading and specifically to stop smaller fires from growing into larger fires on the scale of the entire spatial domain (compare to Fig. (9)CD). The color legend is shown in Fig. (9) with time t in years. Parameter values: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _G = 50000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _F = 10000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q = 20000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 500$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _A=\varphi _G = 1.5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0 = g_0 = 0.01$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_1 = 0.05$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_1 = 0.1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=2000$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L=1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _B=0.03$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _F=\sigma _G=0.05$$\end{document}

The top panel of Fig. (11) highlights the dynamic nature of the open-canopy savanna-like state of the system, where frequent small to medium-scale fires are required to maintain the system state. Thus, our model in this setup could be used to distinguish between woody encroachment and natural fluctuations in the open canopy state. In contrast, we observe relatively small fluctuations in the forest-dominated state in the bottom panel of Fig. (11) and hence our initial modeling suggests that empirically observed fluctuations in forest cover are likely to be of greater concern. It may be possible to quantify this effect with more detailed parameter fitting and comparison with empirical data. However, other factors such as rainfall seasonality, which will impact fire patterns, would need to be added to the modeling structure to see how they amplify or dampen the aforementioned fluctuations (Accatino and De Michele 2013).

Conclusions

We have developed and analyzed a spatial stochastic model appropriate for studying fire-mediated forest-savanna bistability on short timescales where stochasticity and transient dynamics are expected to play crucial roles. Our work bridges the gap between highly resolved global vegetation models and minimalist models of forest-savanna bistability, such as the original Staver-Levin model. While the extant modeling literature has often relied on mean-field approximations (Touboul et al. 2018; Staver and Levin 2012), our work highlights the need to consider appropriately detailed models when estimating ecosystem resilience and stability from empirical data. The mean-field (spatially implicit) version of our FGBA model provides a reasonably simple method to study the qualitative behavior of the system, but does not necessarily provide a full description of the dynamics of the spatial stochastic model. For instance, in the absence of forest trees, the stability of the grassland steady state in parameter regimes with sufficiently high forest mortality is observed in both the mean-field and spatial models (Figs. (5F) and (9)). However, the mean-field model predicts only a single grassland steady state at unrealistically high ash coverage. In contrast, the stochastic spatial version of the model switches between a spatially homogeneous grassland regime with high ash coverage, mirroring the steady state predicted by the mean-field GBA model, and a new grass-dominated regime with occasional fire-spreading events well separated in time (Fig. (9)) not present in the mean-field model and not previously observed in comparable spatial stochastic models of forest-savanna ecosystems (Hébert-Dufresne et al. 2018; Wuyts et al. 2017). We hypothesize that these new transient dynamics are a product of the system’s high dimensionality and the ability of the spatial extent to block large fires until the grass layer has become sufficiently dense; these phenomena are beyond the scope of the mean-field model and hence require a separate, dedicated analysis to fully quantify.

Simulations of our spatial FGBA model explicitly demonstrate the hypothesized mechanisms underpinning bistability in forest tree cover. In particular, occasional ignitions followed by rapid fire spreading can maintain low forest cover even at high forest spreading rates. Meanwhile, the inability of fire to spread in regions of dense tree cover maintains high forest cover even at high fire ignition and spread rates (Fig. (12)), consistent with other recent spatial stochastic models of forest-savanna ecosystems (Hébert-Dufresne et al. 2018; Wuyts et al. 2017). In further simulations (see Appendix D.2 and Fig. (15)), the stability of the forest as a function of the parameter values showed qualitative agreement in both the mean-field and spatial models. In particular, increases in parameters such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} (the rate of fire spread over grass or ash) allow the system to have a finite probability of reaching a high forest state.

Our proposed FGBA model is highly versatile and can be used to study a wide range of possible forest and grassland setups. For example, differences in soil quality could be modeled by choosing vegetation sites within \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} according to a non-uniform probability distribution. One could also model forest spread via heavy-tailed, non-Gaussian spreading kernels, which may be a more accurate model of forest spread than our current assumption that forest trees only spread locally (Nathan et al. 2012). One could also investigate the impact of non-spatially uniform forest distributions, for example, if the forest were distributed into distinct regions of high tree density separated by areas of low tree density or into shapes with varying perimeter-area ratios. An essential addition to consider in future work and developments of this framework will be rainfall seasonality and its associated consequences for fire ignition. Forest-savanna ecosystems are subject to yearly wet and dry seasons that significantly impact tree growth and the relative humidity of the ecosystem, hence making the system more or less fire-prone at certain times of the year (Accatino and De Michele 2013). This will likely alter the dynamics presented in Section 4, but we believe the present analysis provides a crucial baseline against which to benchmark results that incorporate seasonality, isolating the new effects introduced by seasonal forcing.

Our model aims to alleviate several limitations of similar mathematical models of forest-savanna ecosystems, while retaining a reasonable degree of mathematical tractability. In particular, our model is spatially explicit, models dynamics on all time scales (fire, grass, and forest), allows an arbitrary distribution of vegetation sites, and incorporates arbitrary spreading kernels, enabling the study of non-isotropic or long-distance spreading. However, to avoid excessive mathematical and computational complexity, our model has some limitations. For instance, our model considers only fire and vegetation spread as transition processes in this paper, although our overarching framework can incorporate other processes. Other models consider additional factors, such as herbivory, soil layer quality, and competition between species (Holdo et al. 2022; Abades et al. 2014), which could be added to our mathematical framework at the cost of additional complexity. Moreover, we could make our model more realistic by expanding the state-space to include different savanna tree types at different life stages, as in other related models (Touboul et al. 2018). The key role of transient dynamics and stochasticity will likely remain in such extensions. Still, it would undoubtedly be important to extend the present work to these more realistic settings to relate more closely to empirical observations and hence further our understanding of these precious and endangered ecosystems.