Mathematical analysis of a size structured tree-grass competition model for savanna ecosystems

—Several continuous-time tree-grass competition models have been developed to study conditions of long-lasting coexistence of trees and grass in savanna ecosystems according to environmental parameters such as climate or ﬁre regime. In those models, ﬁre intensity is a ﬁxed parameter while the relationship between woody plant size and ﬁre-sensitivity is not systematically considered. In this paper, we propose a mathematical model for the tree-grass interaction that takes into account both ﬁre intensity and size-dependent sensitivity. The ﬁre intensity is modeled by an increasing function of grass biomass and ﬁre return time is a function of climate. We carry out a qualitative analysis that highlights ecological thresholds that summarize the dynamics of the system. Finally, we develop a non-standard numerical scheme and present some simulations to illustrate our analytical results.


I. Introduction
Savannas are tropical ecosystems characterized by the durable co-occurrence of trees and grasses (Scholes 2003, Sankaran et al. 2005) that have been the focus of researches since many years. Savanna-like vegetations cover extensive areas, especially in Africa and understanding savannas history and dynamics is important both to understand the contribution of those areas to biosphereclimate interactions and to sustainably manage the natural resources provided by savanna ecosystems. At biome scale, vegetation cover is known to display complex interactions with climate that often feature delays and feed-backs. For instance any shift from savanna to forest vegetation not only means increase in vegetation biomass and carbon sequestration but also may translate into changes in the regional patterns of rainfall (Scheffer et al. 2003, Bond et al. 2005). In the face of the ongoing global change, it is therefore important to understand how climate along with local factors drive the dynamics of savannas ecosystems. In many temperate and humid tropical biomes, forest vegetation in known to recover quickly from disturbances and woody species are expected to take over herbaceous species. Yet in the dry tropics, it is well-known that grassy and woody species may coexist over decades although their relative proportion may show strong variations (Scholes 2003, Sankaran et al. 2005, 2008. Savanna-like ecosystems are diverse and explanations found in the literature about the longlasting coexistence of woody and grassy vegetation components therefore relate to diverse factors and processes depending on the location and the ecological context. Several studies have pointed towards the role of stable ecological factors in shaping the tree to grass ratio along large-scale gradients of rainfall or soil fertility (Sankaran et al. 2005(Sankaran et al. , 2008. Other studies have rather emphasized the reaction of vegetation to recurrent disturbances such as herbivory or fire (Langevelde et al. 2003, D'Odorico et al. 2006, Sankaran et al. 2008, Smit et al. 2010, Favier et al. 2012 and references therein). Those two points of view are not mutually-exclusive since both environmental control and disturbances may co-occur in a given area, although their relative importance generally varies among ecosystems. Bond et al. (2003) proposed the name of climate-dependent for ecosystems that are highly dependent on climatic conditions (rainfall, soil moisture) and firedependent or herbivore-dependent for ecosystems which evolution are strongly dependent on fires or herbivores. In a synthesis gathering data from 854 sites across Africa, Sankaran et al. (2005) showed that the maximal observed woody cover appears as water-controlled in arid to semi-arid sites since it directly increase with mean annual precipitation (MAP) while it shows no obvious dependence on rainfall in wetter locations, say above c. 650 mm MAP where it is probably controlled by disturbance regimes. Above this threshold, fire, grazing and browsing are therefore required to prevent tree canopy closure and allow the coexistence of trees and grasses.
Several models using a system of ordinary differential equations (ODES) have been proposed to depict and understand the dynamics of woody and herbaceous components in savanna-like vegetation. A first attempt (Walker et al. 1981) was orientated towards semiarid savannas and analyzed the effect of herbivory and drought on the balance between woody and herbaceous biomass. This model refers to ecosystems immune to fire due to insufficient annual rainfall. Indeed, fires in savanna-like ecosystems mostly rely on herbaceous biomass that has dried up during the dry season. As long as rainfall is sufficient, fire can thus indirectly increase the inhibition of grass on tree establishment in a way far more pervasive than the direct competition between grass tufts and woody seedlings.
More recently, several attempts have been made (see Accatino et al. 2010, De Michele et al. 2011 and references therein) to model the dynamics of fire-prone savannas on the basis of the initial framework of Tilman (1994) that used coupled ODES to model the competitive interactions between two kinds of plants. On analogous grounds, Langevelde et al. (2003) have developed a model taking into account fires, browsers, grazers and Walter's (1971) hypothesis of niche separation by rooting zone depth. Models relying on stochastic differential equations have also been used (Baudena et al. 2010). Notably, Accatino et al. (2010) and De Michele et al. (2011) focused on the domain of stability of tree-grass coexistence with respect to influencing "biophysical" variables (climate, herbivory). However, fire was considered as a forcing factor independent of climate and vegetation, while woody cover was treated as a single variable with no distinction between seedling/saplings which are highly fire sensitive and mature trees which are largely immune to fire damages. The way in which the fundamental, indirect retroaction of grass onto tree dynamics is modeled is therefore to be questioned. In the present paper we therefore a model that differs in this respect.
Thus, to take into account the role of fire in savanna dynamics, we consider a tree-grass compartmental model with one compartment for grass and two for trees, namely fire-sensitive individuals (like seedlings, saplings, shrubs) and non-sensitive mature trees. Based on field observations and experiments reported by Scholes and Archer (1997) and by Scholes (2003), we develop a system of three coupled non-linear ordinary differential equations (ODES), one equation per vegetation compartment that describes savanna dynamics. In addition, we model fire intensity (i.e. impact on sensitive woody plants) as an increasing function of grass biomass. Compared to existing models, our model aims to properly acknowledge two major phenomena, namely the fire-mediated negative feedback of grasses onto sensitive trees and the negative feed-back of grown-up, fire insensitive trees on grasses. We therefore explicitly model the asymmetric nature of tree-grass competitive interactions in fire-prone savannas.
After some theoretical results of the continuous fire model, though which we highlighted some ecological thresholds that summarize savanna dynamics and some interesting bistability, we present an appropriate non-standard numerical scheme (see Anguelov et al. 2012and Dumont et al. 2010 for the model considered and we end with numerical simulations. We show that the fire frequency and the competition parameters are bifurcation parameters which allow the continuous fire model of asymmetric tree-grass competition to converge to different steady states.

II. THE CONTINUOUS FIRE MODEL OF ASYMMETRIC TREE-GRASS COMPETITION (COFAC)
As we have mentioned before, we consider the class of sensitive tree biomass (T S ), the class of non-sensitive tree biomass (T NS ) and the class of grass biomass (G). We model the fire intensity by an increasing function of grass biomass w(G).
To built up our model, we consider the following assumptions.
1) The grass vs. sensitive-tree competition has a negative feedback on sensitive tree dynamics.
2) The grass vs. non sensitive-tree competition has a negative feedback on grass dynamics. 3) After an average time expressed in years, the sensitive tree biomass becomes non sensitive to fire. 4) Fire only impacts grass and sensitive Tree. We also consider the following parameters.
• There exists a carrying capacity K T for tree biomass (in tons per hectare, t.ha −1 ). • There exists a carrying capacity K G for grass biomass (in tons per hectare, t.ha −1 ). • Sensitive tree biomass is made up from non sensitive tree biomass with the rate γ NS (in yr −1 ) and from existing sensitive tree biomass with the rate γ S (in yr −1 ). • Sensitive tree biomass has a natural death rate µ S (in yr −1 ). • Non sensitive tree biomass has a natural death rate µ NS (in yr −1 ). • f is the fire frequency (in yr −1 ). • Grass biomass has a natural death rate µ G (in yr −1 ).
• 1 ω S is the average time, expressed in year, that a sensitive tree takes to become non sensitive to fire.
is the average time that a tree spends in the sensitive tree class without competition and fires. • σ G is the competition rate, for light or/and nutrients, between sensitive tree and grass (in ha.t −1 .yr −1 ). • σ NS is the competition rate, for light or/and nutrients, between non sensitive tree and grass (in ha.t −1 .yr −1 ). • η S is the proportion of sensitive tree biomass that is consumed by fire. • η G is the proportion of grass biomass that is consumed by fire.
Remark 1. Competition parameters σ G and σ NS are asymmetric, indeed σ G inhibits sensitive tree (T S ) growth and there is no reciprocal inhibition; likewise, σ NS inhibits grass (G) growth.
Based on these ecological premises, and taking into account the effect of fire as a forcing continuous in time, which is the classical approach, we propose a model for the savanna vegetation dynamics through a system of three interrelate non-linear equations. The COFAC is given by (2) For this continuous fire model, the fire intensity function w is chosen as a sigmoidal function of grass biomass because we want first to investigate the ecological consequences of the nonlinear response of fire intensity to grass biomass, while nearly all published models using differential equations so far assumed a linear response. Non linearity is justified since whenever grass biomass is low fires are virtually absent while fire impact increases rapidly with grass biomass before reaching saturation. Thus, where G 0 = g 2 0 is the value of grass biomass at which fire intensity reaches its half saturation (g 0 in tons per hectare, t.ha −1 ).
The feasible region for system (1) is the set Ω defined by

A. Existence of equilibria, ecological thresholds and stability analysis
We set

1) Existence of equilibria:
Setting the right hand-side of system (1) to zero, straightforward computations lead to the following proposition is ecologically meaningful when R 0 2 > 1 The point E G when it exists is the grassland equilibrium.
• The savanna equilibrium point E TG = (T * S , T * NS , G * ), with T * S , T * NS and G * given in Appendix A, has an ecological significance whenever Remark 2. The number of savanna equilibria depends on the form of the function w.
• If w(G) = G, then the COFAC has at most one savanna equilibrium.  Holling type III function), then the COFAC has at most three savanna equilibria.

