*Article* **Continuum Modeling of Discrete Plant Communities: Why Does It Work and Why Is It Advantageous?**

**Ehud Meron 1,2,\*, Jamie J. R. Bennett 1, Cristian Fernandez-Oto 3, Omer Tzuk 1, Yuval R. Zelnik <sup>4</sup> and Gideon Grafi <sup>5</sup>**


Received: 16 September 2019; Accepted: 15 October 2019; Published: 17 October 2019

**Abstract:** Understanding ecosystem response to drier climates calls for modeling the dynamics of dryland plant populations, which are crucial determinants of ecosystem function, as they constitute the basal level of whole food webs. Two modeling approaches are widely used in population dynamics, individual (agent)-based models and continuum partial-differential-equation (PDE) models. The latter are advantageous in lending themselves to powerful methodologies of mathematical analysis, but the question of whether they are suitable to describe small discrete plant populations, as is often found in dryland ecosystems, has remained largely unaddressed. In this paper, we first draw attention to two aspects of plants that distinguish them from most other organisms—high phenotypic plasticity and dispersal of stress-tolerant seeds—and argue in favor of PDE modeling, where the state variables that describe population sizes are not discrete number densities, but rather continuous biomass densities. We then discuss a few examples that demonstrate the utility of PDE models in providing deep insights into landscape-scale behaviors, such as the onset of pattern forming instabilities, multiplicity of stable ecosystem states, regular and irregular, and the possible roles of front instabilities in reversing desertification. We briefly mention a few additional examples, and conclude by outlining the nature of the information we should and should not expect to gain from PDE model studies.

**Keywords:** continuum models; partial differential equations; individual based models; plant populations; phenotypic plasticity; vegetation pattern formation; desertification; homoclinic snaking; front instabilities

#### **1. Introduction**

Global warming and the concomitant increased frequency of intense droughts threaten the viability of plant populations and communities throughout the world [1–3]. As plants are primary producers that constitute the basal levels of whole food webs, that threat extends to ecosystem function as well. Understanding ecosystem response to drier climates, therefore, calls for unraveling mechanisms by which plant populations tolerate water stress. Such mechanisms operate at the organism level through various forms of phenotypic changes, but also at higher organization levels

through self-organization in spatial patterns. These response forms are particularly relevant in dryland ecosystems, where both phenotypic changes and spatial self-organization are highly significant.

Empirical studies of plant-community dynamics are hampered by the long time-scales of plant growth and the yet longer time scales of self-organization in space. A powerful methodology that complements empirical studies and compensates for the time-scale limitation is the construction and study of dynamic mathematical models. These models can be divided into two major groups. The first group consists of individual-based models (IBM) (also called agent-based models), which are stochastic computational algorithms for interacting individual organisms. Each individual is described by a set of time-dependent attributes, often including its spatial location on a grid and various physiological and behavioral traits [4–6]. The second group of models consists of continuum partial differential equations (PDE) that do not address discrete individual organisms, but rather deterministic processes at small spatial scales. The population is then described by a variable that represents the population size, such as the population number density, and is considered to be continuous in time and space [7–9].

The advantage of PDE models over IBM is that they lend themselves to the powerful methods of dynamical-system and pattern-formation theories, and therefore can provide deeper insights into the mechanisms that drive ecosystem dynamics and ecosystem response to environmental variability [8,10,11]. However, the application of the PDE modeling approach to plant populations in drylands should first be justified, as these populations are typically small, especially in arid regions, where the vegetation is sparse. This inherent character of drylands questions the use of continuous variables to describe population sizes; variables, such as number densities, become highly discrete, and demographic noise and extinction may become important aspects of the dynamics [12,13].

In this paper, we first argue that the discrete nature of small plant populations may still be discarded in modeling their dynamics owing to the high phenotypic plasticity of plants. That calls for a different description of the population size; rather than describing it by a discrete number density, we describe it by a biomass variable, which remains continuous even at the level of a single plant. We further demonstrate with a few examples the advantage of continuum PDE modeling over discrete IBM and discuss the ecological significance of this advantage. The examples include the use of linear stability analysis to calculate thresholds for the emergence of periodic patterns, the use of numerical continuation to calculate bifurcation diagrams, and a study of a transverse instability of desertification fronts. We conclude the paper with a brief discussion of additional examples that demonstrate the utility of continuum PDE models in gaining mechanistic information and insights about large, landscape-scale behaviors.

#### **2. Distinctive Aspects of Plant Populations**

Population size is often described by the number density of the individuals that constitute the population. We describe here two distinctive properties of plants, not shared by most other organisms: high phenotypic plasticity and seed dispersal, that will be used in Section 3 to motivate a better way of describing plant population size, especially in drylands where the vegetation is often sparse.

#### *2.1. Phenotypic Plasticity and Plant Adaptation to Variable Environments*

Phenotypic plasticity is defined as the ability of an organism (or a given genotype) to give rise to distinct observable traits (phenotypes) when exposed to variable environmental conditions [14]. Two basic, distinct developmental approaches have been evolved in plants and animals reflecting different degrees of phenotypic plasticity. Animals have a relatively short period of embryonic phase prior to birth, whereby all major organs are formed and postnatal development is essentially restricted to expansion and growth of existing organs. By contrast, plants spend most of their life, from weeks and months to hundreds and often thousands of years, at the embryonic phase, in the sense that they can produce new organs (stems, leaves, flowers) throughout their life, depending on environmental conditions. The higher phenotypic plasticity of plants can be attributed to their sessile nature and their inability to migrate to less stressful conditions as mobile animals do.

A plant embryo has only a small fraction of the organs an adult body has, consisting of two apical meristems located at the tip of the shoot and the root, which are responsible for growth and expansion of the shoot and the root systems. Both meristems consist of pluripotent cells (cells capable of differentiation into multiple cell types that make up the plant body) that are functionally equivalent to animal stem cells [15]. The shoot apical meristem (SAM) gives rise to the production of the aerial part of the plant including leaves and axillary buds (from which new stems emerge) as well as flowers. The root apical meristem (RAM) allows the growth and development of the root system and its expansion below ground. Thus, while the SAM activity leads to the expansion of the photosynthetic activity, which is carbon fixation by light energy, the RAM activity enhances uptake and mobilization of water and minerals to the canopy.

The proportion between aboveground (shoot system) biomass and belowground (root system) biomass, the so-called root-to-shoot ratio, is dynamic and needs to be balanced to achieve optimal performance (and therefore survival) under variable environmental conditions. Commonly there is a positive correlation between shoot growth and root growth, and both are interconnected [16]. Accordingly, shoot growth provides ample energy for the expansion of the root system, which in turn increases water and nutrient uptake from a larger belowground domain to feed and enable the expansion of the shoot system, a mechanism that plays an important role in the development of vegetation patterns (see root-augmentation feedback in Section 3.1.2).

The capacity of plants to dynamically change the allocation of biomass to different organs is central to plant response to variable environments and crucial for plant survival. A well-known example of periodic change in above-ground (shoot) biomass is the shedding of leaves in temperate deciduous forests during the winter (when leaves are susceptible to cold and are not photosynthetically active) as a mechanism to tolerate low winter temperature, e.g., saving energy via remobilization and storage of leaf nutritional constituents in stems. Shedding of leaves may also occur during the dry season in tropical and subtropical deciduous forests as a mechanism to tolerate seasonal drought by reducing the loss of water through transpiration.

