# Regime shifts in models of dryland vegetation

###### Abstract

Drylands are pattern-forming systems showing self-organized vegetation patchiness, multiplicity of stable states and fronts separating domains of alternative stable states. Pattern dynamics, induced by droughts or disturbances, can result in desertification shifts from patterned vegetation to bare soil. Pattern-formation theory suggests various scenarios for such dynamics; an abrupt global shift involving a fast collapse to bare soil, a gradual global shift involving the expansion and coalescence of bare-soil domains, and an incipient shift to a hybrid state consisting of stationary bare-soil domains in an otherwise periodic pattern. Using models of dryland vegetation we address the question which of these scenarios can be realized. We found that the models can be split into two groups: models that exhibit multiplicity of periodic-pattern and bare-soil states, and models that exhibit, in addition, multiplicity of hybrid states. Furthermore, in all models we could not identify parameter regimes in which bare-soil domains expand into vegetated domains. The significance of these findings is that while models belonging to the first group can only exhibit abrupt shifts, model belonging to the second group can also exhibit gradual and incipient shifts. A discussion of open problems concludes the paper.

Keywords: Models of vegetation pattern formation, multiplicity of stable states, localized patterns, fronts, homoclinic snaking, desertification.

## I Introduction

Water-limited landscapes can generally be described as mosaics of vegetation and bare-soil patches of various forms. Increasing empirical evidence supports the view that this type of vegetation patchiness is a self-organization phenomenon that would have occurred even in perfectly homogeneous physical environments Valentin1999catena (); Deblauwe2008geb (). Much insight into the mechanisms that drive self-organized vegetation patchiness has been achieved using mathematical models of water-limited landscapes Gilad2007jtb (); Rietkerk2008tree (); Borgogno2009rev_geo (). These models first demonstrate that uniform vegetation states can go through spatial instabilities to periodic vegetation patterns upon increasing environmental-stress parameters. They further highlight two main feedbacks that are capable of producing such instabilities Meron2012eco_mod (). The first is a positive feedback between biomass and water that develops as a result of an infiltration contrast between bare and vegetated areas (infiltration feedback). The second is a positive feedback between above-ground and below-ground biomass, related to the root-to-shoot ratio, a characteristic trait of any plant species (root-augmentation feedback). Model studies of vegetation pattern formation along a rainfall gradient have revealed five basic vegetation states vonHardenberg2001prl (); Rietkerk2002an (); Lejeune2004ijqc (); Gilad2004prl (): uniform vegetation, gap patterns, stripe (labyrinth) patterns, spot patterns and uniform bare-soil. Another significant result is the existence of precipitation ranges where alternative stable vegetation states coexist. These are generally bistability ranges of any consecutive pair of basic states: bare-soil and spots, spots and stripes, stripes and gaps and gaps and uniform vegetation. Within any bistability range, spatial mixtures of the two alternative stable states can form long transient patterns that culminate in one of the two alternative states, or stable asymptotic hybrid patterns Meron2012eco_mod ().

The mathematical theory of hybrid patterns is far from being complete. Much progress, however, has been made for the simpler case of bistability of uniform and periodic-pattern states, using simple pattern formation models such as the Swift-Hohenberg equation Knobloch2008nonlinearity (). According to this theory a bistability range of uniform and patterned states may contain a subrange (or an overlapping range) of stable localized patterns, coexisting with the two alternative stable states. For bistability of bare-soil and vegetation spot patterns these localized patterns would correspond to isolated spot-pattern domains in an otherwise bare-soil, or conversely, to isolated bare-soil domains in an otherwise periodic spot pattern. The appearance of these mixed-pattern or hybrid states can be understood intuitively by focusing on the dynamics of the transition zones that separate the two alternative stable states. These zones,
are fronts that can be stationary or propagating. In the case of bistability of two uniform states, isolated fronts always propagate, except for a singular control-parameter value, the so called Maxwell point, at which the direction of propagation changes Bel2012theo_ecol ()^{1}^{1}1A pair of fronts
propagating towards one another, however, can slow down and become stationary due to repulsive front interactions Hagberg1994nonlinearity ().. Bistability of uniform and pattern states, on the other hand, allow for an additional behavior; isolated fronts can be stationary or pinned in a range of the control parameter Pomeau1986pd (). Such a range can give rise to many hybrid states, because the fronts that constitute the boundaries of the alternative-state domains are stationary. In a diagram that shows the various states as functions of the control parameter, the hybrid states often appear as solution branches that ”snake” down from the periodic-pattern branch towards the uniform (zero) state as Fig. 1 illustrates. The control-parameter range where these solutions exists is often called the snaking range and the appearance of such solutions is described as homoclinic snaking Knobloch2008nonlinearity (). In the following we will refer to this range as the ”hybrid-state range” to allow for multistability of hybrid states that is not associated with homoclinic snaking.