2) Ecological thresholds interpretation:
The qualitative behaviors of the COFAC depend on the following thresholds where • R 0 1 is the sum of the average amount of biomass produced by a sensitive/young plant, without fires and competition with grass, and the average amount of biomass produced by a mature plant multiplied by the proportion of young plants which reach the mature stage.
• R G 1 is the sum of the average amount of biomass produced by a sensitive/young plant, in presence of fires and competition with grass, and the average amount of biomass produced by a mature plant multiplied by the proportion of young plants which reach the mature stage.
• R 0 2 is the average amount of biomass produced per unit of grass biomass during its whole lifespan in presence of fires and and free from competition with non-sensitive trees.
is the average biomass produced per unit of grass biomass during its whole lifespan in presence of fires and experiencing competition from non-sensitive trees.

3) Stability analysis: Let
We have the following result: Theorem 1. If R 0 1 < 1 and R 0 2 < 1, then the desert equilibrium E 0 is globally asymptotically stable.
Proof: See Appendix B.
Theorem 2. If R 0 1 > 1, then the forest equilibrium E T exists.
, then the forest equilibrium E T is globally asymptotically stable.
Proof: See Appendix C. Furthermore, using the same approach as in the proof of Theorem 2, we derive the following results > 1, RḠ 1 < 1 and R < 1, then the grassland equilibrium E G is globally asymptotically stable.
Theorem 4. Suppose that R 0 1 > 1, R 0 2 > 1 and R > 1. We have the following three cases: • The savanna equilibrium E TG is locally asymptotically stable (LAS) when it is unique. • When there exists two savanna equilibria, one is LAS and the other is unstable. • When there exists three savanna equilibria, two are LAS and one is unstable. Thus System (1) will converges to one of the two stable savanna equilibria depending on initial conditions.
Proof: See Appendix D.

