Spatial modeling of forest-savanna bistability: Impacts of fire dynamics and timescale separation
Kimberly Shen, Simon Levin, and Denis D. Patterson [∗]
1 Department of Physics, Princeton University, USA 2 Peabody Institute, Johns Hopkins University, USA 3 Department of Mathematical Sciences, Durham University, UK 4 Department of Ecology and Evolutionary Biology, Princeton University, USA 5 High Meadows Environmental Institute, Princeton University, USA
May 5, 2025
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. 1 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. Savannas and forests are frequently observed in the tropics,
∗ Corresponding author: denis.d.patterson@durham.ac.uk while intermediate tree cover is rarely observed. 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. 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. Therefore, a model that accurately predicts the responses of forest-savanna mosaics to various environmental disturbances would usefully inform various climate- and biodiversityrelated interventions, such as selecting and maintaining protected areas and regulations of human ecological disturbances. As such, a maximally valuable model must encompass shortterm 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. Thus, regions with sufficiently high grass cover may burn frequently enough to maintain savanna conditions despite climatological conditions supporting forest growth. The original Staver-Levin model and related models, which incorporate fire-vegetation interactions in a spatially implicit mean-field approximation, have been widely used to explain and study forest-savanna bistability. 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 keep them interpretable and analytically tractable and only indirectly incorporate spatial processes. For example, 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. 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, 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, but such models cannot transparently demonstrate the underlying dynamics of the system nor be applied to novel situations. Other studies based on field sampling of species diversity, fire frequency and/or soil quality have similar drawbacks. Some mean-field models have been extended to incorporate spatial effects by adding diffusion terms to the nonspatial dynamics. However, these models are limited by an inability to model anisotropic forest seed dispersal and do not explicitly model fire dynamics. Explicit fire dynamics are essential for examining the behavior of grass-forest systems on shorter timescales and studying the effects of seasonalitydependent flammability on forest dynamics. 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. Most spatially explicit forest fire models employ a square grid and use nearest neighbordependent transition probabilities to model spatial processes such as fire spread and seed dispersal. These models have already led to notable findings. For example, explicitly demonstrated that fire-vegetation feedback can maintain bistability in tree cover states, found that forest-savanna mosaics occur under limited conditions, suggesting that such mosaics are vulnerable to loss under climate change, and 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. The paper is organized as follows: Section 2 introduces the spatial stochastic modeling framework that forms the basis for this study, including 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. 2 A Spatial Stochastic Forest Fire Model
Ecologically, tropical forests differ markedly from forests of coniferous trees 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 or have thickened bark to prevent stem death. Second, fire does not propagate readily through tropical forests but does propagate rapidly through savanna grassland. Mechanisti cally, this is due to forest understory shade excluding flammable C4 grasses in addition to reduced wind speeds and increased moisture in the forest microclimate. Lastly, tropical forests have relatively open canopies that do not entirely shade out grass. For brevity, we will henceforth refer to tropical forests as “forest” and savanna regions with lower tree density as “grassland”.
2.1 Mathematical Description of the Model
We use a spatially extended Markov jump process to model the dynamics of fires in forestgrass systems using the mathematical framework introduced in. The set Ω ⊂ R
denotes the spatial domain and the random variables {r i ∈ Ω: i = 1 , . . ., N } are sites (locations) where vegetation is present or has the potential to grow. Each site i, with respective location r i, undergoes a Markov jump process with possible states F (forest), G (grass), B (burning), and A (ash).
We thus refer to models with this state space as the “FGBA model” for brevity. The transition rates between these states depend on the spatial locations and states of all the other sites at locations {r j ∈ Ω: j = 1 , . . ., N, j ̸ = i} . This allows us to model both spatial spreading processes (e.g. forest and fire spreading) as well as non-spatial spontaneous transitions. The spatial locations of the vegetation can be drawn from any sufficiently regular distribution on Ω, but, in this paper, we draw from a uniform distribution on Ω to reflect a homogeneous domain scenario. In the limit as N →∞, 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. 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. The parameters in the transition rates are defined as follows:
- φ G , φ A are constants controlling the rate of forest seeding into grass and into ash, respectively. - β F , β G are constants controlling the rate of fire spread within forest and within grass, respectively. - W F , W B : R [+] → R [+] 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 C [∞] . Transition Transition rate at site i Ecological process
G → F φN G ~~�~~ Nj =1 [W] [F] [ (] [r] [i] [ −] [r] [j] [)] {X [j] ( t )= F } forest spreading into grass
A → F φN A ~~�~~ Nj =1 [W] [F] [ (] [r] [i] [ −] [r] [j] [)] {X [j] ( t )= F } forest spreading into ash
β F N F → B N ~~�~~ j =1 [W] [B] [(] [r] [i] [ −] [r] [j] [)] {X [j] ( t )= B} fire spreading into forest
β G N G → B N ~~�~~ j =1 [W] [B] [(] [r] [i] [ −] [r] [j] [)] {X [j] ( t )= B} fire spreading into grass
F → G µ non-fire forest mortality
B → A q fire quenching
A → G γ grass regrowth from ash
Table 1: Transition rates between states and their corresponding ecological processes. For the neighbor (local) transition rate parameters for fire and forest invasion ( φ G, φ A, β F, β G ), 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 ( µ, q, γ ). We normalize the neighbor spread rates by [area] N [(][Ω][)] to ensure that the spreading rates remain bounded in the limit N →∞ . In addition to the diffusive mode of fire spread described by the W B ( · ) kernel and β G , β F parameters, our modeling framework also includes a “cascade” mode of vegetation burning motivated by percolation theory. Assuming that fire spread through the grass-forest land is well-approximated as a percolation process, the flammability of any grass site will increase markedly once the fraction of grass sites in the nearby vicinity of the grass site surpasses the critical percolation threshold; this mechanism is also well supported by empirical data. We incorporate this effect by adding the terms
1 � N
N � W G ( r i − r j ) 1 {X j ( t )= G}
i =1
�
N � W G ( r i − r j ) 1 {X j ( t )= G}
i =1
and Φ F �
Φ G
1 � N
to the G → Band F → B transition rates of site i at time t . Here W G : R [+] → R [+] is another smooth spreading kernel and Φ G , Φ F : R [+] → R [+] are smooth sigmoidal functions which output flammability as a function of the local grass coverand the kernel W G . More explicitly, we assume that Φ G ( · ) and Φ F ( · ) will take the form g 1 − g 0 f 1 − f 0 Φ G ( x ) = g 0 + and Φ F ( x ) = f 0 + 1 + e [−] [(] [x][−][θ] [G] [)] [/s] [G] 1 + e [−] [(] [x][−][θ] [F] [ )] [/s] [F]
where the parameters are defined as follows: - f 0 , g 0 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. - f 1 , g 1 are the total flammability of a grass or forest site, respectively in the limit when the land patch is purely grassland. - θ F , θ G are the percolation thresholds in forest and grass, respectively. - s F , s G 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 s F = s G = 0 . 05. 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 yr [−], and we assume that the vegetation sites are separated by an average distance of 10 m. The vegetation sites are chosen randomly within a compact square domain ×
⊂ R for some fixed L > 0. We use periodic boundary conditions to reduce boundary effects and to model an infinite domain.
The distance between two points r, s ∈ × is then computed as
|r − s| = [L]
2 π
[−][s][x] [)] 2 [π] [(] [r][y] [−][s][y] [)]
L + arg e [i] L �� � �
[π] [(] [r][x][−][s][x] [)] arg e [i] L �� �
[y] [−][s][y] [)] 2
L �� �
where the argument function has range[r] [j] [ −] [r] [i] [|] 2 πσ � 2 σ
-
�
where □ = G, F, or B . To preserve the physical meaning of the three parameters σ - under changes in Land/or N, we choose the values of σ - to be in units of the average spacing L between the vegetation sites in the domain which we define as ∆ x ≡ ~~√~~ N [. The normalization]
constant of W - was chosen that
W - ( r j , r i ) dr i ≈ Ω �
R [ W] [□] [(] [r] [j] [, r] [i] [)] [dr] [i]
1 lim N →∞ N
N �
� W - ( r j , r i ) =
� i =1
where the approximation is valid when σ - ≪ L . The model described in this section is similar to previous models proposed by HebertDufresne et al. and Wuyts and Sieber, particularly the choice of the state space, but has several distinguishing features. Firstly, our model allows the sites to be distributed arbitrarily in Ω, while the related models use a square lattice of sites. Secondly, and include additional spontaneous transitions G, A, → F to model the homogeneous longdistance dispersal of forest seeds. Our model more realistically incorporates long-distance forest seed dispersal by allowing an appropriately heavy-tailed kernel for W F ( · ). 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. 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. 2.1.1 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 L distribution, is ∆ x ≡ ~~√~~ N [. We now estimate the parameter values by considering the case]
where ∆ x = 1 m and L = 100 mare fixed, requiring N = 10 . - A square patch of grassland of side length ∆ x surrounded by fire will burn after several minutes. Using Table 1 and Eq., the burning rate can be approximated as β G then β G ≈ 10 yr [−] . - A square patch of forest of side length ∆ x, which is surrounded by fire, will burn after about an hour. The burning rate can be approximated as β F so β F ≈ 10 yr [−] . - A burning site is expected to burn for several hours before turning into ash, i.e. q ≈ 10 yr [−]
- An ash site is expected to regrow grass in several months i.e. γ ≈ 10 yr [−] . - A square patch of grass or ash of side length ∆ x surrounded by forest will grow a forest tree after several years to a decade i.e. φ A ≈ φ G ≈ 10 to 10 [−] yr [−] . - A tree at a square forest site of side length ∆ x will spontaneously die from non-fire related causes after about 100 years i.e. µ ≈ 10 [−] yr [−] . - A square grass or forest vegetation site of side length ∆ x completely surrounded by forest will spontaneously catch fire once every 100 years i.e. f 0 , g 0 ≈ 10 [−] . - A grass or forest vegetation site will have increased flammability when completely surrounded by grass compared to a forest. We set f 0 = 10 [−] yr [−] and g 0 = 10 yr [−] . - The standard deviation of afire patch spreading throughout grass or forest vegetation is approximately 5 m, so set σ F ≈ 5 m. - The standard deviation of a forest spreading throughout grass or ash is approximately 5 m, so set σ B ≈ 5 m. - The flammability of a vegetation site is dependent on grass proportions in a 50 m radius up to one standard deviation, so σ G ≈ 50 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. Parameter Value Units Ecological Interpretation
β G 10 yr [−] magnitude of fire spread rate over grass
β F 10 yr [−] magnitude of fire spread rate over forest q 10 yr [−] rate at which fire is quenched
γ 10 yr [−] rate at which grass regrows from ash
φ A 1 yr [−] magnitude of forest spread rate over ash
φ G 1 yr [−] magnitude of forest spread rate over grass
µ 10 [−] yr [−] tree mortality rate due to non-fire-related causes f 0 10 [−] yr [−] spontaneous fire rate for each forest site when G ( t ) ≈ 0
f 1 10 [−] yr [−] spontaneous fire rate for each forest site when G ( t ) ≈ 1
g 0 10 [−] yr [−] spontaneous fire rate for each grass site when G ( t ) ≈ 0
g 1 10 [−] yr [−] spontaneous fire rate for each grass site when G ( t ) ≈ 1
σ F 5 m width of forest spreading kernel
σ B 5 m width of fire spreading kernel
σ G 50 m width of spontaneous fire kernel
Table 2: Estimates and representative values of the parameters for the FGBA model in units of yr [−]
There are three levels of time separation in the rate parameters given by β G , β F , q > γ ≫ φ A , φ G , µ, f 0 , g 0 . This concludes our mathematical description of the spatial FGBA model. We summarize the model in the state transition diagram in Fig., showing the relevant parameters and functions governing transitions between states at each vegetation
site. 2.2 Mean-Field Approximations of the Spatial FGBA Model
Theorem 2.2.1 in shows that the spatial FGBA model with uniformly randomly distributed sites on Ω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 µ
F G
φ G
φ A Φ F , β F Φ G , β G
γ
AB
q
Figure 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.
all other locations. Let X ( r, t ) indicate the probability that the limiting McKean-Vlasov process at location rand 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):
d
W F ( r − r [′] ) F ( r [′] , t ) dr [′] − µF ( r, t )
dt [F] [(] [r, t] [) = (] [φ] [G] [G] [(] [r, t] [) +] [ φ] [A] [A] [(] [r, t] [))] � Ω
− β FF ( r, t ) W B ( r − r [′] ) B ( r [′] , t ) dr [′] − Φ F � Ω
W G ( r − r [′] ) G ( r [′] , t ) dr [′] F ( r, t ) ,
�� Ω �
d
W F ( r − r [′] ) F ( r [′] , t ) dr [′] + µF ( r, t )
dt [G] [(] [r, t] [) =] [ γA] [(] [r, t] [)] [ −] [φ] [G] [G] [(] [r, t] [)] � Ω
− Φ G
��
W G ( r − r [′] ) G ( r [′] , t ) dr [′] G ( r, t ) − β G G ( r, t ) Ω � �
W B ( r − r [′] ) B ( r [′] , t ) dr [′] , Ω
d
W G ( r − r [′] ) G ( r [′] , t ) dr [′] − qB ( r, t )
dt [B] [(] [r, t] [) = (Φ] [F] [ F] [(] [r, t] [) + Φ] [G] [G] [(] [r, t] [))] �� Ω �
+ ( β G G ( r, t ) + β FF ( r, t )) W B ( r − r [′] ) B ( r [′] , t ) dr [′] , � Ω
d
W F ( r − r [′] ) F ( r [′] , t ) dr [′] ,
dt [A] [(] [r, t] [) =] [ qB] [(] [r, t] [)] [ −] [γA] [(] [r, t] [)] [ −] [φ] [A] [A] [(] [r, t] [)] � Ω
for each r ∈ Ω. Moreover, F ( r, t ) + G ( r, t ) + B ( r, t ) + A ( r, t ) = 1 for each ( r, t ) ∈ Ω × R + . 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 that takes the form of a system of ordinary differential equations. The general structure of the solution space of the spatial model 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. 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, dis crete and finite-size effects, which we expect to be especially important for short-time and transient dynamics, which we can contrast against the nonspatial mean-field model. To obtain a nonspatial (spatially implicit) mean-field model from, we assume that the types of states are well-mixed throughout the spatial domain Ω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 Ωand then define the land cover proportions F ( t ) := � Ω [F] [(] [r] [′] [, t] [)] [dr] [′] [/] [area(Ω),] [ G] [(] [t] [) :=] � Ω [G] [(] [r] [′] [, t] [)] [dr] [′] [/] [area(Ω), and so on.
The system]
Ω [F] [(] [r] [′] [, t] [)] [dr] [′] [/] [area(Ω),] [ G] [(] [t] [) :=] �
F ( t ) := � Ω [F] [(] [r] [′] [, t] [)] [dr] [′] [/] [area(Ω),] [ G] [(] [t] [) :=] � Ω [G] [(] [r] [′] [, t] [)] [dr] [′] [/] [area(Ω), and so on. The system]
then simplifies to the system of ordinary differential equations:
 d  dt [F] [(] [t] [) = (] [φ] [G] [G] [ +] [ φ] [A] [A] [)] [F][ −] [β] [F] [ BF][ −] [Φ] [F] [ (] [G] [)] [F][ −] [µF]   d