Bistability of alternative stable states has been studied extensively in the context of ecosystem regime shifts, i.e. sudden transitions to a contrasting state in response to gradual changes in environmental conditions Scheffer2001nature (); Scheffer2003tree (). Such shifts have been observed in lakes, coral reefs, oceans, forests and arid lands. Global shifts from one stable state to another, however, may not necessarily be abrupt. Ecosystems are continuously subjected to local disturbances whose effects are spatially limited. Examples of such disturbances in the context of water-limited vegetation include clear cutting, grazing, infestation and limited fires. These disturbances can induce fast local transitions to the alternative stable state, but, according to pattern formation theory, the subsequent dynamics may proceed slowly by the expansion and coalescence of the domains of the alternative stable state through front propagation and front collisions. Such a succession of processes eventually leads to a global regime shift, but in a gradual manner Bel2012theo_ecol ().

How slow can gradual shifts be? When the two alternative stable states are spatially uniform the pace of a gradual shift depends on the value of the control parameter relative to the Maxwell point; the larger the distance from that point the faster the gradual shift. This result often holds for bistability of uniform and patterned states too, except for one important difference - the value of the control parameter should be outside the hybrid-state range (but still within the bistability range) Pomeau1986pd (); Knobloch2008nonlinearity (). The difference between abrupt and gradual shifts can be dramatic, as Fig. 2 illustrates. For systems whose spatial extent is much larger than the size of a spot, gradual shifts can occur on time scales that are orders of magnitude longer than those of abrupt shifts.

Within the hybrid-state range global regime shifts are not expected to occur in steady environments. The system rather shows spatial plasticity; any spatial disturbance pattern shifts the system to the closest hybrid pattern, which is a stable stationary state and therefore involves no further dynamics. It is worth noting that transitions from the periodic pattern state to hybrid patterns, within the hybrid-state range, can also occur as a result of global uniform environmental changes, such as a precipitation drop or a uniform disturbance, provided the initial pattern is not perfectly periodic, e.g. hexagonal spot pattern containing penta-hepta defects Meron2012eco_mod ().

Bistability of uniform and patterned states is most relevant to desertification, a regime shift involving a transition from a productive vegetation-pattern state to an unproductive uniform bare-soil state vonHardenberg2001prl (); Rietkerk2004science (). To what extent are the general results of pattern formation theory displayed in Figs. 1 and 2 applicable to the specific context of desertification? We address this question by studying specific models of vegetation pattern formation of various degrees of complexity. The manuscript is organized as follows. In section II we briefly review the models of water-limited vegetation considered here. In section III we present numerical results for these models, distinguishing between models for which we found indications for hybrid states (homoclinic snaking) and models for which we have not found such indications. These results are discussed and summarized in section IV.

## Ii Models for spatial vegetation dynamics

We chose to study several representative models of increasing complexity. All models are deterministic and specifically constructed to describe vegetation patchiness in water-limited flat terrains (unlike the variant of the Swift-Hohenberg equation used to produce Figs. 1 and 2). The degree of complexity is reflected by the number of dynamical variables and by the number of pattern-forming feedbacks the model captures. The models consist of partial differential equations (PDEs) for a continuous biomass variable and possibly for additional water variables, depending on the model. All models capture an instability of a uniform vegetation state to a periodic-pattern state and a bistability range of periodic patterns and bare soil.

### ii.1 Lefever-Lejeune (LL) model

The simplest model we consider is a single-variable model for a vegetation biomass density, , introduced by Lefever and Lejune Lefever1997bmb (); Lejeune1999 (). We chose to study a simplified version of this model Lejeune2002pre (); Lefever2000 (); Tlidi2008lnp () whose form in terms of non-dimensional variables and parameters is

(1) |

