## 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 of 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, models belonging to the second group can also exhibit gradual and incipient shifts. A discussion of open problems concludes the paper.

## 1. 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 [1,2]. Many insights into the mechanisms that drive self-organized vegetation patchiness have been achieved using mathematical models of water-limited landscapes [3–5]. 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 [6]. The first is also 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 [7–10]: 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 [6].

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 [11]. 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 [12].^{1} Bistability of uniform and pattern states, on the other hand, allows for an additional behaviour: isolated fronts can be stationary or pinned in a *range* of the control parameter [14]. 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 figure 1 illustrates. The control parameter range where these solutions exist is often called the snaking range, and the appearance of such solutions is described as homoclinic snaking [11]. In the following, we will refer to this range as the ‘hybrid-state range’ to allow for multi-stability 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 [15,16]. 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 [12].

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) [14,11]. The difference between abrupt and gradual shifts can be dramatic, as figure 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 [6].

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 [7,17]. To what extent are the general results of pattern formation theory displayed in figures 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 paper is organized as follows. In §2, we briefly review the models of water-limited vegetation considered here. In §3, 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 §4.

## 2. 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 figures 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.

### (a) Lefever–Lejeune model

The simplest model we consider is a single-variable model for a vegetation biomass density, *b*(**r**,*t*), introduced by Lefever & Lejeune [18,19]. We chose to study a simplified version of this model [20–22] whose form in terms of non-dimensional variables and parameters is
2.1In equation (2.1), the parameter *μ* is the mortality to growth ratio, *Λ* is the degree of facilitative, relative to competitive, local interactions experienced by the plants and *L* 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 [23]. The agents responsible for this mechanism in actual dryland landscapes are non-local feedbacks involving water transport towards growing vegetation patches. Explicit modelling 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 Lefever–Lejeune (LL) model.

### (b) Modified Klausmeier model

Next in degree of complexity is a modified version [24] of a model introduced by Klausmeier [25]; hereafter, the K model. In addition to a biomass density variable, *b*, this model contains a water variable, *w*, which we regard as representing soil-water content. The model equations, expressed in terms of non-dimensional quantities, are
2.2aand
2.2bwhere *G*=*wb*. According to equation (2.2a), the biomass growth rate, *G*, 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 ∇^{2}*b*, acts to distribute the biomass to adjacent areas. The water dynamics (equation (2.2b)) is affected by precipitation with a rate *p*, evaporation and drainage (−*w*), the biomass-dependent water uptake rate (−*Gb*=−*b*^{2}*w*) and by soil-water diffusion. The pattern-forming feedback in this model is induced by the combined effect of a 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 *D*_{w} 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 K model [25] does not include a water diffusion term, but rather an advection term to describe run-off 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 [24]. As we focus on plane terrains, we do not need an advection term and therefore omitted it.

### (c) Rietkerk *et al.* model

The third model we consider, the Rietkerk *et al.* (R) model, distinguishes between below-ground and above-ground water dynamics by introducing two water variables: *w*, representing soil water, and *h*, representing surface water. This three-variable model was introduced by Rietkerk and co-workers [8,26] and consists of the following non-dimensional equations:^{2}
2.3a
2.3b
and
2.3cwhere
2.4In equation (2.3a), the biomass growth rate, *G*=*G*(*w*), 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 (−*μb*) and by seed dispersal or clonal growth (∇^{2}*b*). Soil-water content (equation (2.3b)) is increased by the infiltration of surface water (*Ih*). The biomass dependence of the infiltration rate, *I*=*I*(*b*), captures the infiltration contrast that exists between bare soil (low infiltration rate) and vegetated soil (high infiltration rate) for *f*<1. The other terms affecting the dynamics of the soil water represent loss of water owing to evaporation and drainage (−*νw*), water uptake by the plants (−*γGb*) and moisture diffusion within the soil. The dynamics of surface water (equation (2.3c)) are affected by precipitation at a rate *p*, by water infiltration into the soil and by overland flow modelled as a diffusion process.