dt [G] [(] [t] [) =] [ γA] [ +] [ µF][ −] [φ] [G] [FG][ −] [Φ] [G] [(] [G] [)] [G][ −] [β] [G] [BG]
 d
[B] [F] [B]
d dt [F] [(] [t] [) = (] [φ] [G] [G] [ +] [ φ] [A] [A] [)] [F][ −] [β] [F] [ BF][ −] [Φ] [F] [ (] [G] [)] [F][ −] [µF]
d dt [G] [(] [t] [) =] [ γA] [ +] [ µF][ −] [φ] [G] [FG][ −] [Φ] [G] [(] [G] [)] [G][ −] [β] [G] [BG]
    
d dt [B] [(] [t] [) = Φ] [G] [(] [G] [)] [G] [ + Φ] [F] [ (] [G] [)] [F] [ + (] [β] [G] [G] [ +] [ β] [F] [ F] [)] [B][ −] [qB]
d dt [A] [(] [t] [) =] [ qB][ −] [γA][ −] [φ] [A] [FA.]
The state transition diagram for the spatial FGBA model in the mean field limit is
Φ( G ) G qB γA illustrated in Fig.. First notice that due to the G B A G cycle then if any one of G, B, or A is nonzero then all three must be nonzero. Next, since there is an
µFF G transition then F > 0 implies G > 0. Thus, we can classify all steady states as either GBA or FGBA. µF
F G
φ G FG
Φ F ( G ) F + β F BF
φ A FA FF Φ G ( G ) G + β G BG
γA
AB
qB
Figure 2: 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. 3 Analysis of the nonspatial FGBA model
Before studying the full spatial stochastic FGBA model introduced above, we perform some qualitative stability analysis of solutions to the nonspatial mean-field version of the model given by. 3.1 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 F = 0 and use the relation A = 1 − G − B to eliminate A from to reduce the system to
dtd [G] [(] [t] [)] = γ − Φ G ( G ) G − β G BG � dtd [B] [(] [t] [)] = Φ G ( G ) G + β G GB − qB. 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 ≥ 0 , B ≥ 0 , G + B ≤ 1 ), the solution ( G ( t ) , B ( t )) for Eqs. remains ecologically relevant for all times t ∈ R . Next, we solve for the steady-state equilibria of the system by setting both equations in System equal to 0. Solving the system then gives the condition
�
Φ G ( G ) = [γq]
γ + q
1 � G
1
G [−] [β] q [G]
q
for a steady state. Unfortunately, closed-form expressions of the equilibria cannot be obtained except in special cases due to the sigmoidal properties of the function Φ G, but we can readily deduce the uniqueness of the steady state:
Proposition 2. For any positive parameter values γ, q, and β G there exists a unique GBA steady state. The proof of Proposition 2 is deferred to Appendix A in the interests of brevity. 3.1.1 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 1. The GBA steady state is stable to perturbations within the F = 0 subspace for any parameter values. Proof. We compute the linearization matrix at F = 0 of System 4 which gives
J =