In this equation the parameter is the mortality to growth ratio, is the degree of facilitative, relative to competitive, local interactions experienced by the plants, and is the ratio between the spatial range of facilitative interactions and the range of competitive interactions. The spatial derivative terms represent short range facilitation and long range competition, a well known pattern formation mechanism Gierer1972kibernetic (). The agents responsible for this mechanism in actual dryland landscapes are nonlocal feedbacks involving water transport towards growing vegetation patches. Explicit modeling of these feedbacks requires the addition of water variables. Although the model does not include a precipitation parameter, water stress can be accounted for by increasing the mortality parameter . In what follows we will refer to this model as the LL model.

### ii.2 Modified Klausmeier (K) model

Next in degree of complexity is a modified version vanderStelt2012nonl_sci () of a model introduced by Klausmeier Klausmeier1999science (), hereafter the K model. In addition to a biomass density variable, , this model contains a water variable, , which we regard as representing soil-water content. The model equations, expressed in terms of non-dimensional quantities, are

(2a) | ||||

(2b) |

where . According to equation (2a) the biomass growth rate, , increases with the biomass density, reflecting a positive local facilitation feedback. Natural mortality at a rate , acts to reduce the biomass, and local seed dispersal or clonal growth, represented by the diffusion term , act to distribute the biomass to adjacent areas. The water dynamics (equation (2b)) is affected by precipitation with a rate , evaporation and drainage (), biomass-dependent water-uptake rate (), and by soil-water diffusion. The pattern-forming feedback in this model is induced by the combined effect of higher local water-uptake rate in denser vegetation patches and fast water diffusion towards these patches, which inhibits the growth in the patch surroundings. This mechanism may be applicable to sandy soils for which is relatively large. This third type of pattern-forming feedback (besides the infiltration and root-augmentation feedbacks) has not been stressed in earlier studies.

The original Klausmeier model Klausmeier1999science () does not include a water diffusion term, but rather an advection term to describe runoff on a slope. While accounting for banded vegetation on a slope, the original model does not produce stationary vegetation patterns in flat terrains. To capture the latter we added the soil-water diffusion term vanderStelt2012nonl_sci (). Since we focus on plane terrains we do not need an advection term and therefore omitted it.

### ii.3 Rietkerk et al. (R) model

The third model we consider, the R model, distinguishes between below-ground and above-ground water dynamics by introducing two water variables; , representing soil water, and , representing surface water. This three-variable model has been introduced by Rietkerk et al. Rietkerk2002an (); HilleRisLambers2001ecology () and consists of the following non-dimensional equations ^{2}^{2}2The non-dimensional form of the equations was derived from the dimensional model, presented in Rietkerk2002an (), using the following transformation (the denotes the variables in Rietkerk2002an ()):
, , , , , , , , , , , , where , .:

(3a) | ||||

(3b) | ||||

(3c) |

where

(4) |

In equation (3a) the biomass growth rate, , depends on the soil-water variable only (no biomass dependence as in the K model); the dependence is linear at small soil-water contents and approaches a constant value at high contents, representing full plant turgor. Biomass growth is also affected by mortality () and by seed dispersal or clonal growth (). Soil-water content (equation (3b)) is increased by the infiltration of surface water (). The biomass dependence of the infiltration rate, , captures the infiltration contrast that exists between bare soil (low infiltration rate) and vegetated soil (high infiltration rate) for . The other terms affecting the dynamics of the soil water represent loss of water due to evaporation and drainage (), water uptake by the plants (), and moisture diffusion within the soil. The dynamics of surface water (equation (3c)) are affected by precipitation at a rate , by water infiltration into the soil, and by overland flow modeled as a diffusion process.

The R model captures an important pattern-forming feedback - the infiltration feedback. When the infiltration contrast is high () patches with growing vegetation act as sinks for runoff water. This accelerates the vegetation growth, sharpens the infiltration contrast and increase even further the soil water content in the patch areas. The water flow towards vegetation patches inhibits the growth in the patch surroundings thereby promoting vegetation pattern formation. The infiltration feedback allows vegetation pattern formation at lower, more realistic, values of the soil-water diffusion constant in comparsion to the K model.

### ii.4 Gilad et al. (G) model

The fourth model to be studied, the G model, was introduced by Gilad et al. Gilad2004prl (); Gilad2007jtb () and contains the same three dynamical variables, , and as the R model, but with the interpretation of as representing the above-ground biomass. This is because the G model explicitly considers the root system and the relation between the root-zone size and the above-ground biomass. This additional element allows the introduction of another important pattern-forming feedback besides the infiltration feedback, the root-augmentation feedback. The model equations, in non-dimensional forms, read

(5a) | ||||

(5b) | ||||

