Tropical tree cover in a heterogeneous environment: A reaction-diffusion model

Observed bimodal tree cover distributions at particular environmental conditions and theoretical models indicate that some areas in the tropics can be in either of the alternative stable vegetation states forest or savanna. However, when including spatial interaction in nonspatial differential equation models of a bistable quantity, only the state with the lowest potential energy remains stable. Our recent reaction-diffusion model of Amazonian tree cover confirmed this and was able to reproduce the observed spatial distribution of forest versus savanna satisfactorily when forced by heterogeneous environmental and anthropogenic variables, even though bistability was underestimated. These conclusions were solely based on simulation results for one set of parameters. Here, we perform an analytical and numerical analysis of the model. We derive the Maxwell point (MP) of the homogeneous reaction-diffusion equation without savanna trees as a function of rainfall and human impact and show that the front between forest and nonforest settles at this point as long as savanna tree cover near the front remains sufficiently low. For parameters resulting in higher savanna tree cover near the front, we also find irregular forest-savanna cycles and woodland-savanna bistability, which can both explain the remaining observed bimodality.


Model construction
General overview 2 The model in space and time can be written as a system of stochastic partial differential equations of the reaction-diffusion type. In ecological context, one describes the 4 dynamics of a species density as the sum of a reaction term, representing local demography, and a diffusion term, representing migration of the species through space 6 from areas with high to areas with low density. Our model can be most compactly expressed in vector form where Y = (S, T, F, G) with S representing savanna sapling cover density, T savanna tree cover density, F forest tree cover density, and G grass cover density. A is a vector 10 of exogenous environmental variables or parameters that force the system, such as mean rainfall, rainfall seasonality and soils [8]. These forcing variables are in general all 12 heterogeneous in space, but in the main text we kept rainfall seasonality and soils constant, and took mean annual rainfall as the only heterogeneous forcing. J is a vector 14 of reaction terms representing local population dynamics involving both gains and losses and contains nonlinear terms in both Y and A. Note that one other forcing variable is 16 contained in J, being human impact. With regards to qualitative dynamics and steady state distributions, the particular choices of the functional forms of J are arbitrary to 18 some extent. As long as we choose the right shape, the phase portrait should be topologically equivalent to the true functional form. D is a diagonal matrix with 20 diffusion constants. We take the forcing variables as constant in time. This is done by replacing A(x, t) by its long-term mean, which is only a function of space. We denote it 22 further as A(x) = A. We only consider 1D space here. Hence, ∇ 2 = ∂ 2 x , x = x and y = y. In 2D, front dynamics will be influenced by front curvature but this is minimal 24 for the spatial scales considered [19].
Local rates of change 26 Here, we show how the local rates of change J Y for Y any of the cover types (S, T, F, G) are chosen. As in any population model, we have 28 change = gain − loss Below, the set of gain processes P G contains recruitment and growth, while that of the loss processes P L contains mortality from competition for resources, drought, fire and 30 human impact. Each of those processes can be captured with a different term, such that the equations of the cover types S, T, F and G take the form As external climatic/edaphic forcing we choose where M AR stands for the observed multi-annual mean of rainfall, M SI Markham's 34 seasonality index and EF S the edaphic suitability for forest [this lets our model agree with that of [8] but using a more compact notation]. Y refers to any of the cover types. 36 Functions G Y,i and L Y,i are respectively, total gains and total losses per time of species Y by process i. The functional forms will be chosen inspired by an understanding of the 38 effect of all relevant processes.
Gain functions consist of recruitment, competition, drought, fire, humans, other They can occur due to: recruitment to an older stage (i = r), interspecific competition (i = c), drought (i = d), human impact (i = h) or other causes (i = o). Below, each of 46 the gain/loss terms are discussed.
• Cover gains and losses due to expansion. Gains due to expansion involve increase 48 of cover area of a more competitive cover type at cost of the cover area of a less competitive type. Hence, with Y the species that expands its cover and v c,Y the competitiveness vector for cover type Y , which has ones at its elements where the corresponding cover type 52 is less competitive than, and can be colonized by, Y . The rationale is that colonization by a species occurs due to an interaction of the space that is occupied 54 by that species with the space that can be colonized by it. Having chosen the competitive hierarchy F > T > S > G [10] in absence of water limitation, we have These are removed from the areas of the less competitive types. Hence, losses due to expansion of competing types are The only age-structured part of our model is that of savanna saplings S and May 31, 2019 4/15 savanna (adult) trees T. Note therefore that G S,e has T G instead of SG because 60 expansion only comes from interaction of adult species T with places that saplings S can colonize (G). The dependence of R Y on A captures the effect of water 62 availability on growth, which we choose as with r Y the maximal growth rate, k R Y as growth rate increase for every 64 component of A and a Y fixing the rate for A = 0. This function captures saturation of growth rate where water limitation is less severe, which is supported 66 by empirical work (e.g. saturation of NDVI as a function of rainfall and high temporal correlations between rainfall and NDVI below saturation; [35,36]].