Plants thriving in variable desert environments show many additional mechanisms to cope with seasonal climate variations that involve changes in above-ground biomass. A good example of such mechanisms is found in *Zygophyllum dumosum* Boiss (bushy bean caper), which is well adapted to a variable, desert environment [17,18]. Its root system is composed mostly of lateral, rapidly growing roots in the upper soil layer that constitute the major active elements in absorbing water [19], but also a few roots that can extend several meters in depth [20]. The shrub develops new stems with compound leaves during the wet season, where each leaf consists of two leaflets carried on a thick, fleshy petiole (Figure 1A). On entry into the summer, the above-ground biomass is altered due to the shedding of leaflets (to reduce whole plant transpiration) while the fleshy, wax-covered petioles remain alive (Figure 1B,C) [18].

Another strategy employed by *Z. dumosum*, particularly in successive drought years, is splitting the main axis (Figure 1D,E), a common phenomenon in desert shrubs [21] that enables the survival of certain units at the expense of others [22]. Also, certain desert shrubs such as *Artemisia sieberi* Besser (synonym name: *A. herba-alba*) and *Encelia farinosa* can change their mode of development and produce new leaf types on the transition from the wet to dry season, which is better adapted to dry conditions [23,24].

**Figure 1.** Phenology of *Zygophyllum dumosum* Boiss growing on a southeast-facing slope at Sede Boqer research area (30o 51 N 34o 46 E; elevation 498 m). (**A**) A typical branch of *Zygophyllum* plant during the wet season (March 2008). Note the coumpond leaf composed of two leaflets (LL) carried on a fleshy, cylindrical petiole (P). 1YOP, 1-year-old petiole. (**B**) A typical *Z. dumosum* branch during the dry season carrying petioles (P). (**C**) A closer look at emerging new bud from the axil of a 1YOP at the beginning of the winter. (**D**) Axis splitting resulting in a dual appearance of a *Z. dumosum* plant showing healthy branches carrying new leaves and flowers (right from the broken line) and unhealthy branches (left from the broken line). (**E**) A *Z. dumosum* plant with the main axis divided into several distinct units.

#### *2.2. Seed Persistence in the Soil*

Seeds persist in the soil until they germinate or die as a result of aging, predation or decay by fungi or bacteria [25]. Based on their longevity in soil, soil seed banks are divided into transient and persistent types, whereby the latter refers to seeds that remain viable for more than one year (Thompson and Grim, 1979) and often persist for decades [26] and even centuries [27]. There are several attributes that assist in seed persistence and longevity in the soil including dormancy, the capacity to repair accumulated damages and to possess active and/or passive defense mechanisms against potential predators and pathogens (reviewed in Ref. [28]).

Seed persistence is an important factor controlling the survival of a species long after the death of the mother plant. It avoids germination under unfavorable conditions (bet-hedging strategy, [29]) and allows for genetic preservation and distribution in time and space [30,31]. The capacity of seed persistence in the soil has implications for weed management, flora restoration and for understanding plant community dynamics, particularly in light of global climate change [31,32]. Persistence in the soil varies significantly between plant species and is dependent on the physical, physiological, chemical and biochemical properties of the dispersal unit and its capacity to withstand variable biotic and abiotic conditions [28].

These persistence characteristics may be altered upon exposure of mother plants to various biotic and abiotic stress conditions during seed development and prior to dispersal. Some dispersal unit characteristics underlying seed persistence in the soil are of maternal origin, embedded within the dead organs enclosing the embryos (DOEEs) including the seed coat (dry, dehiscent fruits), the pericarp (dry, indehiscent fruits) and dead floral bracts (glumes, lemmas, paleas) in grasses [33]. DOEEs are thought to function in seed dispersal and in protecting the embryo from mechanical and physical damages. However, detailed study of DOEEs revealed their capacity to store substances such as proteins (e.g., hydrolytic enzymes), growth factors (e.g., phytohormones) and various metabolites that might affect various aspects of seed biology including longevity, germination, and seedling vigor [33].

#### **3. Modeling Dryland Vegetation**

The high phenotypic plasticity of plants and the consequent wide range of above-ground biomass values a single plant can assume, suggest a description of plant-population sizes in terms of a biomass density variable, continuously varying in time and space, rather than in terms of a discrete number-density variable, as is often done in studies of animal populations. Since overland water flow and soil–water diffusion are continuous processes too, a natural way of describing dryland-vegetation dynamics is in terms of systems of partial differential equations (PDE) for biomass and water variables. This continuum modeling approach avoids the introduction of factitious demographic noise, which is a concern in small, low-plasticity animal populations. It is also consistent with the little relevance of extinction events in small plant populations, unlike animal populations. Such events are unlikely because of dispersed seeds that can remain viable long after plant mortality takes place, and are capable of reviving the plant population once favorable environmental conditions resume. The residual exponentially-small biomass that follows vegetation decay and convergence to bare soil in model solutions can be viewed as representing long-lived seeds.

PDE models of dryland vegetation describe the size of a plant population by a biomass-density variable, *B*(*X*,*Y*, *T*), which stands for the above-ground biomass of plants per unit area. Here, *X*,*Y* are the spatial coordinates in the plane (in units of meters), *T* is time (in units of years), and *B* has units of kg/m2. The biomass *BdXdY* in a small area element *dXdY* may represent the contribution of a single, several or many plants, depending on the particular plant species and on the spatial scale over which *B* varies. Several PDE models for dryland vegetation have been proposed. The simplest of which is a single-variable model for the vegetation biomass [34], while more detailed models also include a water variable [35–37], or two water variables representing below-ground water per unit area, *W*(*X*,*Y*, *T*), and above-ground water per unit area, *H*(*X*,*Y*, *T*) (both in units of kg/m2) [38,39].

In order to account for the emergence of vegetation patterns from uniform vegetation, the models should capture positive feedback loops that are capable of inducing nonuniform instabilities of a uniform-vegetation state. The more detailed models are advantageous in that they capture several pattern-forming feedbacks of this kind and thus allow us to study the interplay between different feedbacks [40–42]. These models also introduce better defined and measurable parameters. When applied to particular ecological contexts these models often simplify considerably and allow further mathematical analysis [43]. We refer the reader to Ref. [43] for a detailed description of a three-variable vegetation model that captures three distinct pattern-forming feedbacks. Here, we briefly describe this feedbacks and presents two simplified versions of the model that will be used later on to demonstrate the advantages of PDE models in studies of dryland ecosystems.

#### *3.1. Pattern-Forming Feedbacks*

Pattern-forming feedbacks in flat terrains follow the general scheme illustrated in Figure 2, that is, a positive feedback loop between local vegetation growth and water transport towards the growing vegetation [11,44]. They differ from one another in the mechanism of water transport: overland water flows, conduction of water by laterally extended root systems, and soil–water transport, as explained below. The transport of water towards denser vegetation patches accelerates vegetation growth there, and, at the same time, inhibits the growth in their surroundings. It, therefore, favors the growth of nonuniform perturbations. This is a short-range activation—long-range inhibition mechanism of

pattern formation [45], also termed "scale-dependent feedback" [46], where biomass is the activator and lack of water is the inhibitor. Associated with the three water-transport mechanisms are three different positive feedback loops as described below.

**Figure 2.** A general positive feedback loop that drives the formation of vegetation patterns in drylands. The feedback concomitantly accelerates vegetation growth in patches of denser vegetation and inhibits the growth in adjacent sparser patches. The combined processes favor the growth of nonuniform perturbations and the formation of vegetation patterns. From [11].

#### 3.1.1. Infiltration Feedback

Infiltration rates of surface-water into the soil in sparsely vegetated areas are typically lower than those in densely vegetated areas. Two main factors contribute to that effect: soil crusts in bare-soil that reduce the infiltration rate [47,48] and denser roots in denser vegetation that make the soil more porous and increase the infiltration rate. The infiltration contrast that builds up as the vegetation becomes denser in a given location induces overland water flow towards that location, which accounts for the upper arrow in Figure 2, i.e., enhancement of water transport by local vegetation growth. The increased soil moisture in the growth location further increases the rate of vegetation growth (lower arrow in Figure 2), which completes the positive feedback loop and sets the ground for a yet stronger positive-feedback loop. This pattern-forming feedback is referred to as the infiltration feedback.