4) Summary table of the qualitative analysis:
The qualitative behavior of system (1) is summarized in the following In Table I, the notations U, L and G stand for unstable, locally asymptotically stable, globally asymptotically stable, respectively, while the notation L means that we have the global stability if there are no periodic solutions.  Table I are interesting because in these cases, one ton of grass biomass will produce during it lifespan at least one ton of grass biomass (R 0 2 > 1) and simultaneously, one ton of tree biomass (sensitive and non sensitive) will produce during it lifespan at least one ton of tree biomass (R 0 1 > 1). Moreover, it is also in these cases that we have the most interesting situations of savanna dynamics, namely bistability cases (Lines 2, 4, 7 in Table 1) and a tristability case (Line 6 in Table  1).

A. A nonstandard scheme for the COFAC
System (1) is discretized as follows: can be written in the following matrix form: Using (5), the numerical scheme (4) can be rewritten as follows: In particular, choosing φ such that lead to positive diagonal terms and nonpositive off diagonal terms. We need to show that B(X k ) is invertible. Obviously 1 − φ(h)A k 22 is a positive eigenvalue. Let us define N, a submatrix of matrix B(X k ), as follows We already have trace(N) > 0. Then, a direct computation shows that det Thus we have α(N) > 0, i.e. the eigenvalues have positive real parts, which implies that B(X k ) is invertible. Finally, choosing with matrix B(X k ) is an M-matrix. Furthermore, assuming X k ≥ 0, we deduce Lemma 1. Using the expression of φ defined in (8), the numerical scheme (4) is positively stable ( i.e for X k ≥ 0, we obtain X k+1 ≥ 0).
An equilibrium X e of the continuous model (1) verifies A(X e )X e = 0. Multiplying the above expression by φ(h) and summing with X e yields (Id 3 − φ(h)A(X e ))X e = X e , Thus, we deduce that the numerical scheme (4) and the continuous model (1) have the same equilibria which are (assumed to be) hyperbolic.
The dynamics of model (1) can be captured by any number Q satisfying where λ ∈ sp(J) with J i j = ∂A i ∂X j . We also have the following result:

B. NUMERICAL SIMULATIONS AND BIFUR-CATION PARAMETERS
In literature we found the following parameters values We now provide some numerical simulations to illustrate the theoretical results and for discussions.
We choose With the chosen parameters, the savanna equilibrium is stable, i.e. sensitive trees, nonsensitive trees and grasses coexist. Figure  1 presents the 3D plot of the trajectories of system (1). It illustrates that the savanna equilibrium point is stable. Figure 1 also illustrates the monostability situation presented in Ligne 1 of Table I. • Bistability -Bistability involving forest and grassland equilibria. The state trajectories of the model will converge to a state depending of initial quantity. We choose The 3D plot of the trajectories of system (1) is depicted in Figure 2. It clearly appears that the forest and grassland equilibria are stables. Figure 2 illustrates the bistability situation presented in Ligne 7 of Table I. -Bistability involving forest and savanna equilibria. The state trajectories of the model will converge to a state depending on initial quantity. We choose For these parameters, there exist two savanna equilibria but only one is stable as shown in Figure 3. Figure 3 also illustrates the bistability situation presented in Ligne 2 of Table I  Remark 5. Note also that we didn't observed periodic behaviors in the previous simulations, considering the set of parameters presented in Table 2, while their existence cannot be completely ruled out by the analytical analysis.