G [(] [G] [)] [G][ −] [Φ] [G] [(] [G] [)] [ −] [β] [G] [B] −γ − β G G  [−][γ] Φ [ −] [′] G [Φ][(] [G] [′] [)] [G] [ + Φ] [G] [(] [G] [) +] [ β] [G] [B] −q + β G G
Φ [′] G [(] [G] [)] [G] [ + Φ] [G] [(] [G] [) +] [ β] [G] [B] −q + β G G
 . The fixed point ( G, B ) is stable to perturbations within the F = 0 boundary exactly when J F =0 ( G, B ) has eigenvalues with strictly negative real parts. By the Routh-Hurwitz stability criteria this occurs if and only if Tr) < 0 and det) > 0. Note that since G < βq G [then] [ −][q] [ +] [ β] [G] [G <] [J] [F] [(] [G, B] [))] [ <] condition det) > 0 is equivalent to
~~2~~ G
γq Φ [′] ( G ) > ( γ + q )
β � q
β
q [−] G
= F [′] ( G ) �
which clearly holds from graphical inspection or more rigorously by recalling from our earlier argument that Φ [′] G [(] [G] [)] [ ≥] [ F] [′] [(] [G] [)] [ <]
Proposition 3. The GBA steady state is stable to invasions by F exactly when
β F B + Φ F ( G ) + µ > φ G G + φ AA. Proof. We compute the linearization matrix for the full FGBA system and evaluate it at F = 0. We use A = 1 − F − G − B to obtain a dynamical system in the three variables F, G, and B which gives the matrix
J
φ G G + φ AA − β F B − Φ F ( G ) − µ 0 0