To capture the infiltration feedback in the model equations we assume a monotonically increasing dependence of the infiltration rate *I* on the biomass variable *B* [38,39,49],

$$I = A \frac{B + Qf}{B + Q}.\tag{1}$$

The parameter *Q* controls how fast the asymptote *I* = *A* is approached. The parameter *f* ∈ [0, 1] controls the infiltration contrast, and, thus, the strength of the infiltration feedback. The value *f* = 1 corresponds to a constant infiltration rate, *I* = *A*, or no infiltration contrast, while values *f* 1 correspond to high infiltration contrasts, for which the infiltration rate in bare soil is significantly higher than in vegetated soil. This feedback may apply to ecosystems in which bare soil is covered by soil crusts, as these crusts typically reduce the infiltration rate relative to vegetation patches [48]. It is not expected to apply to uncrusted sandy soil, where the infiltration rate is high everywhere.

#### 3.1.2. Root-Augmentation Feedback

As noted in Section 2.1, there is a positive correlation between root growth and shoot growth that is quantified by the so-called root-to-shoot ratio. This property constitutes another mechanism by which water transport is enhanced by vegetation growth, as the lateral extension of the roots as the shoot grows enables water uptake and conduction from a larger volume. These processes and the consequent accelerated vegetation growth define the root-augmentation feedback.

The root-augmentation feedback is modeled by introducing a biomass-dependent kernel function,

$$G(\mathbf{X}, \mathbf{X}') = G\left(\frac{|\mathbf{X} - \mathbf{X}'|}{S[B(\mathbf{X})]}\right). \tag{2}$$

This function describes the spatial distribution of the roots in the horizontal directions *X* and *Y*, where **X** = (*X*,*Y*) represents the plant (shoot) location and **X** a distant point. The function *S*(*B*) determines the width of the kernel function and provides a measure for the horizontal extension of the roots. The root-augmentation feedback is captured by letting *S* increase monotonically with the biomass variable *B*(**X**), which represents the shoot mass. Figure 3 shows an example of a root kernel and its lateral extension as the plant grows and the above-ground biomass increases. The mathematical form of the kernel function is given by Equations (A1) and (A2) in Appendix A.

**Figure 3.** Root system growth as shoot biomass increases. A plot of the axisymmetric root kernel, *G* (see (A1) in Appendix A), as a function of the space coordinate *X* when (**A**) *B*(*X*, *T*) = 0, (**B**) *B*(*X*, *T*) = 0.2, (**C**) *B*(*X*, *T*) = 0.3, (**D**) *B*(*X*, *T*) = 0.4 =: *K*. The kernel is a Gaussian function multiplied by a polynomial factor, that mimics adventitious root branching observed in some plant species [50]. Parameters in (A1) and (A2): *<sup>S</sup>*<sup>0</sup> <sup>=</sup> 0.15m, *<sup>c</sup>*<sup>0</sup> <sup>=</sup> 1, *<sup>c</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>0.135, *<sup>c</sup>*<sup>4</sup> <sup>=</sup> 7.09 <sup>×</sup> <sup>10</sup>−3, *<sup>c</sup>*<sup>6</sup> <sup>=</sup> <sup>−</sup>9.11 <sup>×</sup> <sup>10</sup><sup>−</sup>5, *<sup>c</sup>*<sup>8</sup> <sup>=</sup> 3.84 <sup>×</sup> <sup>10</sup><sup>−</sup>7.

#### 3.1.3. Soil–Water Diffusion Feedback

Water depletion due to water uptake by the plants' roots, followed by soil–water diffusion from the patch surroundings, is a third mechanism of water transported that is enhanced by vegetation growth [51]. The associated feedback loop is referred to as the lsoil—water diffusion feedback. Like the root-augmentation feedback, this feedback relies on the root-to-shoot property of plants, except that here the role of the roots is to create soil–water gradients by local water uptake. This pattern-forming feedback can possibly apply to plants with vertical roots and strong water uptake, and to soil types for which lateral water diffusion is fast relative to the rate of vegetation spread.

#### *3.2. Mathematical Models*

For the sake of simplicity, we confine ourselves in this paper to ecosystems with sandy soil for which the infiltration rate is high also in bare-soil areas and, therefore, overland water flow can be neglected. Mathematically, this amount to assuming *f* = 1 (no infiltration contrast) and to the elimination of the equation for surface water. The full model then reduces to the following pair of equations for *B* and *W* [43]:

$$
\partial\_T B = G\_B B \left(1 - B/K\right) - MB + D\_B \Delta B,
$$

$$
\partial\_T W = P - L\_W W - G\_W W + D\_W \Delta W.\tag{3}
$$

Here, Δ is the Laplacian operator in the *X*,*Y* plane, *LW* is a biomass dependent evaporation rate,

$$L\_W = \frac{N}{\left(1 + RR/K\right)}\tag{4}$$

where *N* is the evaporation rate in bare soil and *R* quantifies the reduction in evaporation rate in vegetation patches, and *GB* and *GW* are the rates of biomass growth and water uptake by plants' roots, respectively. These rates are given by the integrals

$$G\_{B} = \Lambda \int\_{\Omega} G(\mathbf{X}, \mathbf{X}', T) \mathcal{W}(\mathbf{X}', t) \, \mathbf{d} \mathbf{X}' \,\, \,\tag{5a}$$

$$G\_W = \Gamma \int\_{\Omega} G(\mathbf{X}', \mathbf{X}, T) B(\mathbf{X}', t) d\mathbf{X}',\tag{5b}$$

over the system domain Ω, where the kernel function *G*(**X**, **X** , *t*) satisfies the general form (2), and Λ and Γ are rate constants in units of (kg/m2)−1y−1. These nonlocal forms represent the effects of laterally extended roots; the growth of a plant at point **X** depends on water availability at all points **X** within the reach of the plant's roots, and likewise, the water uptake at point **X** is due to plants at all points **X** whose roots extend to **X**. Two forms for the root kernel, *G*(**X**, **X** , *T*), will be used in this study. The first form represents vertical roots with negligible lateral extension. In that case the expressions in Equation (5a) for *GB* and *GW* simplify to algebraic expressions as shown in Appendix B. The second form represents laterally extended roots, as described in Appendix A.

Equation (3) capture several processes that affect biomass dynamics: water-dependent plant growth (*GBB*), mortality (−*MB*), and short-distance seed dispersal (*DB*Δ*B*). In addition, the late biomass growth phase is slowed down by species-specific constraints that dictate a maximal standing biomass per unit ground area *K*. These may represent self-shading, limited stem strength of a woody plant, etc. The processes that affect soil–water dynamics, according to Equation (3) are precipitation (*P*), water loss due to biomass dependent evaporation (−*LW*), water uptake by plants' roots (−*GWW*), and soil–water diffusion (*DW*Δ*W*). The simplified model equations do not capture the infiltration feedback, as no overland water flow takes place, but still capture the root-augmentation feedback and the soil–water diffusion feedback.

#### **4. Advantages of Continuum PDE Models**

Continuum PDE models are advantageous over discrete models, such as individual based models (IBM), in that they lend themselves to mathematical analysis, which results in deeper insights into the phenomena in question and the ecological implications they have. We chose to demonstrate this aspect using three examples, as described in the following subsections.

#### *4.1. Instability Thresholds*