(5c) |

Like in the K model, the biomass growth rate, , depends both on and but in a non-local way that accounts for the contribution of soil-water availability at point to biomass growth at point through a biomass-dependent root system that extends from point to point . Similarly, the water-uptake rate, , by the plants’ roots depends on and in a nonlocal manner to account for the uptake at a point by a plant located at whose roots extend to . Specifically,

(6a) | ||||

(6b) | ||||

(6c) |

The root-augmentation feedback is captured by letting the width of the root kernel , which represents the lateral root-zone size, to linearly increase with the above-ground biomass. As a plant grows its root zone extends to new soil regions. As a result the amount of water available to the plant increases and the plant can grow even further. While accelerating the local plant growth, this process also depletes the soil-water content in the plant surroundings, thereby inhibiting the growth there and promoting vegetation pattern formation. The proportionality parameter appearing in equation (6c) controls the strength of the root-augmentation feedback. It is a measure of the root-to-shoot ratio, a characteristic plant trait. Note that the soil-water dependence of the biomass growth term in equation (5a) and of the water uptake term in equation (5b) is linear. Nonlinear forms, including that used in the R model, have been studied in Ref. Kletter2009jtb (). Like in the R model, the infiltration feedback appears through the biomass-dependent form of the infiltration rate :

(7) |

Other differences with respect to the R model involve the introduction of (i) the logistic growth form in equation (5a), which represents genetic growth limitations at high biomass densities (e.g. stem strength), (ii) the biomass-dependent evaporation rate in the soil-water equation (5b) (second term on right side) which accounts for reduced evaporation by canopy shading and introduces a local positive water-biomass feedback, and (iii) nonlinear overland flow term in the surface-water equation (5c) motivated by shallow water theory Gilad2004prl (); santillana2010comp-geo (), rather than a diffusion term as in the R model.

### ii.5 Simplified Gilad et al. (SG) model

The fifth model is a simplified version of the G model, in which the root kernel is assumed to vary sharply in comparison to and , and therefore can be approximated by a Dirac delta function. This approximation is suitable for plant species that grow deep roots with small lateral dimensions. The simplified model, denoted SG, reads

(8a) | ||||

(8b) | ||||

(8c) |

This version of the model includes the same pattern-forming infiltration feedback as the original model ( is defined the same way it was defined in the G model), but the root-augmentation feedback is modified; water transport towards growing vegetation patches is no longer a result of uptake by the laterally spread roots, but rather a result of soil-water diffusion.

## Iii Results of numerical model studies

The ecological context we consider is water-limited ecosystems in flat terrains exhibiting bistability of a periodic vegetation pattern and bare soil. We will mostly be concerned with initial states consisting of periodic patterns that are locally disturbed to form bare-soil domains. The numerical studies described below are based on numerical continuation methods, used to identify spatially periodic solutions, and on PDE solvers, used to identify stable branches of localized patterns and to follow the dynamics of bare-soil domains. As we will shortly argue, these dynamics crucially depend on the additional stable pattern states, periodic or localized, that the system supports.

There are several properties that all models appear to share: (i) the coexistence of a family of stable periodic solutions, describing vegetation patterns of different wavelengths, with a stable uniform solution that describes the bare-soil state; (ii) bare-soil domains do not expand into patterned domains; (iii) the existence of a stable localized solution describing a single vegetation spot in an otherwise bare soil state. An additional property that is most significant for regime shifts is not shared by all models - multiplicity of stable hybrid states. We use this property to divide the models into two groups, models that do not show multiplicity of stable hybrid states and models that do show such a multiplicity of states. The two groups display different forms of regime shifts as described below.

### iii.1 Models lacking multiplicity of hybrid states

The models that belong to this group are the K model (Section II.2), the R model (Section II.3) and the SG model (Section II.5). These models have wide bands of periodic solutions with stable branches that coexist with the stable branch of the bare-soil solution Dijkstra2011IJBC (); vanderStelt2012nonl_sci (). Figure 3 shows bifurcation diagrams for the R and SG models in 1D, computed by a numerical continuation method auto (). The bifurcation parameter was chosen to be the precipitation rate . The diagrams show overlapping periodic solutions whose wavelengths increase as decreases. The last periodic solution to exist corresponds to a single hump. We have not been able to identify (by numerical continuation) solution branches that describe hybrid patterns, either groups of humps in an otherwise bare soil state or holes in an otherwise periodic pattern. To further test whether such solutions can exist in these models or, if they exist, whether they are stable, we solved the models’ equations numerically using initial conditions that describe fronts separating the patterned and the bare-soil states. Convergence to front solutions that are stationary over a range of values would indicate the possible existence of hybrid solutions Pomeau1986pd (); Knobloch2008nonlinearity (). Such front pinning, however, has not been observed; in all simulations the patterned state propagated into the bare-soil state. We conclude that stable hybrid solutions, apart from a single hump solution, do not exist in these models, or if they do, their existence range is extremely small.