The R model captures an important pattern-forming feedback—the infiltration feedback. When the infiltration contrast is high ( *f*≪1), patches with growing vegetation act as sinks for run-off water. This accelerates the vegetation growth, sharpens the infiltration contrast and increases 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 comparison with the K model.

### (d) Gilad *et al.* model

The fourth model to be studied, the Gilad *et al.* (G) model, was introduced by Gilad *et al.* [10,3] and contains the same three dynamical variables—*b*, *w* and *h*—as the R model, but with the interpretation of *b* 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
2.5a
2.5b
and
2.5cAs in the K model, the biomass growth rate, *G*_{b}, depends on both *w* and *b* but in a non-local way that accounts for the contribution of soil-water availability at point **x**′ to biomass growth at point **x** through a biomass-dependent root system that extends from point **x** to point **x**′. Similarly, the water uptake rate, *G*_{w}, by the plants’ roots depends on *b* and *w* in a non-local manner to account for the uptake at a point **x** by a plant located at **x**′ whose roots extend to **x**. Specifically,
2.6a
2.6b
and
2.6cThe root-augmentation feedback is captured by allowing the width of the root kernel *g*, 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 (2.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 (2.5a) and of the water uptake term in equation (2.5b) is linear. Nonlinear forms, including that used in the R model, have been studied in [27]. As in the R model, the infiltration feedback appears through the biomass-dependent form of the infiltration rate *I*,
2.7Other differences with respect to the R model involve the introduction of (i) the logistic growth form *b*(1−*b*) in equation (2.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 (2.5b) (second term on the right-hand side), which accounts for reduced evaporation by canopy shading and introduces a local positive water–biomass feedback, and (iii) the nonlinear overland flow term in the surface-water equation (2.5c) motivated by shallow water theory [10,28], rather than a diffusion term as in the R model.

### (e) Simplified Gilad *et al.* model

The fifth model is a simplified version of the G model, in which the root kernel *g* is assumed to vary sharply in comparison with *b* and *w*, 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 simplified Gilad (SG), reads
2.8a
2.8b
and
2.8cThis version of the model includes the same pattern-forming infiltration feedback as the original model (*I* is defined in 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.

## 3. 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 (WLs), with a stable uniform solution that describes the bare-soil state; (ii) bare-soil domains do not expand into patterned domains; and (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.

### (a) Models lacking multiplicity of hybrid states

The models that belong to this group are the K model (§2*b*), the R model (§2*c*) and the SG model (§2*e*). These models have wide bands of periodic solutions with stable branches that coexist with the stable branch of the bare-soil solution [24,29]. Figure 3 shows bifurcation diagrams for the R and SG models in one dimension, computed by a numerical continuation method [30]. The bifurcation parameter was chosen to be the precipitation rate *p*. The diagrams show overlapping periodic solutions whose WLs increase as *p* 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 *p* values would indicate the possible existence of hybrid solutions [11,14]. 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. As the patterned state was always found to propagate into the bare-soil state, such initial bare-soil domains contract and disappear. This behaviour rules out the occurrence of a gradual regime shift to the bare-soil state (similar to that shown in figure 2*f*–*h*). The final pattern, however, can differ from the initial one in its WL as the one-dimensional simulations of the R model displayed in figure 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 WL (*a*), or, at higher precipitation, by hump splitting, which results in a decrease in the pattern’s WL (*b*). Similar responses to local disturbances were found in the K and SG models. Figure 5 displays the results of two-dimensional 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 figure 6 shows.

### (b) Models exhibiting multiplicity of hybrid states

Numerical solutions of the LL and G models (§§2*a*,*d*) using PDE solvers point towards the existence of stable hybrid states in addition to periodic-pattern states.^{3} Figure 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, whereas the lowest branch corresponds to the bare-soil state. The branches in between correspond to stable hybrid states describing localized patterns, a few examples of which are shown in the right-hand panels. Solutions of this kind in one dimension and two dimensions have been found previously [20]. Figure 8 shows a partial bifurcation diagram for the G model in two dimensions. The upper line corresponds to a spot-pattern state,^{4} whereas the lower lines correspond to hybrid patterns with decreasing numbers of spots as the right-hand panels show. Note the difference between the hybrid solution branches in the two models: whereas, in the LL model, they all terminate at the same control parameter value *μ*_{f}, which coincides with the fold–bifurcation point of the periodic-pattern solution, in the G model the hybrid solution branches are slanted [32]—solutions with smaller numbers of spots terminate at lower *p* 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 behaviour 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 (*p* or *μ*). In the LL model, all localized pattern solutions terminate at the fold–bifurcation point *μ*_{f} (figure 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 figure 7 and the diagram obtained with the Swift–Hohenberg equation in figure 1. In the latter, there is a subrange (*λ*_{f}<*λ*<*λ*_{1}) 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 figure 8 was obtained from the periodic state (*a*) by an incremental decrease in *p*. Likewise, the hybrid states (*c*) and (*d*) were obtained from the states (*b*) and (*c*) by further incremental decreases of *p*. The degree of slanting increases as the root-to-shoot parameter *η* is increased.

## 4. 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} 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 behaviours. 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 figure 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 figure 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 the possible scenarios that 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 behaviour shown in figure 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, whereas, 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 WL 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, owing 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 figure 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 behaviour can be attributed to the positive pattern-forming infiltration and root-augmentation feedbacks. Both give advantage to plants at the rim of a patterned domain as compared with inner plants; the rim plants receive more run-off from the surrounding bare soil and experience weaker competition for soil water. These factors act against the retreat of vegetated domains. Processes that may favour such a retreat include soil erosion and roots exposure in sandy soils under conditions of high wind power [33], or insect outbreak [34]. 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 favour the formation of localized patterns alone. The results for the G model and its simplified version (SG) hint towards the possible role of the non-local water uptake by laterally extended root systems in inducing hybrid states. This non-local 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 non-local 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 [35], 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 [36].

## Funding statement

The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant no. (293825).

## Acknowledgements

We thank Arjen Doelman and Moshe Shachak for helpful discussions and Arik Yochelis for helping with the numerical continuation analysis.

## Footnotes

One contribution of 13 to a Theme Issue ‘Pattern formation in the geosciences’.

↵1 A pair of fronts propagating towards one another, however, can slow down and become stationary owing to repulsive front interactions [13].

↵2 The non-dimensional form of the equations was derived from the dimensional model, presented in [8], using the following transformation (the ∼ denotes the variables in [8]): , ,

*b*=*P*/*k*_{2},*w*=*W*/*k*_{1},*h*=*O*/*k*_{1}, , ,*μ*=d*t*_{0}, ,*f*=*W*_{0},*γ*=*k*_{2}/(*k*_{1}*c*),*ν*=*r*_{W}*t*_{0},*p*=*Rt*_{0}/*k*_{1}, where , .↵3 The application of numerical continuation methods for these models is harder than in the K, R and SG models and has not been pursued in this study.

↵4 The spot pattern is rhombic rather than hexagonal. Such patterns apparently exist in the G model, as in other pattern formation models [31]. In the present case, the system converged to the rhombic pattern following a disturbance of a hexagonal pattern. As it contains a direction of denser spots (vertical) along which the water stress is higher, transitions to hybrid states can be induced upon decreasing

*p*in the absence of local disturbances.↵5 We refer here to the Swift–Hohenberg equation as a prototype of pattern formation behaviours in systems showing bistability of uniform and patterned states, and to reference [12] for the implications of these general behaviours to the context of regime shifts. Similar pattern formation behaviours have been found with many other pattern formation models [11].

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.