. 
=

 
−γ − φ G G + µ −γ − Φ [′] G [(] [G] [)] [G][ −] [Φ] [G] [(] [G] [)] [ −] [β] [G] [B] −γ − β G G
Φ F ( G ) + β F B Φ [′] G [(] [G] [)] [G] [ + Φ] [G] [(] [G] [) +] [ β] [G] [B] −q + β G G
The GBA steady state is stable to invasion by F if all the eigenvalues of J have negative real parts. We have already shown that the eigenvalues of the 2 × 2 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. An intuitive interpretation of equation 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. 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 β G, β F, φ = φ G = φ A and µ using partially timescale-separated parameter values are given in Fig.. As expected, when φ is increased past a critical value or when β G, β F, and µ are decreased past a critical value, the proportion of forest sites within Ω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 ashdominated, spatially homogeneous GBA steady state. 3.1.2 Estimates of the GBA Steady State
We will now show that using our estimates of the parameter values, the steady state ( G, B, A ) can be easily estimated. Notice that F ( G ) becomes asymptotic to γγq + q []
G [−] [β] [G]
( G, B, A ) can be easily estimated.
Notice that F ( G ) becomes asymptotic to γγq + q [(] G [−] q [G] [)]
γq as G → 0 and the coefficient γ + q [G] [(] [G] G ∈ then we can approximate the condition Φ G ( G ) = F ( G ) as 0 = F ( G ) which gives q G ≈ β G [. 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 β G, i.e. decreasing the rate of fire spread through grass, both of which make sense intuitively. By solving System 4 and using G + B + A = 1 we can compute an estimate of the GBA steady state as well as the ratio B/A :
q q γ]
β G , q + γ
q + γ
⇐⇒ [B] � A
q ( G, B, A ) ≈ � β
q q − β G [)] q]
q + γ , q + γ
A [≈] [γ] q
q [.]
Based on the estimate, the ratio between the portions of land in burning and ash states is determined by γ, the rate of grass regrowth from ash, and q, the rate of fire quenching in a manner that agrees with intuition. Then Eq. can be used to give Eq. purely in terms of the model parameters:
q β F γ]
− β G q
+ Φ F q + γ � β
q + γ . β G
� + µ > [qφ] [G]
q
[G] qφ A]
+ β G q + γ
In particular, if we vary each of the parameters β F , γ, f 0 , µ, and φ = φ A = φ Gone 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 Φ F as a Heaviside step function, that is
Φ F ( x ) = f 0 + H ( x − θ ) , where H ( x ) =
which can be viewed as Φ Fin the limit s F → 0. 3.2 FGBA Model Steady States
1 x ≥ 0 , �0 x < 0 ,
We next solve for the FGBA steady states. We first use the relation A = 1 − F − G − B to eliminate the ash land cover fraction from, which gives the following system of equations:
d dt [F] [(] [t] [ −] [F][ −] [B] [)] [φF][ −] [β] [F] [ BF][ −] [Φ] [F] [ (] [G] [)] [F][ −] [µF]
d dt [G] [(] [t] [) =] [ γ] [ −] [F][ −] [G][ −] [B] [) +] [ µF][ −] [φFG][ −] [Φ] [G] [(] [G] [)] [G][ −] [β] [G] [BG]
d dt [B] [(] [t] [) = Φ] [G] [(] [G] [)] [G] [ + Φ] [F] [ (] [G] [)] [F] [ + (] [β] [G] [G] [ +] [ β] [F] [ F] [)] [B][ −] [qB.]
From this point forward, we will set φ = φ A = φ G 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. Then after setting G [˙] = B [˙] = A [˙] = 0 we can solve for roots ( F, G, B, A ). 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 4. For any ecologically relevant initial value ≥ 0 , G ≥ 0 , B ≥ 0 , F + G + B ≤ 1 ), the solution ( F ( t ) , G ( t ) , B ( t )) remains ecologically relevant for all times t ∈ R . 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 B ( F ) be a function that outputs the equilibrium burning proportion(s) B for any given input F . The function B ( F ) can be obtained from Eqs. to by solving for Bin terms of F after setting G [˙] = B [˙] = 0 and enforcing F + G + B + A = 1. In general, B ( · ) is a complicated function we will analyse later.
The vector field for forest cover can be expressed as
Q ( F ) =) Fφ − β F B ( F ) F − F
˜ = −B ( F ) F ( φ + β F ) − φF + F ( φ − µ )
where we introduce the parameter ˜ µ := µ + f 1 to reduce notational clutter. The normal form for a transcritical bifurcation is ˙ x = a 1 x + x . To verify that the genericity conditions for a transcritical bifurcation hold, we note first that Q = 0 as expected since F = 0 is an equilibrium. Furthermore,
Q [′] = −B( φ + β F ) + φ − µ ˜
which indicates the stability of the F = 0 equilibrium (i.e. the GBA equilibrium). For instance, note that the inequality Q [′] ( F ) < 0 is exactly the same as Eq. after making the simplifications φ = φ G = φ A and Φ F ( G ) = f 1 . Lastly, we have
~~′~~ Q [′′] = − 2 B( φ + β F ) − 2 φ
and it follows that as long as Q [′′] ̸ = 0 a non-degenerate transcritical bifurcation occurs at parameter values where Q [′] i.e. when
B = [φ][ −] [µ] [˜]
φ + β F
Next, we know that after passing the transcritical bifurcation, the system displays two distinct equilibria, one stable and the other unstable. Since F = 0 is always an equilibrium, it remains to be determined whether or not the second equilibrium is ecologically plausible, i.e. contained in FFGBA . To do this, we first note that the nonzero forest equilibrium can be calculated from Eq. by dividing out the factor of F and setting F [˙] = 0. This gives an implicit equation
[β] [F] [)] [B] [(] [F] [)] F ( B ) = [(] [φ][ −] [µ] [˜][)] [ −] [ δF = δφ = δβ F = δµ. � φ + β F � � ( φ + β F ) � � ( φ + β F ) � � ( φ + β F ) �
˜ µ − φ δφ = � � ( φ + β F )
− 1 δβ F = � � ( φ + β F )
δµ. ˜ �
Then, the new equilibria emerging at the transcritical point after perturbation of a parameter value is ecologically plausible exactly when δF > 0. Thus, whether or not the non-zero forest equilibria are ecologically plausible can be determined based on the parameter values µ, β ˜ F , φ as well as B ~~′~~. 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 g 1 :
˜ φ − µ qg 1 [γ] [)] [g] 1 = [γ] [(] [β] [G] [ −] [q] [)] [q] [] φ + β F β G ( q + γ ) [+] β G ( β G − q ) [+] ( q − β G ) γ [+] [ O] [ plotted to various orders in g 0 and f 0 is given in Fig.. 3.2.1 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.
Unfortunately, numerical computations with timescale-separated parameter values is extremely slow. Thus, we performed a numerical analysis with somewhat less timescale-separated parameter values: β G = 50, β F = 10, q = 30, γ = 10, φ = 1, µ = 0 . 01, f 0 = g 0 = 0 . 01, f 1 = 0 . 5, and g 1 = 1. As expected, the zero forest cover state is stable at low values of φ and high values of β F . Notably, a stable high forest state is also present in the same parameter ranges as the low forest states. Further, as expected, a transcritical bifurcation is present at the critical parameter values of both φ and β F where the GBA steady state becomes unstable. In the two-parameter bifurcation diagram in Fig., the regions enclosed by the blue saddle-node curves indicate the parameter regimes in which we can expect bistability in forest cover. 3.2.2 Bifurcation diagrams (timescale separated)
Several one-parameter bifurcation diagrams were plotted using fully timescale-separated parameter values using MatCont. The values used were β G = 50000, β F = 10000, q = 20000, γ = 500, φ = 1, µ = 0 . 01, f 0 = g 0 = 0 . 01, f 1 = 0 .
05, and g 1 = 0 . 1 for the bifurcation diagrams in Fig.. Due to prohibitively slow computation time, twoparameter bifurcation diagrams were not computed for fully timescale-separated parameter values. The same parameter values were used for the µ bifurcation diagrams in Fig. except that β F = 25 was set to ensure that the µ branch point occurs at a reasonable positive value. The µ bifurcation notably has two highly linear regimes which can be computed to high accuracy as demonstrated in Appendix (D). Figure 3: A plot of the computed branch points governing the stable to unstable transition for the GBA steady state alongside the analytically predicted branch points in Eq.. Branch points are plotted as a function of β F, the rate of fire spread over forest, and φ, the rate of forest spreading over grass and ash (assuming φ = φ A = φ G ). Parameter values: β G = 50, β F = 10, q = 30, γ = 10, φ = 0 . 1, f 0 = g 0 = 0 .
01, f 1 = 0 . 5, g 1 = 1
4 Dynamics of the spatial stochastic FGBA model
4.1 Simulation algorithms and codes
The spatial FGBA model was simulated in MATLAB using the Gillespie algorithm. All code used in this project can be accessed at `` `[firemodels]` . Testing and validation of the simulation codes is outlined in Appendix C. 4.2 Grassland without forest
Unless otherwise noted, all simulations in this section and the next section were run with L = 1 and N = 2000. Parameter values were identical to those used in the timescale-separated bifurcation diagram analysis, with the widths of all spreading kernels set to σ F = σ B = σ G = 0 . 05. The mean-field analysis predicts that the GBA model will exhibit a single ashdominated 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 Figure 4: One parameter bifurcation diagrams for partially timescale separated parameter values. 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. 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. 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 GBA simulation 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 afire scale vs a grass time scale are shown in Fig.AB. Fig.A shows the profile of a single ignition event over a day, while Fig.B shows multiple small fires and grass regrowth phases taking place over several months. Fig.C shows a montage of the spatial progression of atypical fire ignition and spreading event occurring over a very short time (several hours). On the fire timescale, the spatial structure of afire 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.D. In Fig.B, several ignition events occur well separated in time on the grass timescale. This suggests that, in this parameter regime, the spatial GBA model possesses a largely grass quasi-steady state, distinct from the 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 Figure 5: Equilibrium land state proportions reached starting from a low forest state for partially timescale separated parameter values near branch points. also demonstrate the mean-field grass steady state under certain conditions. For example, the spatial model can reach 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 20 , 000 instead of 15 , 000 produces the behavior shown in Fig.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.B for t = 0 .
02, while Fig.B with t = 0 . 042 shows a grass-dominated quasi-steady state 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, as observed here. In particular, the grass-dominated quasi-steady state in Fig.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 meanfield. We hypothesize that this mechanism plays a key role in determining the duration of these extended transients. 4.3 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 Figure 6: Two parameter bifurcation diagrams of the limit points in Fig. for partially timescale separated parameter values. CP indicates a codimension two cusp point where two curves of saddle-node bifurcations intersect. Figure 7: One-parameter bifurcation diagrams in β F (fire spreading rate in forest), γ (grass regrowth rate from ash) and ϕ (fire spreading rate in ash and grass) for fully timescaleseparated parameter values. 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. 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. and confirm that our model produces the desired forest-savanna bistability in its complete spatial stochastic form. As forest cover decreases, fires become largerand more frequent. 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 noiseinduced switching 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., corresponding to the fire occurring at the start of the simulation. Figure 8: A. Bifurcations in µ (non-fire forest mortality rate) for fully timescale separated parameter values with β F = 25. Red lines indicate stable equilibria while black lines indicate unstable equilibria. BP indicates a branch point. The linear regimes near F = 0 and F = 1 are plotted in blue in panels Band C, respectively. Figure 12: A simulation 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 domainCD). The color legend is shown in Fig.. The top panel of Fig. highlights the dynamic nature of the open-canopy savannalike 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. 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. 5 Conclusions
We have developed and analyzed a spatial stochastic model appropriate for studying firemediated 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. The mean-field (spatially implicit) version of our FGBA CD
Figure 9: A. System dynamics on the fire timescale 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 big fire events followed by grass recovery. C.
Montage of fire ignition and spreading in grassland on afire timescale as shown in A. D. Montage of grass regrowth following afire ignition event in grassland on a grass timescale as shown in A. 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 under parameter regimes with sufficiently high forest mortality is observed in both the mean-field and spatial models. 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 steady state with high ash coverage and an additional quasi-steady state characterized by high grass coverage and occasional fire-spreading events well separated in time). We hypothesize that these new transient dynamics are a product of the high-dimensional nature of the system 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 quantify fully. In a related example, the stability of the forest as a function of the parameter values showed Figure 10: 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. B. Snapshot of the spatial system in the stochastic analogue of the ash-dominated mean-field GBA steady state at time t = 0 . 02 and a snapshot of the system in the grass-dominated quasi-steady state at time t = 0 . 042.
The color legend used is shown in Fig.. qualitative agreement in both the mean-field and spatial models with increases in parameters such as φ (the rate of fire spread over grass or ash), allowing the system to have a finite probability of reaching a high forest state. While the extant modeling literature has often relied on mean-field approximations, our work highlights the need to consider appropriately detailed models when estimating ecosystem resilience and stability from empirical data. In addition, 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. 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 Ω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. 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. This will likely alter the dynamics presented in section 4, but we believe the present analysis is a crucial baseline against which to benchmark results incorporating seasonality to isolate 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 in the spatialstochastic FGBA model; the proportion of space covered by each state is on the y-axis. Both simulations were run with the same parameter values and initial conditions. Forest was randomly distributed and occupied 40% of sites. The remaining sites were all grass. The color legend is shown in Fig.. 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, 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. 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. Statements and Declarations
Acknowledgements
Denis Patterson and Simon Levin are grateful for the support of NSF DMS-1951358.
Competing Interests
The authors have no competing interests to declare. Data Availability
The codes used in this paper are available at ``
`[firemodels]` . Appendix A Mathematical Proofs
Proof of Proposition (forward invariance of the GBA subspace). Equivalently ( G ( t ) , B ( T )) remains contained in the closed triangle T GBA = { ( G, B ) ∈ R : G ≥ 0 , B ≥ 0 , G + B ≤ 1 } at all times. This is easily shown by observing that the vector field ( G, [˙] B [˙] ) points towards the interior of T 1 everywhere on ∂T 1 :
 G = 0 ⇒ G [˙] = γ ≥ 0