In order to study regime shifts in the K, R and SG models we simulated the model equations within the bistability range of periodic patterns and bare soil, starting with periodic patterns that contain bare-soil domains. Since the patterned state was always found to propagate into the bare-soil state, such initial bare-soil domains contract and disappear. This behavior rules out the occurrence of a gradual regime shift to the bare-soil state (similar to that shown in panels e-h of Fig. 2). The final pattern, however, can differ from the initial one in its wavelength as the 1D simulations of the R model displayed in Fig. 4 show. The system can respond by mere readjustment of the spacings between individual humps without a change in their number, which leads to an increase in the pattern’s wavelength (left panel), or, at higher precipitation, by hump splitting, which results in a decrease of the pattern’s wavelength (right panel). Similar responses to local disturbances were found in the K and SG models. Figure 5 displays results of 2D simulations of the SG model showing that the two response forms, spacing readjustments and spot splitting, can occur at the same precipitation by changing the size of the initial bare-soil domain. Reducing the precipitation rate to values below the bistability range of periodic patterns and bare soil leads to an abrupt global transition to the bare-soil state as Fig. 6 shows.

### iii.2 Models exhibiting multiplicity of hybrid states

Numerical solutions of the LL and G models (Sections II.1 and II.4) using PDE solvers point towards the existence of stable hybrid states in addition to periodic-pattern states ^{3}^{3}3The application of numerical continuation methods for these models is harder in comparison to the K, R and SG models and has not been pursued in this study..
Fig. 7 shows a bifurcation diagram for the LL model, using the mortality rate as the bifurcation parameter. The upper solution branch corresponds to a periodic-pattern state, while the lowest branch corresponds to the bare-soil state. The red branches in between correspond to stable hybrid states describing localized patterns, a few examples of which are shown in the right panels. Solutions of this kind in 1D and 2D have been found earlier Lejeune2002pre (). Figure 8 shows a partial bifurcation diagram for the G model in 2D. The upper line corresponds to a spot-pattern state ^{4}^{4}4The spot pattern is rhombic rather than hexagonal. Such patterns apparently exist in the G model, like in other pattern-formation models Bachir2001epl (). In the present case, the system converged to the rhombic pattern following a disturbance of an hexagonal pattern. Since it contains a direction of denser spots (vertical) along which the water stress is higher,
transitions to hybrid states can be induced upon decreasing in the absence of local disturbances., while the lower lines correspond to hybrid patterns with decreasing number of spots as the right panels show. Note the difference between the hybrid solution branches in the two models; while in the LL model they all terminate at the same control-parameter value , which coincides with the fold-bifurcation point of the periodic pattern solution, in the G model the hybrid solution branches are slanted Dawes2008sjads () - solutions with smaller numbers of spots terminate at lower values.

The multitude of stable hybrid patterns, i.e. patterns consisting of groups of spots in an otherwise bare soil, groups of holes in otherwise periodic patterns and various combinations thereof, suggests a form of spatial plasticity. That is, any pattern of local disturbances shifts the system to the closest hybrid pattern with no further dynamics. This behavior rules out the occurrence of a gradual regime shift as a result of initial local disturbances, but unlike the K, R and SG models the system does not recover from the disturbances. This suggests the possible occurrence of a gradual regime shift in a continuously disturbed system.

While the two models share spatial plasticity in response to local disturbances, they differ in the response to gradual parameter changes ( or ). In the LL model all localized pattern solutions terminate at the fold bifurcation point (see Fig. 7). Above that point the only stable state is bare soil and, therefore, any hybrid state must collapse to this state. Note the difference between the bifurcation diagram in Fig. 7 and the diagram obtained with the Swift-Hohenberg equation in Fig. 1. In the latter there is a subrange () outside the hybrid-state range which is still within the bistability range, where disturbed patterns go through gradual shifts. No such subrange has been found in the LL model. Contrary to the LL model, the slanted structure of localized pattern solutions in the G model, allows for a gradual response. In fact, the hybrid state (b) in Fig. 8 was obtained from the periodic state (a) by an incremental decrease of . Likewise, the hybrid states (c) and (d) were obtained from the states (b) and (c) by further incremental decreases of . The degree of slanting increases as the root-to-shoot parameter is increased.