The simplest solutions of Equation (3) are stationary uniform solutions obtained by solving the algebraic equations that result by setting the time and space derivatives of the state variables to zero. Two such solutions are found, (*B*, *W*)=(0, *W*0) representing bare soil and existing for all positive precipitation values, *P*, and (*B*, *W*)=(*B*1, *W*1) representing uniform vegetation and existing only beyond some precipitation threshold (see Figure 4). The stability of these solutions to infinitesimal perturbations of all forms can be calculated using linear stability analysis [11,52]. Such an analysis reveals that the bare-soil state is stable in the precipitation range 0 < *P* < *P*0, where *P*<sup>0</sup> is an instability threshold at which a uniform mode (characterized by zero wavenumber) begins to grow, as the dispersion curves in Figure A1A show. This is a uniform stationary instability [43] where the first mode to grow is spatially uniform and the growth is monotonic in time. Spatial patterning, induced by the root-augmentation feedback or the soil–water diffusion feedback, cannot arise from a nonuniform instability of the bare-soil state as this would involve the growth of a periodic mode from a zero-biomass state and therefore imply negative biomass values, which are unphysical. Failing to obtain a uniform instability of the bare-soil state would indicate a modeling flaw. Such a flaw would be harder to detect in discrete models, such as IBM.

Stationary periodic vegetation patterns may arise from a nonuniform stationary instability of the uniform vegetation state. A linear stability analysis of that state indeed reveals such instability as the precipitation rate *P* decreases below a threshold value *PT*. At that threshold a periodic mode with a characteristic wavenumber *kT* begins to grow monotonically in time while all other modes still decay, as the dispersion curves in Figure A1B shows. This analysis can be performed even for quite complicated laterally extended root kernels *G* as shown in Figure 3, see Appendix A. The analysis further confirms the expectation that each of the pattern-forming feedbacks alone, root-augmentation and soil–water diffusion, can induce a nonuniform stationary instability that culminates in stationary periodic patterns.

Linear stability analysis not only unravels possible instabilities and the nature of the modes that grow beyond instability points, but it also provides information about the parameters that control the instability threshold. For example, the instability of the bare-soil state occurs at *P* = *P*<sup>0</sup> = *MN*/Λ, indicating that species with low ratios of growth rates, Λ, to mortality rates, *M*, will grow from bare soil at higher precipitation thresholds, or that high evaporation rates, *N*, act to stabilize the bare soil state. Such relations are less apparent in IBM.

#### *4.2. Bifurcation Diagrams*

Linear stability analysis does not provide information about the new state that the system converges to following an instability. This kind of information can be obtained by nonlinear analysis of the model equations. A common approach, valid close enough to the instability point, involves the derivation of amplitude equations [11,53–55]. These are nonlinear differential equations for the amplitudes of the modes that grow beyond the instability point. Amplitude equations are, in general, much easier to analyze than the original model equations. Moreover, they are dictated by the instability type, rather than by the specific model equations, and therefore are universal; different systems that go through the same instability type (e.g., nonuniform stationary instability) will behave similarly near the instability point. Since all pattern-forming feedbacks lead to the same instability type they will also induce similar dynamics and patterns, although fine details, such as the relative spatial distribution of biomass and water, can reflect differences in the these feedbacks [56].

A convenient way to summarize the outcomes of nonlinear analysis is to draw a bifurcation diagram. Such a diagram shows the existence and stability ranges of various model solutions. In the context of dryland vegetation it often shows graphs of spatial biomass average (or L2 norm) vs. precipitation rate for selected model solutions, using the convention that solid (dashed) lines represent stable (unstable) solutions. If amplitude equations are known, a bifurcation diagram can often be calculated analytically [11,57,58]. When the derivation of amplitude equations is not easy to perform, or when there is an interest in dynamical behaviors far from instability points, bifurcation diagrams can be calculated by numerical continuation methods using software packages such as AUTO [59], pde2path [60], and others. Since these packages solve the equations using iterative solution methods, they converge to unstable solutions as well and, thereby, can provide complete solution branches in the bifurcation diagram.

The utility of bifurcation diagrams can be demonstrated with the example shown in Figure 4, a bifurcation diagram for a vegetation model (see Appendix B) in one space dimension. We focus here on one property that bifurcation diagrams often reveal, namely, parameter ranges where multiple stable states coexist. The (partial) bifurcation diagram of Figure 4 shows a bistability range of bare soil and periodic patterns, a bistability range of bare soil and uniform vegetation, and a tristability subrange of bare soil, periodic patterns and uniform vegetation. The stability of the bare soil state at relative high precipitations values, where uniform vegetation is stable too, can be achieved assuming high evaporation rates (*N*) in bare soil. Within the tristability range there exists many more stable states, which describe spatial hybrids of periodic pattern and uniform vegetation. These hybrid states appear as confined domains of an increasing size of the patterned state in an otherwise uniform vegetation state, as the insets in Figure 4 show. Hybrid states in bistability ranges of uniform and patterned states have been studied extensively, both mathematically and in particular physical systems. The confined pattern domains are homoclinic solutions that snake back and forth in the bifurcation

diagram, a behavior termed "homoclinic snaking" [61–66]. They are tightly related to front pinning, that is, to the existence of stationary front solutions between uniform and patterned states over a range of the control parameter. This is in contrast to fronts between distinct uniform states that propagate in general, except, possibly, for a particular parameter value, the so-called Maxwell point at which the fronts are stationary [67,68], and unlike fronts that are pinned by external heterogeneities [69].

**Figure 4.** Bifurcation diagram for the model presented in Appendix B (Equation (A17)) in 1D (panel **a**). The vertical axis is the spatial biomass average, while the horizontal axis is the precipitation rate. Solid (dashed) lines represent stable (unstable) solutions. The diagram shows bistability precipitation ranges of bare soil and periodic patterns, and of bare soil and uniform vegetation, as well as a tristability range of bare soil, uniform vegetation and periodic patterns. Within the latter range there exists a subrange of hybrid patterns, consisting of confined pattern domains of increasing size in an otherwise uniform vegetation, examples of which are shown in the three panels (**b**–**d**). Parameters: Λ = 0.1 [(m2/kg)/y], Γ = 4 [(m2/kg)/y], *E* = 7 [m2/kg], *K* = 1 [kg/m2], *M* = 4.5 [1/y], *N* = 8 [1/y], *R* = 0.1, *DB* = 0.05 [m2/y], *DW* = 30 [m2/y].

Homoclinic snaking has been used to explain fairy-circle dynamics in Namibian grasslands [56,70–72]. Fairy circles are circular bare-soil gaps in grasslands that often form nearly periodic hexagonal patterns, where each gap is surrounded by nearly equidistant six other gaps on average [73]. Empirical studies indicate that fairy circles are occasionally "born" or "die" locally [74]. Taking into account the high rainfall inter-annual variability in fairy-circles sites, such dynamical events have been interpreted as transitions between different hybrid states (homoclinic solutions) induced by rainfall fluctuations that are strong enough to drive the ecosystem temporarily outside the existence range of these states (snaking range) [51]. The existence of a multitude of stable hybrid states makes dryland landscapes more plastic in the sense that their response to varying environmental conditions and disturbances can involve temporal convergence to different hybrid states, rather than direct convergence to a single alternative stable state [75].

In order for homoclinic snaking to occur in models of dryland landscapes bistability of uniform vegetation and periodic vegetation patterns are generally sufficient [75,76]. The snaking range in the bifurcation diagram of Figure 4, however, lies within a tristability range of uniform vegetation, periodic pattern and (uniform) bare soil. The existence of that tristability range may have interesting consequences. In particular, the organization of hybrid states within a homoclinic snaking structure can break down as it meets a Maxwell point, where fronts connecting the bare soil and uniform vegetation states are stationary [77]. Another implication of the tristability of two uniform states and a periodic pattern state is the emergence of a family of complex front structures that involve all three states [77].