68
• Gains and losses due to recruitment. Saplings will recruit into adults, such that, where recruitment rate Q is a function of burnt area fraction Φ, Q 0 is the recruitment rate in absence of fire and Q 0 (1 − h) the recruitment rate in presence of fire for a particular year, with 0 ≤ h ≤ 1. Hence, in agreement with 72 previous empirical work, fire affects the establishment rather than the mortality of savanna trees.

74
• Base mortality. This is the base mortality in absence of fire, drought, competition or human impact. As is common in ecology, we choose this mortality linear in Y , with m Y,o the base mortality rate of Y .
• Losses due to drought. Drought-related mortality rate M Y (A) will be chosen such 78 that drought effects cause increased mortality below a certain threshold of A. We or when absorbing the base mortality rate into this term, Such nonlinear increase of mortality with dryness is assumed to be a consequence of exceedance of thresholds related to tree water availability.

84
• Losses due to fire. Φ is burnt area fraction per year. Grasses and savanna tree saplings resprout very rapidly so they are assumed to be unaffected by fire on the 86 considered time scale. Savanna trees are fire adapted so also they are also assumed to be unaffected by fire [10]. Hence, in this model, only forest trees 88 experience direct mortality due to fire. This mortality is chosen proportional to burnt area and burning is assumed to occur in a homogeneously distributed feedbacks [4]. A fire feedback operating on small spatial scales is crucial for producing the two stable states, savanna and forest [7,10,37]. We parameterized further also depend on climate. The double-striked notation is used to distinguish it from the functions and variables related to the cover types. The fundamental 104 process responsible for the fire feedback is small-scale spatial fire percolation over a fire prone layer. Simulations have shown [37] that this process induces a sharp 106 increase of fire-related mortality around the percolation threshold, which occurs when about 60% of the landscape is fire-prone. We do not intend to model this Here, k c is a constant vector and Y c,0 a constant scalar. The elements of k c represent the sensitivity of Y c (A) to the different components of A. Y c has a value 120 of about 40% for common conditions, which is the tree cover value at which fire has been observed to increase [8,38,39].

122
• Losses due to human impact. Deforestation of forest trees is chosen as where C(z) is the deforestation rate and z distance from anthropogenically 124 May 31, 2019 7/15 impacted areas. We choose such that the deforestation rate decays with distance from impacted areas. c is 126 the the maximum deforestation rate, which occurs in agricultural areas (z = 0).
• Gains due to mortality of other cover types. When any cover type loses space, it 128 makes place for other cover types. When this does not occur due to competition or recruitment, grass is the default cover type that gains ground. This agrees with 130 taking the assumption that grass grows back instantly, which was also taken in [10]. Therefore, of which the terms are defined above.
The system of equations locally obeys the aforementioned mathematical constraint for The conservation equation S2 implies that the sum of all loss and gain terms should be zero. We can see that this is the case because: (i) total expansion due to successful 138 competition is at the cost of total loss due to unsuccessful competition, (ii) recruitment lost by S is gained by T , (iii) tree cover losses due to fire, drought and human impact , we need to take into account not only the spatial heterogeneity of these variables but also the relevant spatial interactions. 146 1. Spatial heterogeneity. Climatic, edaphic and anthropogenic spatial heterogeneity can be included by taking A and z as functions of space A(x) and z(x). 148 2. Spatial interaction. We assume that diffusion of cover types only occurs due to spread of seeds. We do not model seed dispersal but approximate it by dispersal 150 of saplings. Hence the diffusion coefficient of savanna adult tree cover is zero.
That of forest cover is not zero because part of its population is in the sapling 152 stage. Hence, The cover types that diffuse from neighboring areas settle in the areas that are 154 taken by grasses. Therefore, in the spatial model, the diffusion terms are also grass cover loss terms due to unsuccessful competition, or Spatial interaction can also occur due to spread of fire. While analysis of the model with fire spread is a bit more involved [19], the conclusions are the same.

Forest-savanna model
Here, we develop the forest-savanna model with all previously mentioned cover types.
Filling in the gains and losses and making use of Forest growth rate r F Here, we will show how we derived the maximum forest growth rate. We do this to set We will further set m F,o ≡ m and r F ≡ r. The ODE can be solved by separation of 178 variables using partial fractions such that the time that forest needs to grow from F 0 to If we take as initial tree cover a small value that could result from noise and as final tree cover the carrying capacity, we have 182 F 0 = 0.01, F 1 = 0.8.
At carrying capacity, we have r(1 − F )F − mF = 0,such that of which only the first is stable. As the data shows that F * has to be equal to 0.8, such 184 that m = 0.2r. Substituting this and the values for F 0 , F 1 , we obtain