## Iv Discussion

All models considered in this study predict the same basic vegetation states and stability properties along a rainfall gradient, including a bistability range of bare soil and periodic spot patterns. We may therefore expect these models to depict similar scenarios for desertification shifts, i.e. transitions from productive spot patterns to the unproductive bare-soil state. Pattern-formation theory,
represented here by results obtained with the Swift-Hohenberg equation^{5}^{5}5We refer here to the Swift-Hohenberg equation as a prototype of pattern-formation behaviors in systems showing bistability of uniform and patterned states, and to Ref. Bel2012theo_ecol () for the implications of these general behaviors to the context of regime shifts. Similar pattern-formation behaviors have been found with many other pattern-formation models Knobloch2008nonlinearity ()., suggests various possible forms for such scenarios; abrupt, gradual or incipient, induced by environmental changes, by disturbances or both. Underlying these forms are several nonlinear behaviors. The first and simplest is a global transition from the spot-pattern state to the bare-soil state, induced by a slow change of a control parameter past a fold bifurcation, or by a disturbance that shifts the system as a whole to the attraction basin of the bare-soil state. Such processes induce global abrupt shifts to the bare-soil state as Fig. 2(a-d) illustrates. Local disturbances, on the other hand, can lead to partial shifts that result in spatially-limited domains of the bare-soil state in an otherwise periodic-pattern state. The subsequent course of events depends on the dynamics of the fronts that bound these domains. When the fronts propagate, a slow process of expansion and coalescence of bare-soil domains can eventually culminate in a global gradual shift, as Fig. 2(e-h) illustrates. When the fronts are pinned, the domains remain fixed in size, after some small adjustments, in which case the shift is incomplete or incipient - the system converges to one of the many hybrid states it supports. To our surprise the models we studied do not capture all possible scenarios pattern-formation theory allows. Moreover, scenarios that are captured by some models are not captured by others.

Our studies first suggest that in all five vegetation models (K, R, SG, LL, G) the bare-soil state never grows at the expense of the periodic-pattern state (unlike the behavior shown in Fig. 2(e-h)) through the entire bistability range; bare-soil domains either stay fixed in size or contract and disappear. Furthermore, the K, R, and SG models do not show hybrid states at all, while the models that do show hybrid states, LL and G, differ in the existence ranges of these states. In the LL model the branches of all hybrid states terminate at the same threshold which coincides with that of the periodic pattern state, while in the G model the termination points are aligned on a slanted line. The results for the K, R and SG models suggest that shifts to the bare-soil state can only occur outside the bistability range of vegetation patterns and bare soil, and are therefore abrupt. Within the bistability range, bare-soil domains induced by local disturbances contract and disappear, thus restoring the vegetation-pattern state, although a wavelength change is likely to occur. Both the LL and G models predict the possible occurrence of incipient regime shifts within the bistability range of periodic vegetation patterns and bare soil. These shifts can be induced by local-disturbance regimes and culminate in one of the stable hybrid states when the disturbance regimes are over. Complete shifts to the bare-soil state, due to increased stress, are abrupt in the LL model but can be gradual in the G model because of the slanted structure of the hybrid solution branches; incremental precipitation decrease in the G model can result in step-like transitions to hybrid states of lower bioproductivity as Fig. 8 shows.

These results raise several open questions. The first is related to the finding that bare-soil domains do not expand into patterned domains in the entire bistability range. This behavior can be attributed to the positive pattern-forming infiltration and root-augmentaiton feedbacks. Both give advantage to plants at the rim of a patterned domain as compared with inner plants; the rim plants receive more runoff from the surrounding bare soil and experience weaker competition for soil water. These factors act against the retreat of vegetated domains. Processes that may favor such a retreat include soil erosion and roots exposure in sandy soils under conditions of high wind power Okin2001jae (), or insect outbreak Allen2010fem (). Whether bare-soil expansion can be explained by water-biomass interactions alone, or additional processes must be considered, is still an open question that calls for both empirical and further model studies. From the perspective of pattern formation theory the finding that the bare-soil state never expands into vegetation-pattern states questions the utility of the Maxwell-point concept far from the instability of uniform vegetation to periodic patterns and calls for further mathematical analysis.