An additional important aspect that bifurcation diagrams uncover is related to the existence of unstable solution branches. Tracking unstable solutions in bifurcation diagrams can be highly significant when the response of ecosystems to disturbances or human intervention is studied, as the appearance or disappearance of unstable solutions can dramatically affect the flow in phase space and, thus, the response [78]. These subtle aspects, which may have significant ecological implications, become apparent once a bifurcation diagram is calculated. Bifurcation diagrams showing stable and unstable uniform, periodic and localized solutions can be calculated using PDE models [60] but, practically, not with IBM. Bifurcation diagrams have been calculated using IBM but the information they contain is very limited [79,80].

#### *4.3. Front Dynamics*

Bistable ecosystems can go through state transitions, or regime shifts [81], in various ways, including a passage through a bifurcation point (B-tipping), as a result of environmental fluctuations (N-tipping), or as a result of fastly varying environmental conditions (R-tipping) [82,83]. Such transitions are generally discussed as whole-system responses, occurring simultaneously at all points in space. Spatially confined disturbances, however, can induce local state transitions. The dynamics that follow such disturbances are dictated by the motion of the fronts that bound the transition area. Of particular interest are degradation fronts, where a dysfunctional state gradually displaces a functional state by front propagation. An example of such a front in the context of dryland vegetation is a desertification front, where a bare-soil domain displaces a vegetated domain [68,84].

PDE models are highly valuable in analyzing desertification fronts and addressing questions such as how to reverse the process of desertification, as we briefly discuss below. We consider again the bifurcation diagram shown in Figure 4, and focus on a bistability precipitation range of bare soil and uniform vegetation, which may or may not include the tristability range where periodic patterns are stable too. In that range we consider precipitation values below the Maxwell point, where desertification fronts exist (bare soil displacing uniform vegetation). In the vicinity of the bare-soil instability to uniform vegetation (*P* = *PC* in Figure 4) the two-variable model presented in Appendix B can be reduced to an amplitude equation for the uniform mode that begins to grow at this instability [85]. Analysis of this equation reveals a transverse front instability [86,87], whereby small bulges along the front line are first enhanced, then develop into growing fingers that avoid one another, and ultimately fill up the system domain with a stationary labyrinthine pattern.

According to the analysis of the amplitude equation, the instability occurs as the soil–water diffusion coefficient, *DW*, exceeds a threshold value, and that threshold is inversely related to the root-to-shoot ratio, *E*, which controls water uptake by plants' roots [85]. This result uncovers the mechanism of the instability—fast soil–water diffusion (relative to biomass expansion) towards incidental bulges along the front line that locally deplete the soil–water content. The fast diffusion acts to enhance the growth of these bulges and, at the same time, to inhibit vegetation growth on both sides of any bulge. Fast diffusion can be obtained by increasing the diffusion coefficient, *DW*, which is consistent with the finding of a threshold value above which the instability develops. Fast diffusion is also obtained by steepening the soil–water gradients, which can be achieved by strong water uptake. This is consistent with the finding that the threshold value decreases as the parameter *E* that controls water uptake, is increased.

Figure 5 shows a desertification front that develops a transverse instability. While the bare soil area initially expands into the vegetated area, the instability results in vegetation fingers that grow backward into the bare-soil area and thereby reverse the desertification process. The resulting state is a

productive vegetation pattern that prevents further irreversible degradation processes (not captured by the model), such as soil erosion, and maintains the ecosystem in a reversible state capable of forming uniform vegetation when favorable rainfall conditions resume. But how can a transverse instability be induced in stable desertification fronts? One possible answer to this question is the introduction of a species with sufficiently high root-to-shoot ratio so as to reduce the threshold value of *DW* above which the instability sets in, and thereby induce a transverse instability. There is, however, another possibility associated with the tristability range of uniform vegetation, bare soil and periodic patterns. In this range, desertification fronts that are stable to small perturbations (linearly stable) may still be unstable to larger perturbations (nonlinearly unstable), which drive the system to the periodic-pattern state through finger growth [88] as Figure 6 demonstrates.

These results are very appealing from the point of view of ecosystem management. First, they imply local manipulations in the front zone only, rather than extensive intervention across the whole ecosystem. The manipulations may involve planting a different species capable of inducing a linear front instability, or modulating the front line strongly enough in order to induce a nonlinear front instability. The latter intervention form may involve periodic grazing management, clear-cutting or irrigation, and is relevant in cases when indications for the existence of periodic patterns exist (e.g., scattered patches of periodic vegetation). Second, such manipulations are limited in time, as they are needed only to trigger the instability. Once the instability sets in, a process of self-recovery begins. These conclusions, which are based on the mathematical analysis of front solutions, could not have been obtained in studies of IBM.

**Figure 5.** Snapshots of model solutions (see Appendix B) that demonstrate a linear front instability of a desertification front. Following a short phase in which the bare-soil domain (yellow) expands into the uniform-vegetation domain (green), transverse perturbations begin to grow and form vegetation fingers that grow back into the bare-soil domain. The time indicated in every snapshot is in units of years. From [85] .

**Figure 6.** Nonlinear front instability in a tristability range. Snapshots of numerical model solutions showing (**a**) the stability of a planar desertification front to small transverse modulation, (**b**) instability to transverse modulations that are sufficiently large, and the development of vegetation fingers that grow back into bare soil. The time indicated in every snapshot is in units of years. From [85].

#### **5. Discussion**

Despite their typically small size, discrete plant populations in dryland ecosystems can still be described by continuum PDE models because of the high phenotypic plasticity of plants and the dispersal of stress-tolerant seeds, as discussed in Sections 2 and 3. Since the state variables that describe plant population sizes should reflect that plasticity, biomass densities are more appropriate choices for plant populations than the often used number densities in population dynamics studies. Unlike IBM, where the smallest entity is a single individual, in PDE models with biomass densities as state variables, the smallest entity is a small area element and the processes that take place there. That area element can be significantly smaller than the scale of a single plant and its roots.

We presented here several examples that demonstrate the utility of continuum PDE models in gaining mechanistic information and insights about large, landscape-scale behaviors, such as the onset of a nonuniform stationary instability that culminates in periodic vegetation patterns, uncovering precipitation ranges of bistability, tristability and multistability of uniform states, periodic patterns and localized patterns (hybrid states), and dynamics of desertification fronts, where we focused on front instabilities that may reverse gradual desertification. Many more examples exist, a few of them are briefly discussed below.

On a slope, a linear stability analysis of uniform vegetation has revealed a nonuniform oscillatory instability that results in traveling-wave solutions. These solutions describe periodic vegetation stripes migrating uphill [89,90], as observed in empirical studies [91]. The migration mechanism is easy to understand; while plants at the top part of a vegetation stripe receive runoff and grow, plants at the bottom part loose runoff and die. PDE models have been instrumental in clarifying the relations between migration speed, pattern wavelength, and slope [92,93].

Studies of PDE models in one space dimension have identified many periodic solutions along the rainfall gradient, the wavenumbers of which decrease as precipitation decreases, down to zero in a solution that represents a single vegetation patch [94–96]. These solutions reflect two manners by which patchy vegetation responds to water stress. The first is an increase of bare soil areas at the expense of vegetation-patch areas, keeping the wavenumber unchanged. This response occurs along each solution branch [97]. The second response is a transition to a periodic solution with a lower wavenumber, quite often half the original wavenumber [95–97]. In this response, the area of each vegetation patch does not necessarily decrease, but the number of vegetation patches decreases. In both response forms, the increase in bare-soil area compensates for the reduced rainfall through increased water transport to adjacent vegetation patches by one of the transport mechanisms discussed in Section 3.1. Empirical indications for the existence of periodic patterns with different wavenumbers have been found in studies of banded vegetation in Somalia [98].