B = 0 ⇒ B [˙] = Φ G ( G ) G ≥ 0

G = 0 ⇒ G [˙] = γ ≥ 0 B = 0 ⇒ B [˙] = Φ G ( G ) G ≥ 0 G + B = 1 ⇒ G [˙] + B [˙] = −qB ≤ 0 .  
Proof of Proposition 2 (existence and uniqueness of the GBA steady state). We define for convenience
. �
F ( G ) ≡ [γq]
γ + q
1 � G
1
G [−] [β] [G]
q
Then G is a steady state if and only if F ( G ) = Φ G ( G ). Let us first consider the case where q < β G . Then in the interval, F has exactly two roots located at βq G
Furthermore, F > 0 on and F < 0 on and F < 0 on. Since Φ G > 0 on ≤∞ then any roots of F − Φ Gin must occur in the interval. Next note that the derivative of F is
β G � q
γq F [′] ( G ) = γ + q
G
q [−] G
G
�
q/β G . Furthermore, F [′] < 0 on. Since
⊂ it follows that F is strictly monotonically decreasing on. ⊂ it follows that F is strictly monotonically decreasing on. Next since Φ G is monotonically increasing on then F− Φ G is strictly monotonically decreasing on. Since F − Φ G is continuous on R [+] and lim G→∞ ( F − Φ G )( G ) = + ∞ while ( F − Φ G )( q/β G ) = − Φ G ( q/β G ) < 0 then F − Φ G must have exactly one root in. It follows that there is a unique steady state G ∈. In the case where q ≥ β G then F has exactly one root at 1 in the interval and F > 0 on. By the same argument as before F [′] < 0 on so F is strictly monotonically decreasing on. The rest of the argument is the same as the q < β G case with q/β G replaced by 1.
Thus, a unique GBA steady state always exists for any choice of parameters for the system. See Fig. for a graphical illustration. Figure 13: Plots of F and Φ Gin the q < β G case (left) and the q ≥ β G case (right)
Graphically, we note that Gincreases as γγq + q [increases. Noting that] [ ∂] [γ] [(] γ [γq] +
Graphically, we note that Gincreases as γγq + q [increases. Noting that] [ ∂] [γ] [(] γ [γq] + q [)] [ >]
∂ q ( [γq] [)] [ >] [ γ, q >] [ γ] [ or] [ q] [ will increase the proportion]
∂ q ( γ [γq] + q [)] [ >] [ γ, q >] [ γ] [ or] [ q] [ will increase the proportion]
of grass in the steady state. This is logical since increasing γ allows faster regrowth of grass from ash, and increasing q allows faster fire quenching into ash from which grass can regrow. Similarly, increasing βq G [also increases] [ G] [. This also follows intuition: the case of increasing] q has already been discussed, while decreasing β G reduces the rate at which grass burns due to fire spread. We also note that increasing g 0 or g 1 reduces Gas expected since doing so increases the flammability of grass. Lastly, we note that the presence of the GBA steady state depends on the non-decreasing non-tonicity of Φ G . Thu,s if Φ G were constant, the result of the system containing a single unique GBA steady state would still hold. Proof of Proposition 4 (forward invariance for the nonspatial FGBA model). Analogously to the GBA system, we check for flow invariance of the tetrahedron T FGBA = { ( F, G, B ) | F ≥ 0 , G ≥ 0 , B ≥ 0 , F + G + B ≤ 1 } to ensure that any trajectory starting at an ecologically relevant condition remains ecologically relevant at all times.  F = 0 ⇒ F [˙] = 0   
G = 0 ⇒ G [˙] = γ + µF ≥ 0
 BB [˙] = Φ F
F = 0 ⇒ F [˙] = 0 G = 0 ⇒ G [˙] = γ + µF ≥ 0 B = 0 ⇒ B [˙] = Φ G ( G ) G + Φ F ( G ) F ≥ 0 F + G + B = 1 ⇒ F [˙] + G [˙] + B [˙] = −qB ≤ 0 ,
    
so T FGBA is invariant, as desired. Appendix B Computation of equilibrium burning cover value
We now examine the function B ( F ), which generates an equilibrium burning cover value for any input forest value. To find B ( F ), set the four equations in system to 0, forming an essentially three-dimensional system. From the equation for A [˙] we obtain
γ + φF + q G ( F, B ) = 1 − F − � γ + φF
Plugging this into the equation for B [˙] then gives
B. �
where
0 = c 0 ( F ) + c 1 ( F ) B + c 2 ( F ) B
c 0 ( F ) = g 0 + f 0 F
γ + φF + q c 1 ( F ) = β G − g 0 � γ + φF
+ β FF − q �
γ + φF + q c 2 ( F ) = −β G � γ + φF
�
where we work in the assumption that F ≈ 0. A simple application of the quadratic formula then gives c ~~~~
~~�~~ 1 [(] [F] [)] [ −] [c] [(] [F] [)] [c] [(] [F] [)] B ( F ) = [−][c] [(] [F] [)] [ ±] . 2 c 2 ( F )
Next note that as g 0 → 0 then
B = [(] [β] [G] [ −] [q] [)] [ ∓] [(] [β] [G] [ −] [q] [)]
2 β G ( [γ] [+] [q]
γ [q] [)]
where we assume that β G > q . Comparing this expression to the GBA steady state with g 0 → 0 in Eq. then indicates that we must choose the + sign. Next, using Mathematica we can expand B ( F ) in powers of F :
B ( F ) = d 0 + d 1 F + d 2 F + O
where
q d 0 = [(] [β] [G] [ −] [q] [)] [γ]
β G ( q + γ ) [+] � β G ( β G − q )
g 0 + O] �
d 1 = [(] [β] [f] [ −] [β] [G] [)] [γ] [(] [q] [ +] [γ] [)][ +] [ q] [(] [β] [G] [ −] [q] [)] [φ] + 1
β G ( q + γ ) � β g − q
+ O]
( q − β F )( β F − β G ) f 0 + � � ( q − β G )
g 0 �
[φ] [)][ + ][(] [q][ −] [β] [G] [)] [φ] [)] d 2 = [qφ] [((] [β] [F] [ −] [β] [G] [) . β G ( q + γ ) Using our orders of magnitude approximation β G ≈ q ≫ γ ≫ φ ≈ 1 and β G ≫ β F we find that
[γ]
q [+] q
[γ] q q q [,]
d 0 ≈ [γ]
d 1 ≈ [γ] q [g] [,] q d 2 ≈ [φ] q [] [,] q so d 0 ≈ d 1 ≫ d 2, and we can reasonably approximate
B ( F ) ≈ d 0 + d 1 F,
and it is also apparent that
~~′~~ B = d 0 , B = d 1 .
We can perform a similar analysis for F ≈ 1. In this case, we write
B ( F ) = d [˜] 0 + d [˜] 1 + O )
where
˜ f 0 d 0 = + O] [,] q − β F
˜ q − β Gd 1 = � ( q − β F )
1 f 1 + � � β f − q g 1 + O] [.] �
Appendix C Simulation Testing
Below, we perform some sense checks on the model to ensure that the simulation codes are bug-free and behaves as expected. We test each transition rate separately by introducing artificial situations with readily predictable behaviourand comparing the simulation results with analytic predictions. C.1 Testing non-spatial transitions
The simplest transitions to test are the spontaneous, non-spatial transitions governed by the parameters q ( B → A ), γ ( A → G ) and µ ( F → G ) corresponding to fire quenching, grass regrowth, and tree mortality, respectively. All testing was performed with L = 100 and N = 500. To test the fire quench rate we began the simulation with all states in the burning state and removed grass regrowth by setting γ = 0. The times at which each burning state transitioned into the ash state were recorded at various values of q . The quench times are expected to follow a Poisson distribution with rate q, giving an expected quench time of q with variance q [.] Similarly, to check the grass regrowth rate from ash we began with all sites in the ash state. We set g 0 = g 1 = 0 to remove spontaneous grass fires. The times at wich each ash site transitoned into a grass state was recorded for various values of γ . As expected, he regrowth times followed a Poisson distribution with rate γ and expected regrowth time
γ with variance
γ [.] Figure 14: Testing the fire quench (left), grass regrowth (center), and tree mortality (right) transitions
Lastly, to check forest mortality transitions we began the simulation with all sites in the forest state and removed forest spread by setting φ G = φ A = 0. Grass and forest fires were also eliminated by setting g 0 = g 1 = f 0 = f 1 = 0. The times at which each forest site transitioned into a grass state was recorded at various values of µ . The tree lifetimes are expected to follow a Poisson distribution wtih rate µ giving an expected lifetime of µ [with]
variance
µ [.]
C.2 Testing spatial transitions
We next tested the four directly spatial transitions in the FGBA model: the spread of forest through grass and through ash and the spread of fire through grass and through forest. Fire spread through grass was tested introducing the initial condition where a single randomly chosen vegetation site in the domain would be in a grass state while the remaining sites are all randomly assigned to be a either burning or ash states according to some fixed probability p ash of the state being ash. Further, g 0 = g 1 = γ = q = 0 was set to remove spontaneous grass ignitions and grass regrowth and prevent fire quenching. The burn rate at the grass site is then expected to be approximately β G p ash . The time at which the grass site transitioned to a burning state was then measured at various values of β G and p ash . The average of 600 trials for each pair of values ( β G , p ash ) is plotted in Fig.. Fire spread through forest was tested in an analogous situation where instead of a single grass site, a single forest site was used and f 0 = f 1 = µ = 0 was set to eliminate all possible transitions except the forest to burning transition. The burn rate at the forest site is then expected to be approximately β F. The time at which the grass site transitioned to a burning state was then measured at various values of β F and p ash . The average of 600 trials for each pair of values ( β F , p ash ) is plotted in Fig..
Next, the spread of forest through grass was tested by randomly choosing a single vegetation site to be grass; the remaining sites are randomly chosen to be in the forest state with probability p forest and are otherwise in the ash state. To prevent all other possible transitions besides the grass to forest transition we set g 0 = g 1 = f 0 = f 1 = q = γ = µ = φ A = 0. Then the rate at which the grass site should transition to forest is approximately φ G p forest . Figure 15: Testing the spatial transition of grass burning (left) and forest burning (right) with N = 2000. The time at which the grass site transitioned to a forest state was then measured at various values of φ G and p forest . The average of five trials for each pair of values ( φ G , p forest ) is plotted in Fig.. Lastly, the spread of forest through ash was tested by randomly choosing a single vegetation site to be ash while the remaining sites are randomly chosen to be in the forest state with probability p forest and are otherwise in the grass state. To prevent all other possible transitions besides the ash to forest transition, we set g 0 = g 1 = f 0 = f 1 = q = γ = µ = φ G = 0. Then the rate at which the ash site should transition to forest is approximately φ A p forest . The time at which the grass site transitioned to a forest state was then measured at various values of φ A and p forest . The average of five trials for each pair of values ( φ A , p forest ) is plotted in Fig.. Appendix D Computation of linear regimes in non-fire forest mortality bifurcation plot
Note the largely linear regimes within the non-fire forest mortality ( µ ) bifurcation plots aside from an elbow at F ≈ 0 . 6.
To determine exact expressions for these two regimes, we note that the equation for F [˙] in system can only vanish for F > 0 if
φ − ( φ + β F )( B ( F )) − Φ F ( F ) − µ = 0
For F < 0 . 6 we approximate B ( F ) as a linear function of F (see App. B) so that we obtain a line with x-intercept
( β G − q ) γ µ = φ − ( φ + β F ) � β G ( q + γ )
29
− f 0 Figure 16: Testing of the spatial transition of forest spread into grass (left) and ash (right). and slope
β G ( q + γ ) m =
β F ( β F − β G ) γ ( q + γ ) + qβ F ( β G − q ) φ [.]
Similarly, for F > 0 . 6 an analogous calculation gives a line with a y -intercept at and slope
[β] [F] [)] [ −] [f] F = [φ][ −] [(] [φ] [ +]
φ + ( φ + β F ) d [˜] 1
1
m = .
−φ − ( φ + β F ) d [˜] 1
~~′~~ To understand the origin of the elbow at F ≈ 0 . 6 we note that B ( F ) is undefined when the discriminant c 1 ( F ) − 4 c 0 ( F ) c 2 ( F ) vanishes. To zeroth order in f 0 and g 0 we then find ~~′~~ that B ( F ) is undefined when
F [∗] = [q][ −] [β] [G] . β G − β F
Furthermore, to zeroth order in f 0 and g 0
~~′~~ B ( F ) =
0 F > F [∗]
−q φ + qβ G+( β F −β G )( γ + F φ ) + qβ F
F < F [∗]
� β G ( q + γ + Fφ )
−q φ + qβ G+( β F −β G )( γ + F φ ) + qβ F
F − G FF < F [∗] . β G ( q + γ + Fφ )
~~′~~ Note that all F -dependent terms are suppressed since q ≫ γ ≫ φ so B ( F ) is approximately constant in its two regimes, so B ( F ) is approximately linear.