Another open question is what elements in the LL and G models, and correspondingly what ecological and physical processes, are responsible for the multitude of stable hybrid states. The results for the LL model clearly show that reducing local facilitation, by decreasing the parameter , narrows down the hybrid-state range and can eliminate the hybrid states altogether. However, it also narrows down the bistability range of periodic patterns and bare soil, and therefore does not resolve processes that favor the formation of localized patterns alone. The results for the G model and its simplified version SG hint towards the possible role of the nonlocal water uptake by laterally extended root systems in inducing hybrid states. This nonlocal competition mechanism is absent in the SG model and may possibly be responsible for the absence of hybrid states in this model. Further studies are needed, first to substantiate the existence of stable hybrid states, particularly in the G model, and second to clarify the roles of local and nonlocal facilitation and competition processes in inducing them.

Finally, the models we have studied are all deterministic. Real ecosystems, however, are generally subjected to stochastic fluctuations in time and space, which may affect the bifurcation structure of spatial states. Additive temporal noise, for example, can induce the propagation of pinned fronts Clerc2005prl (), and thereby affect the hybrid-state range. The effect of noise on abrupt, gradual and incipient regime shifts is yet another open problem that calls for further studies.

Studying these questions is significant for identifying the nature of desertification shifts, i.e. whether they are abrupt, gradual or incipient, in various environments and for different plant species, and for assessing the applicability of early warning signals for imminent desertification Scheffer2009nature ().

###### Acknowledgements.

We wish to thank Arjen Doelman and Moshe Shachak for helpful discussions and Arik Yochelis for helping with the numerical continuation analysis. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant number [293825].## References