In two space dimensions, patchy vegetation can respond to water stress by yet another mechanism of increasing bare-soil areas, namely, morphological changes, first from hexagonal gap patterns to stripe patterns, and, at lower precipitation, from stripe patterns to hexagonal spot patterns [36,43,99]. In cases where the instability of bare soil to uniform vegetation is supercritical, weak nonlinear analyses of PDE models in two space dimensions have produced bifurcation diagrams that unfold the basic periodic vegetation patterns along the rainfall gradient: hexagonal gap patterns, stripe patterns and hexagonal spot patterns [57,58]. Empirical indications for the existence of these vegetation patterns in nature are abundant, although not yet with a single plant species along a rainfall gradient [56,73,100–102].

PDE models have also been used to study the interactions between two distinct plant species. One interesting context of such interactions is the bistability of a uniform state of one species and a periodic-pattern state of the other species. In this range, homoclinic snaking can result in the multistability of hybrid states that involve the two species, and therefore constitutes a mechanism for species coexistence [103]. Another interesting context is woody-herbaceous systems, where the woody species is pattern forming. Studies of a PDE model that captures both the infiltration feedback and the root-augmentation feedback reveal changes in the relative importance of the two feedbacks along the rainfall gradient. At high precipitation dominance of the root augmentation feedback results in a strong depletion of the soil–water content and the exclusion of herbaceous vegetation. At low precipitation, the dominance of the infiltration feedback results in soil–water concentration at woody patches and the consequent facilitation of herbaceous-vegetation growth [41]. These results are consistent with the stress-gradient hypothesis in ecology [104].

While these results are hard to obtain using IBM, we should not regard PDE models as substitutes to IBM, but rather as complementary means to gain deeper mechanistic understanding. IBM typically go into much more details, and attempts to include these details in PDE models would render them mathematically intractable too. PDE models should rather be motivated by specific questions of interest and by judicious assessments of the processes that are most relevant to these questions. The vegetation models described in Section 3 have been motivated by empirical observations of regular vegetation patterns in drylands and by the understanding that the most relevant processes to these phenomena are positive feedback loops involving vegetation growth and water availability that can induce pattern-forming instabilities. Various other processes, considered to be of secondary significance, have been left aside, such as the effect of transpiration on the atmosphere, soil erosion and deposition, and various plant physiology processes. Such processes are likely to have quantitative rather than qualitative effects on vegetation pattern formation; they may affect, for example, instability thresholds, but are not likely to affect the occurrence of instabilities or change their nature. While providing deep mechanistic insights into ecological processes, and predicting possible qualitative responses to environmental changes, PDE models should not be regarded as tools for making quantitative forecasts.

**Author Contributions:** Conceptualization, E.M.; methodology, all authors; software, Y.R.Z. and O.T.; formal analysis, J.J.R.B. and C.F.-O.; investigation, all authors; writing—original draft preparation, E.M. and G.G.; writing—review and editing, all authors; visualization, all authors; funding acquisition, E.M.

**Funding:** This study was funded by the Israel Science Foundation under grant number 1053/17, and by the European Union's Horizon 2020 research and innovation programme under grant agreement No 641762.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Model for Laterally Extended Roots**

We consider here the model Equation (3) in two space dimensions for a laterally extended root kernel (2) of the form