C. Some bifurcation parameters
In this section we emphasize on some bifurcation parameters of system (1) which are such that the COFAC can converge to different steady state depending on the variation of these parameters.
• The grass vs. sensitive-tree competition parameter σ G is a bifurcation parameter. Figure  4 presents how the system (1) changes from the savanna state to the grassland state as a function of the grass vs. sensitive-tree competition parameter σ G . We choose For these parameters values, system (1) undergoes a transcritical bifurcation. Indeed, we move from ligne 1 to ligne 5 of Table I Fig. 4. From savanna to grassland as a function of the grass vs. sensitive-tree competition parameter σ G . From left to right, the fire period τ = 1 f is fixed, while the grass vs. sensitive-tree competition parameter σ G increases. In (a) (τ = 5, σ G = 0), in (b) (τ = 5, σ G = 0.02) and in (c) (τ = 5, σ G = 0.04) • The fire period parameter τ = 1 f is a bifurcation parameter. Figure 5 presents a shift of the convergence of system (1) from the forest state to the grassland state as a function of the fire period τ. We choose For these parameters values, system (1) undergoes a forward bifurcation. Indeed, we move from ligne 3 to ligne 7 of Table I Suppose now that fire period is fixed and the grass vs. sensitive-tree competition parameter σ G varies. Figure 6 illustrates a shift of the convergence of system (1) from the forest state to the grassland state through a savanna state as a function of the grass vs. sensitivetree competition parameter σ G . For the parameters values in Table III, system (1) exhibits to bifurcation phenomena: a pitchfork bifurcation and a transcritical bifurcation. We move from ligne 3 (figure 6 (a), -(b)) to ligne 7 (figure 6 (e), -(f)) through ligne 2 (figure 6 (c), -(d)) of Table I. Indeed, in figure 6 (a), -(b) E TG is undefined, E G is unstable, E T is stable. In in figure 6 (c), -(d), E G remains unstable but we have bistability between E TG and E T : it is a case of pitchfork bifurcation. In figure 6 (d), -(f), E TG becomes unstable and we have a bistability between E T and E G : it is a case of transcritical bifurcation.
and R are given in Appendix F.  Fig. 6. From forest to grassland, with a transition through a savanna state, as a function of the sensitive tree-grass competition parameter. From left to right, the fire period τ is fixed at 4, while the grass vs. sensitive-tree competition parameter σ G increased.
• The grass vs. non sensitive-tree competition parameter σ NS is a bifurcation parameter. A shift of the convergence of system (1) from the grassland state to the forest state as a function of the grass vs. non sensitive-tree competition parameter σ NS is depicted in Fig.  7.
We choose For these parameters values, system (1) exhibits a pitchfork bifurcation. We move from ligne 5 ( figure 7 (a), -(b)) to ligne 7 (figure 7 (c), -(d)) of Table I. Indeed, in figure 7 (a), -(b) E TG is undefined, E T is unstable, E G is stable. In figure 7 (c), -(d), E TG exits but it is unstable and we have bistability between E G and E T : it is a case of pitchfork bifurcation. Values of R 0 1 , R 0 2 , R G 1 , R T NS 2 and R are given in Appendix G.  Fig. 7. From grassland to forest as a function of the grass vs. non sensitive-tree competition parameter σ NS . From left to right, the fire period τ is fixed, while the grass vs. non sensitive-tree competition parameter σ NS increases.