- (1) C. Valentine, J.M. d’Herbes, and J. Poesen. Soil and water components of banded vegetation patterns. Catena, 37:1–24, 1999.
- (2) V. Deblauwe, N. Barbier, P. Couteron, O. Lejeune, and J. Bogaert. Annual plant-shrub interactions along an aridity gradient. Global Ecology and Biogeography, 17:715–723, 2008.
- (3) E. Gilad, J. Von Hardenberg, A. Provenzale, M. Shachak, and E. Meron. A Mathematical Model for Plants as Ecosystem Engineers. J. Theor. Biol., 244:680, 2007.
- (4) Max Rietkerk and Johan van de Koppel. Regular pattern formation in real ecosystems. Trends in Ecology & Evolution, 23(3):169–175, 2008.
- (5) F. Borgogno, P. D’Odorico, F. Laio, and F. Laio. Mathematical models of vegetation pattern formation in ecohydrology. Rev. Geophys., 47:RG1005, 2009.
- (6) E. Meron. Pattern-formation approach to modelling spatially extended ecosystems. Ecological Modelling, 234:70–82, 2012.
- (7) J. Von Hardenberg, E. Meron, M. Shachak, and Y. Zarmi. Diversity of Vegitation Patterns and Desertification. Phys. Rev. Lett., 89(198101), 2001.
- (8) M. Rietkerk, M.C. Boerlijst, F. van Langevelde, R. HilleRisLambers, J. van de Koppel, L. Kumar, H.H.T. Prins, and A.M. De Roos. Self-organization of vegetation in arid ecosystems. Am. Nat., 160:524–530, 2002.
- (9) O. Lejeune, M. Tlidi, and R. Lefever. Vegetation spots and stripes: dissipative structures in arid landscapes. Int. J. Quant. Chem., 98:261–271, 2004.
- (10) E. Gilad, J. Von Hardenberg, A. Provenzale, M. Shachak, and E. Meron. Ecosystem Engineers: From Pattern Formation to Habitat Creation. Phys. Rev. Lett., 93(098105), 2004.
- (11) E. Knobloch. Spatially localized structures in dissipative systems: open problems. Nonlinearity, 21:T45, 2008.
- (12) Golan Bel, Aric Hagberg, and Ehud Meron. Gradual regime shifts in spatially extended ecosystems. Theoretical Ecology, 5:591–604, 2012.
- (13) Y. Pomeau. Front motion, metastability and subcritical bifurcations in hydrodynamics. Physica D, 23:3, 1986.
- (14) M. Scheffer, S. Carpenter, J.A. Foley, C. Folke, and B. Walkerk. Catastrophic shifts in ecosystems. Nature, 413:591–596, 2001.
- (15) M. Scheffer and S. R. Carpenter. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol., 18:648–656, 2003.
- (16) M. Rietkerk, S. C. Dekker, P. C. de Ruiter, and J. van de Koppel. Self-Organized Patchiness and Catastrophic Shifts in Ecosystems. Science, 305:1926–1929, 2004.
- (17) R. Lefever and O. Lejeune. On the origin of tiger bush. Bull. Math. Biol., 59:263–294, 1997.
- (18) O. Lejeune. Une théorie champ moyen de l’organisation spatio-temporelle des écosystèmes végétaux. PhD thesis, University of Brussels, 1999.
- (19) O. Lejeune, M. Tlidi, and P. Couteron. Localized vegetation patches: A self-organized response to resource scarcity. Phys. Rev. E, 66:010901, Jul 2002.
- (20) R. Lefever, O. Lejeune, and P. Couteron. Mathematical Models for Biological Pattern Formation, volume IMA Volumes in Mathematics and its Applications, vol 121, chapter Generic modelling of vevetation patterns, pages 83–112. Springer, 2000.
- (21) M. Tlidi, R. Lefever, and A. Vladimirov. On Vegetation Clustering, Localized Bare Soil Spots and Fairy Circles. Lect. Notes Phys., 751:381–402, 2008.
- (22) A. Gierer and H. Meinhardt. A theory of biological pattern formation. Kibernetic, 12:30–39, 1972.
- (23) S. van der Stelt, A. Doelman, G.M. Hek, and J. Rademacher. Rise and fall of periodic patterns for a generalized Klausmeier-Gray-Scott model. J. Nonl. Sc., 2012.
- (24) C.A. Klausmeier. Regular and irregular patterns in semiarid vegetation. Science, 284:1826–1828, 1999.
- (25) R. HilleRisLambers, M. Rietkerk, F. Van den Bosch, H.H.T. Prins, and H. de Kroon. Vegetation pattern formation in semi-arid grazing systems. Ecology, 82:50–61, 2001.
- (26) A. Kletter, J. Von Hardenberg, E. Meron, and A. Provenzale. Patterned Vegetation and Rainfall Intermittency. J. Theo. Biol., 256:574–583, 2009.
- (27) Mauricio Santillana and Clint Dawson. A numerical approach to study the properties of solutions of the diffusive wave approximation of the shallow water equations. Computational Geosciences, 14:31–53, 2010.
- (28) H. A. Dijkstra. Vegetation pattern formation in a semid-arid climate. International Journal of Bifurcation and Chaos, 21(12):3497–3509, 2011.
- (29) E. Doedel, R. C. Paffenroth, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. E. Oldeman, B. Sandstede, and X. Wang. AUTO2000, Technical report, Concordia University, 2002. Technical report.
- (30) J. Dawes. Localized pattern formation with a large-scale mode: Slanted snaking. SIAM Journal on Applied Dynamical Systems, 7(1):186–206, 2008.
- (31) Gregory S Okin, Bruce Murray, and William H Schlesinger. Degradation of sandy arid shrubland environments: observations, process modelling, and management implications. Journal of Arid Environments, 47(2):123 – 144, 2001.
- (32) Craig D. Allen, Alison K. Macalady, Haroun Chenchouni, Dominique Bachelet, Nate McDowell, Michel Vennetier, Thomas Kitzberger, Andreas Rigling, David D. Breshears, E.H. (Ted) Hogg, Patrick Gonzalez, Rod Fensham, Zhen Zhang, Jorge Castro, Natalia Demidova, Jong-Hwan Lim, Gillian Allard, Steven W. Running, Akkin Semerci, and Neil Cobb. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. Forest Ecology and Management, 259(4):660 – 684, 2010.
- (33) M. G. Clerc, C. Falcon, and E. Tirapegui. Additive noise induces front propagation. Phys. Rev. Lett., 94:148302, Apr 2005.
- (34) M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. van Nes, M. Rietkerk, and G. Sugihara. Early-warning signals for critical transitions. Nature, 461:387–393, 2009.
- (35) A. Hagberg and E. Meron. Pattern Formation in Non-Gradient Reaction Diffusion Systems: The Effects of Front Bifurcations. Nonlinearity, 7:805–835, 1994.
- (36) M. Bachir, S. Métens, P. Borckmans, and G. Dewel. Formation of rhombic and superlattice patterns in bistable systems. Europhys. Lett., 54:612–618, 2001.