$$G(\mathbf{X}, \mathbf{X}', T) = \frac{F(\mathbf{X}, \mathbf{X}')}{2\pi S\_0^2} \mathbf{e}^{-\frac{|\mathbf{X} - \mathbf{X}'|^2}{2S\_0^2 (1 + E B(\mathbf{X}, T))^2}},\tag{A1}$$

where *S*<sup>0</sup> is the lateral root length of a seedling and *F* is a polynomial given by

$$F(\mathbf{X}, \mathbf{X}') = \left(c\_0 + c\_2 S\_0^{-2} |\mathbf{X} - \mathbf{X}'|^2 + c\_4 S\_0^{-4} |\mathbf{X} - \mathbf{X}'|^4 + c\_6 S\_0^{-6} |\mathbf{X} - \mathbf{X}'|^6 + c\_8 S\_0^{-8} |\mathbf{X} - \mathbf{X}'|^8\right) / \psi \,. \tag{A2}$$

Here *ψ* = *c*<sup>0</sup> + 2*c*<sup>2</sup> + 8*c*<sup>4</sup> + 48*c*<sup>6</sup> + 384*c*<sup>8</sup> is a normalization factor that ensures the integral of *G*(**X**, **X** , *T*) over the entire domain is unity when *B*(**X**, *T*) = 0. Since *G* has units of of *m*−<sup>2</sup> (see Equation (5a)), *<sup>F</sup>* is a dimensionless quantity. The "shape parameters" *<sup>c</sup>*0, *<sup>c</sup>*2, *<sup>c</sup>*4, *<sup>c</sup>*6, *<sup>c</sup>*<sup>8</sup> <sup>∈</sup> <sup>R</sup> and the normalization factor *ψ* are then dimensionless too. For proper choices of these parameters the root kernel can describe root branching as Figure 3 shows.

**Non-dimensionalisation**: Using the scalings

$$b = \frac{B}{K'}, w = \frac{W\Lambda}{K\Gamma}, \ t = MT, \ x = \frac{\mathbf{x}}{S\_0},\tag{A3}$$

one can non-dimensionalise model (3) to obtain

$$\frac{\partial b}{\partial t} = G\_b b (1 - b) - b + \delta\_b \Delta b,\tag{A4a}$$

$$\frac{\partial w}{\partial t} = p - \nu (1 - \rho b)w - G\_w w + \delta\_w \Delta w. \tag{A4b}$$

The non-dimensional biomass growth rate and water uptake rate are given by

$$\mathbf{G}\_{\mathsf{b}} = \lambda \int\_{\Omega} \mathbf{g}(\mathbf{x}, \mathbf{x}', t) w(\mathbf{x}', t) d\mathbf{x}',\tag{A5a}$$

$$\mathbf{G}\_{\rm w} = \lambda \int\_{\Omega} \mathbf{g}(\mathbf{x}', \mathbf{x}, t) b(\mathbf{x}', t) d\mathbf{x}',\tag{A5b}$$

with

$$\log(\mathbf{x}, \mathbf{x}', t) = \frac{f(\mathbf{x}, \mathbf{x}')}{2\pi} \mathbf{e}^{-\frac{|\mathbf{x} - \mathbf{x}'|^2}{2(1 + \eta b(\mathbf{x}, t))^2}},\tag{A6}$$

where

$$f(\mathbf{x}, \mathbf{x}') = (c\_0 + c\_2|\mathbf{x} - \mathbf{x}'|^2 + c\_4|\mathbf{x} - \mathbf{x}'|^4 + c\_6|\mathbf{x} - \mathbf{x}'|^6 + c\_8|\mathbf{x} - \mathbf{x}'|^8) / \psi. \tag{A7}$$

The relations between non-dimensional and dimensional quantities are given by

$$\lambda = \frac{K\Gamma}{M}, \eta = E\\K, \ p = \frac{\Lambda P}{K\Gamma M}, \nu = \frac{N}{M}, \rho = R, \ \delta\_b = \frac{D\_B}{MS\_0^2}, \ \delta\_w = \frac{D\_W}{MS\_0^2}$$

**Linear stability**: We use here linear stability analysis to study the stability properties of the solution of (3) that represents steady uniform vegetation, which we denote by **U**<sup>0</sup> = (*b*0, *w*0)T. In this method one considers (infinitesimally) small nonuniform perturbations about the state in question and study whether all perturbations decay to zero, or some perturbations grow, indicating an instability. We denote a small perturbation of this kind by **U˜** of **U**0, so that the perturbed state is given by

$$\mathbf{U}(\mathbf{x},t) = \mathbf{U}\_0 + \mathbf{\tilde{U}}(\mathbf{x},t),\tag{A8}$$

.

where **U** = (*b*, *w*)<sup>T</sup> and **U˜** = (˜ *b*, *w*˜)T. In the standard way [11] we let **U**˜ (**x**, *t*) = **a**(*t*)ei**<sup>k</sup>**·**<sup>x</sup>** + *c*.*c* where **a**(*t*)=(*a*1, *a*2)<sup>T</sup> and **k** = (*kx*, *ky*). Substitution of (A8) into (3) gives

$$\frac{\partial \tilde{b}}{\partial t} = \mathcal{G}\_b \big|\_{\mathcal{U}\_0 + \tilde{\mathcal{U}}} (b\_0 + \tilde{b})(1 - (b\_0 + \tilde{b})) - (b\_0 + \tilde{b}) + \delta\_b \Delta \tilde{b},\tag{A9a}$$

$$\frac{\partial \tilde{w}}{\partial t} = p - \nu (1 - \rho (b\mathbf{\hat{v}} + \tilde{b}))(w\mathbf{\hat{v}} + \tilde{w}) - G\_{\mathbf{w}}|\_{\mathbf{U}\_0 + \tilde{\mathbf{U}}}(w\mathbf{\hat{v}} + \tilde{w}) + \delta\_{\mathbf{W}} \Delta \tilde{w}. \tag{A9b}$$

To evaluate *Gb* **<sup>U</sup>**0+**U˜** and *Gw* **<sup>U</sup>**0+**U˜** we expand the kernels (A6) as

$$g(\mathbf{x}, \mathbf{x}', t) = g\_0(\mathbf{x}, \mathbf{x}') + g\_1(\mathbf{x}, \mathbf{x}')\tilde{b}(\mathbf{x}, t) + \mathcal{O}(\tilde{b}^2) \tag{A10a}$$

$$\mathbf{g}(\mathbf{x}',\mathbf{x},t) = \mathbf{g}\_0(\mathbf{x}',\mathbf{x}) + \mathbf{g}\_1(\mathbf{x}',\mathbf{x})\mathbf{\tilde{b}}(\mathbf{x}',t) + \mathcal{O}(\mathbf{\tilde{b}}^2) \tag{A10b}$$

where

$$g\_0(\mathbf{x} - \mathbf{x}') = g\_0(\mathbf{x}' - \mathbf{x}) = g(\mathbf{x} - \mathbf{x}', t)|\_{\mathbf{U} = \mathbf{U}\_0} = \frac{f(\mathbf{x}, \mathbf{x}')}{2\pi} \mathbf{e}^{-\frac{|\mathbf{u} - \mathbf{x}'|^2}{2v\_0^2}}\tag{A11a}$$

$$\log g\_1(\mathbf{x} - \mathbf{x}') = g\_1(\mathbf{x}' - \mathbf{x}) = \frac{\partial g}{\partial b}(\mathbf{x} - \mathbf{x}', t) \Big|\_{\mathbf{U} = \mathbf{U}\_0} = \frac{\eta |\mathbf{x} - \mathbf{x}'|^2 f(\mathbf{x}, \mathbf{x}')}{2\pi \sigma\_0^3} \mathbf{e}^{-\frac{|\mathbf{x} - \mathbf{x}'|^2}{2\sigma\_0^2}},\tag{A11b}$$

with *σ*<sup>0</sup> = 1 + *ηb*0. Considering linear terms only we can then calculate

$$\begin{split} \left. G\_{\mathsf{b}} \right|\_{\mathsf{U}\_{0} + \mathsf{D}} & \approx \lambda \mathsf{u}\_{0} \int\_{\Omega} g\_{0}(\mathbf{x}, \mathbf{x}') d\mathbf{x}' + \lambda \int\_{\Omega} g\_{0}(\mathbf{x}, \mathbf{x}') \vec{w}(\mathbf{x}', t) d\mathbf{x}' + \lambda \mathsf{u}\_{0} \vec{b}(\mathbf{x}, t) \int\_{\Omega} g\_{1}(\mathbf{x}, \mathbf{x}') d\mathbf{x}' \\ & = \frac{\lambda \mathsf{u}\_{0} \sigma\_{0}^{2} \psi\_{0}}{\Psi} + \lambda \mathcal{F}\_{0}(k) \vec{w}(\mathbf{x}, t) + \frac{\lambda \eta \mathsf{u}\_{0} \psi\_{1}}{\sigma\_{0} \psi} \vec{b}(\mathbf{x}, t), \\ \mathsf{C} & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \int\_{\Omega} g\_{1}(\mathbf{x}, \mathbf{x}', t) \mathsf{u}\_{1} \vec{w}(\mathbf{x}, t) \tilde{\mathsf{U}}(\mathbf{x}, t) \tilde{\mathsf{U}}(k) \vec{w}(\mathbf{x}, t) \cdot \int\_{\Omega} g\_{2}(\mathbf{x}, t) \tilde{\mathsf{U}}(\mathbf{x}, t) \tilde{\mathsf{U}}(\mathbf{x}, t) \vec{w}(\mathbf{x}, t) \end{split} \tag{A12a}$$

$$\begin{split} \left. \mathcal{G}\_{w} \right|\_{\mathbf{U} \rightarrow \mathbf{0}} & \approx \lambda b\_{0} \int\_{\Omega} \underline{\mathbf{g}\_{0}(\mathbf{x}', \mathbf{x}) d\mathbf{x}'} + \lambda \int\_{\Omega} \underline{\mathbf{g}\_{0}(\mathbf{x}', \mathbf{x}) \tilde{b}(\mathbf{x}', t) d\mathbf{x}'} + \lambda b\_{0} \int\_{\Omega} \underline{\mathbf{g}\_{1}(\mathbf{x}', \mathbf{x}) \tilde{b}(\mathbf{x}', t) d\mathbf{x}'} \\ &= \frac{\lambda b\_{0} \nu\_{0}^{2} \psi\_{0}}{\psi} + \lambda \mathcal{F}\_{0}(k) \tilde{b}(\mathbf{x}, t) + \lambda b\_{0} \mathcal{F}\_{1}(k) \tilde{b}(\mathbf{x}, t), \end{split} \tag{A12b}$$

where *<sup>k</sup>* = |**k**| is the wavenumber of the perturbation, *<sup>ψ</sup>*<sup>0</sup> = *<sup>c</sup>*<sup>0</sup> + <sup>2</sup>*c*2*σ*<sup>2</sup> <sup>0</sup> + 8*c*4*σ*<sup>4</sup> <sup>0</sup> + 48*c*6*σ*<sup>6</sup> <sup>0</sup> + 384*c*8*σ*<sup>8</sup> 0 , and *ψ*<sup>1</sup> = 2*c*0*σ*<sup>2</sup> <sup>0</sup> + <sup>8</sup>*c*2*σ*<sup>4</sup> <sup>0</sup> + <sup>48</sup>*c*4*σ*<sup>6</sup> <sup>0</sup> + <sup>384</sup>*c*6*σ*<sup>8</sup> <sup>0</sup> + <sup>3840</sup>*c*8*σ*<sup>10</sup> <sup>0</sup> . F0(*k*) and F1(*k*) are the Fourier transforms of *g*0(**x**, **x** ) and *g*1(**x**, **x** ), respectively, and are expressed in terms of the quantities

$$\begin{aligned} F\_0 &= \sigma\_0^2 \exp(-\sigma\_0^2 k^2 / 2), \\ F\_2 &= F\_0 \sigma\_0^2 (2 - \sigma\_0^2 k^2), \\ F\_4 &= F\_0 \sigma\_0^4 (8 - 8\sigma\_0^2 k^2 + \sigma\_0^4 k^4), \\ F\_6 &= F\_0 \sigma\_0^6 (48 - 72\sigma\_0^2 k^2 + 18\sigma\_0^4 k^4 - \sigma\_0^6 k^6), \\ F\_8 &= F\_0 \sigma\_0^8 (384 - 768\sigma\_0^2 k^2 + 288\sigma\_0^4 k^4 - 32\sigma\_0^6 k^6 + \sigma\_0^8 k^8), \\ F\_{10} &= F\_0 \sigma\_0^{10} (3840 - 9600\sigma\_0^2 k^2 + 4800\sigma\_0^4 k^4 - 800\sigma\_0^6 k^6 + 50\sigma\_0^8 k^8 - \sigma\_0^{10} k^{10}), \end{aligned}$$

as

$$\begin{aligned} \mathcal{F}\_0(k) &= \frac{1}{\Psi} (c\_0 F\_0 + c\_2 F\_2 + c\_4 F\_4 + c\_6 F\_6 + c\_8 F\_8), \\ \mathcal{F}\_1(k) &= \frac{\eta}{c\_0^3 \Psi} (c\_0 F\_2 + c\_2 F\_4 + c\_4 F\_6 + c\_6 F\_8 + c\_8 F\_{10}). \end{aligned}$$

The uniform steady states of (3) satisfy

$$\frac{\lambda \mu\_0}{\Psi} w\_0 b\_0 (1 - b\_0) \sigma\_0^2 - b\_0 = 0,\tag{A13a}$$

$$p - \nu w\_0 (1 - \rho b\_0) - \frac{\lambda \psi\_0}{\Psi} w\_0 b\_0 \sigma\_0^2 = 0,\tag{A13b}$$

which yield two homogeneous steady states: the bare soil state (0, *p*/*ν*) and the uniform vegetation state, the expression of which we omit here for brevity. Substitution of (A12) into (A9), retaining only first order terms, and using (A13) yields the system of equations

$$\frac{d\mathbf{a}}{dt} = \mathbf{J}\mathbf{a} \tag{A14}$$

where the components of **<sup>J</sup>** <sup>∈</sup> <sup>R</sup>2×<sup>2</sup> are given by

$$\mathbf{J}\_{11} = \frac{\lambda w\_0}{\psi} \left( \frac{\eta \psi\_1 b\_0 (1 - b\_0)}{\sigma\_0} + \sigma\_0^2 \psi\_0 (1 - 2b\_0) \right) - 1 - \delta\_b k^2,\tag{A15a}$$

$$\mathbf{J}\_{12} = \lambda \mathcal{F}\_0(k) b\_0 (1 - b\_0),\tag{A15b}$$

$$\mathbf{J}\_{21} = \rho v w\_0 - \lambda w\_0 \mathcal{F}\_0(k) - \lambda b\_0 w\_0 \mathcal{F}\_1(k),\tag{A15c}$$

$$\mathcal{J}\_{22} = -\nu(1 - \rho b\_0) - \frac{\lambda b\_0 \sigma\_0^2 \psi\_0}{\psi} - \delta\_w k^2. \tag{A15d}$$

Assuming *a*1, *a*<sup>2</sup> ∝ exp(*μt*) we can calculate the dispersion relation *μ*=*μ*(*k*). Substitution of the steady states into *μ* reveal their stability at a given precipitation rate (Figure A1). For *P* > *P*<sup>0</sup> the bare soil state undergoes a uniform instability (Figure A1A) that drives the system towards a uniform vegetation state, which may not necessarily be stable. For *P* < *PT* the uniform vegetation state undergoes a non-uniform instability (Figure A1B) which generates a periodic vegetation state. Both instabilities are stationary (perturbations grow monotonically in time) as *Im*(*μ*) = 0 for both of them. Furthermore, we verified that the root-augmentation feedback alone can induce a non-uniform instability by setting *δ<sup>w</sup>* = 0 (i.e., no soil–water diffusion feedback).

**Figure A1.** Dispersion relations for stable, marginally stable and unstable steady states. (**A**) Instability of the bare soil state to the growth of a uniform mode at *P*<sup>0</sup> = 175.3 mm/y (or (kg/m2)y<sup>−</sup>1). (**B**) Instability of the uniform vegetation state to the growth of a non-uniform mode (of finite wavenumber) at *PT* = 296.2 mm/y. Precipitation values: (**A**) *P* = 299 > *PT*, *P* = 293.4 < *PT*, (**B**) *P* = 175.8 > *P*0, *P* = 174.2 < *P*0. Other parameters are fixed: *E* = 7[m2/kg], *K* = 0.4[kg/m2], *M* = 10.5[1/y], *N* = 15[1/y], Λ = 0.9[(m2/kg)/y], Γ = 12[(m2/kg)/y], *R* = 0.7, *DB* = 0.1[m2/y], *DW* = <sup>50</sup>[m2/y], *<sup>S</sup>*<sup>0</sup> <sup>=</sup> 0.15m, *<sup>c</sup>*<sup>0</sup> <sup>=</sup> 1, *<sup>c</sup>*<sup>2</sup> <sup>=</sup> 0.0125, *<sup>c</sup>*<sup>4</sup> <sup>=</sup> 5.06 <sup>×</sup> <sup>10</sup><sup>−</sup>4, *<sup>c</sup>*<sup>6</sup> <sup>=</sup> <sup>−</sup>5.70 <sup>×</sup> <sup>10</sup><sup>−</sup>6, *<sup>c</sup>*<sup>8</sup> <sup>=</sup> 3.33 <sup>×</sup> <sup>10</sup><sup>−</sup>7.

#### **Appendix B. Model for Laterally Confined Roots**

Laterally confined roots are modeled by considering the limit *S*<sup>0</sup> → 0 in Equation (2). This amounts to replacing *G*(**X**, **X** , *T*) in Equation (5a) by a delta function, which leads to the following model:

$$\begin{array}{rcl} \frac{\partial B}{\partial T} &=& \Lambda B W (1 + E B)^2 (1 - \frac{B}{K}) - M B + D\_B \Delta B \,, \\\frac{\partial W}{\partial T} &=& P - \frac{N W}{1 + R B / K} - \Gamma B W (1 + E B)^2 + D\_W \Delta W \,. \\\end{array} \tag{A16}$$

This model is used in the studies described in Sections 4.2 and 4.3, with the nondimensional form

$$\begin{array}{rcl}\partial\_t b &=& bw(1+\eta b)^2(1-b) - b + \Delta b \\ \partial\_t w &=& p - \frac{nw}{1+\rho b} - \gamma bw(1+\eta b)^2 + \delta \Delta w \end{array} \tag{A17}$$

obtained by introducing the non-dimensional variables:

$$b = \frac{B}{K}, \text{ } w = \frac{W\Lambda}{M}, \text{ } \mathbf{x} = \mathbf{X}\sqrt{M/D\_{\mathbb{B}}}. \text{ } y = \mathbf{Y}\sqrt{M/D\_{\mathbb{B}}}, \text{ } \mathbf{t} = MT,\tag{A18}$$

and the non-dimensional parameters

$$p = \frac{P\Lambda}{M^2}, \; n = \frac{N}{M}, \; \gamma = \frac{\Gamma K}{M}, \; \rho = R, \; \eta = EK, \; \delta = \frac{D\_W}{D\_B}.\tag{A19}$$

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