V. Conclusion and discussion
In this work, we present and analyze a new mathematical model to study the interaction of tree and grass that explicitly makes fire intensity dependent on the grass biomass and distinguishes two levels of fire sensitivity within the woody biomass (implicitly relating to plant size and bark thickness). Fire was considered as a timecontinuous forcing as in several existing models (Langevelde et al. 2003, Accatino et al. 2010, De Michele et al. 2011 and reference therein) with a constant frequency of fire return that can be interpreted as mainly expressing an external forcing to the tree-grass system from climate and human practices. What is novel in our model is that fire impact on tree biomass is modeled as a non-linear function w of the grass biomass. Using a non-linear function is to our knowledge only found in . But this latter model made peculiar assumptions and does not predict grassland and forest as possible equilibria (only desert and savanna). The advantage of a non-linear function is that it can account for the absence of fire at low biomass. As a consequence and Distinguishing fire sensitive vs. fire insensitive woody biomass lead to three variables expressing fractions of the above ground phytomass, namely grass and both fire-sensitive andinsensitive woody vegetation. It featured three coupled, non-linear ordinary differential equations.As several existing models (Baudena et al. 2010, our model acknowledges two major phenomena that regulate savanna dynamics, namely the fire-mediated negative feedback of grasses onto sensitive trees and the negative feedback of grown-up, fire insensitive trees on grasses. We therefore explicitly model the asymmetric nature of tree-grass competitive interactions in fireprone savannas. The analytical study of the model reveals three possible equilibria excluding tree-grass coexistence (desert, grassland, forest) along with equilibria for which woody and grassy components show durable coexistence (i.e. savanna vegetation). The number of such equilibrium points depends on the function used to model the increase of fire intensity with grass biomass(see Remark 2); for our model, we can have at most three savanna equilibria. We identified four ecologically meaningful thresholds that defined in parameter space regions of monostability, bistability as in Accatino et al. 2010, De Michele et al. 2011 and tristability with respect to the equilibria. Tristability of equilibria may mean that shifts from one stable state to another may often be less spectacular that hypothesized from previous models and that scenarios of vegetation changes may be more complex.
The model features some parameters that have been analytically identified as liable to trigger bifurcations (i.e., the state variables of the model converges to different steady states), notably parameters σ NS and σ G of asymmetric competition that embody the depressing influence of fire insensitive trees on grasses and of grasses on sensitive woody biomass respectively. Since tree-grass asymmetric competition is largely mediated by fire, this finding of the role of those two parameters is not intuitive and is the result of the modeling effort and of the analytical analysis. Since such parameters that quantify direct interactions between woody and grassy components appear crucial to understand the tree-grass dynamics in savanna ecosystems and for enhanced parameter assessment, they could be the focus of straightforward field experiments that would not request manipulating fire regime. Another bifurcation parameter is the fire frequency, f , (or fire period parameter τ = 1 f ) which has been assumed to be an external forcing parameter that integrates both climatic and human influences. Frequent fires preclude tree-grass coexistence and turn savannas into grasslands. In the wettest situations, or under subequatorial climates, very high fire frequencies (above one fire per year) seem to be needed to prevent the progression of forests over savannas (unpublished data of experiments carry out at La Lopé National Park in Gabon).
However, it is questionable to model fire as a continuous forcing that regularly removes fractions of fire sensitive biomass. Indeed, several months can past between two successive fires, such that fire may be considered as an instantaneous perturbation of the savanna ecosystem. Several recent papers have proposed to model fires as stochastic events while keeping the continuoustime differential equation framework (Beckage et al. 2011) or using time discrete matrix models (Accatino & De Michele 2013). But in all those examples, fire characteristics remain mainly a linear function of grass biomass. Another framework that we will explore in a forthcoming work in order to acknowledge the discrete nature of fire events is based on system of impulsive differential equations (Lakshmikantham et al. 1989, Bainov andSimeonov 1993).
Appendix A: Expressions of T * S , T * NS and G * After straightforward but long computation, we show that where G * is solution of with We summarize the problem of existence of solutions of equation (11) in the following Table   TABLE IV Existence of solutions of equation (11) A B Number of solutions > 0 > 0 0 or 2 solutions < 0 1 or 3 solutions < 0 > 0 1 solution < 0 0 solution Note that solutions G * of (11) that give rise to savanna equilibria must satisfy 0 < G * < Appendix B: Proof of Theorem 1 . In a matrical writing, System (1) reads as is a Metzler matrix ( i.e all its offdiagonal terms are nonnegative) and α(A max (X)) ≤ 0 if α(A) ≤ 0 and α(D) ≤ 0 where α denotes the stability modulus. Moreover, for matrix D, Furthermore, Thus, from relations (13), (14), (15) and (16) we deduce that the desert equilibrium (0; 0; 0) is globally asymptotically stable whenever R 0 1 < 1 and R 0 2 < 1.
Appendix C: Proof of Theorem 2 If R 0 1 > 1, then the forest equilibrium E T exists. • Using the Jacobian matrix of system (1) at E T , one can prove that E T is locally asymptotically stable if R T NS 2 < 1. • The solution G of system (1) verify So, if R 0 2 < 1, then lim t→+∞ G(t) = 0.
Moreover, the solutions T S and T NS of system (1) admit as a limit system, the system: S . Then, one has and by the Bendixson-Dulac theorem, we deduce that system (19) don't admits a periodic solution in Ω 2 . Moreover, the equilibrium (T S ,T NS ) exists if R 0 1 > 1 and using the Jacobian matrix of system (19), we deduce that (T S ,T NS ) is locally asymptotically stable and then, globally asymptotically stable since there is no periodic solution. Thus, if R 0 2 < 1, then one has lim (20) From the Jacobian matrix of system (1), the equilibria (0, 0,Ḡ) and (T S , T NS , G ) are unstable if RḠ 2 > 1 and R < 1.
In the sequel, we suppose that RT NS 2 < 1 to process with the discussion.
Thus, C 2 > 0 if and only if R > 1.
Finally, we deduce that the savanna equilibrium E TG , when it is unique, is locally asymptotically stable if R = R(G) > 1. The first part of Theorem 4 holds.
One should note that C 2 > 0 means that the slope of w (the sigmoidal function) is less than the slope of F where F is given by (11).
Furthermore, by using relation (11) we deduce part 2 and part 3 of Theorem 4 graphically as follow Fig. 8. There exist two savanna equilibria but one is stable and the other is unstable. Fig. 9. There exist three savanna equilibria two are stable and one is unstable. Thus system (1) will converge to one of the two stable equilibria depending on initial conditions.