Bistability and regular spatial patterns in arid ecosystems

Abstract

A variety of patterns observed in ecosystems can be explained by resource–concentration mechanisms. A resource–concentration mechanism occurs when organisms increase the lateral flow of a resource toward them, leading to a local concentration of this resource and to its depletion from areas farther away. In resource–concentration systems, it has been proposed that certain spatial patterns could indicate proximity to discontinuous transitions where an ecosystem abruptly shifts from one stable state to another. Here, we test this hypothesis using a model of vegetation dynamics in arid ecosystems. In this model, a resource–concentration mechanism drives a positive feedback between vegetation and soil water availability. We derived the conditions leading to bistability and pattern formation. Our analysis revealed that bistability and regular pattern formation are linked in our model. This means that, when regular vegetation patterns occur, they indicate that the system is along a discontinuous transition to desertification. Yet, in real systems, only observing regular vegetation patterns without identifying the pattern-driving mechanism might not be enough to conclude that an ecosystem is along a discontinuous transition because similar patterns can emerge from different ecological mechanisms.

Introduction

Regular pattern formation occurs in a variety of physical, chemical, and biological systems (Cross and Hohenberg 1993; Levin and Segel 1985). Most illustrative examples are the typical markings found on animal coats, such as spots on leopards or stripes on zebras (Aragon et al. 1998; Murray 1981; Murray 2001; Sekimura et al. 2000). In the context of morphogenesis, Turing (1952) studied reaction–diffusion equations describing a system composed of two substances: an activator, which promotes its own production (positive feedback), and an inhibitor, which slows down the activation process (negative feedback). He demonstrated that, when the activation effect occurs over a relatively short range and the inhibition effect occurs over a relatively long range (which occurs when the inhibitor diffuses much more rapidly than the activator), the system can be unstable to spatially heterogeneous perturbations. As a result, a perturbed system can evolve to regular spatial patterns that are stable in time (so-called stationary patterns). Such activator–inhibitor models, including short-range activation and long-range inhibition, can reproduce a large variety of regular patterns observed in nature (e.g., patterns on seashells: Meinhardt 1995; spatial organization of vegetation: Klausmeier 1999; Lejeune et al. 1999; Okayasu and Aizawa 2001; Rietkerk et al. 2002, 2004; von Hardenberg et al. 2001; spatial organization of mussel beds: van de Koppel et al. 2005).

A particular case of activator–inhibitor systems are the so-called activator-depleted substrate systems. In these systems, the inhibition effect results from the depletion of a substrate consumed during the production of the activator (Meinhardt 1995). In ecosystems, organisms (e.g., plants) can modify their environment in such a way that they increase the flow of a resource (e.g., water) toward them, leading to a local concentration of this resource. The resource is harvested from the surrounding areas, and the movement of the resource toward the organisms consequently leads to its depletion from farther away from the organisms. We refer to this particular case of activator-depleted substrate interactions as resource–concentration mechanism. Rietkerk et al. (2004) and Rietkerk and van de Koppel (2008) proposed that a variety of patterns in ecosystems can be explained by resource–concentration mechanisms.

Interestingly, activator-depleted substrate systems can exhibit bistability (Borckmans et al. 2004; Judd and Silber 2000; Murray 2001; Rietkerk et al. 2002; von Hardenberg et al. 2001), meaning that two stable states can be reached by the population, depending on the initial conditions and given the same parameter values. A bistable system might abruptly shift from one of its two stable states to the other because of gradual changes in external conditions. We refer to these abrupt shifts as discontinuous transitions. Bistability has often been observed in ecological models (e.g., Lotka 1925; Ludwig et al. 1978) and has interested many ecologists because it provides a theoretical explanation for “big effects from small causes” (Ricker 1963) as observed in empirical studies (e.g., Meijer 2000; Belyea and Malmer 2004; Mumby et al. 2007). We call bistability area the set of environmental factors and attributes of the system under which bistability occurs. It has been proposed that bistability and regular pattern formation are linked in activator-depleted substrate systems, meaning that specific shapes of the regular patterns could signal bistability (Rietkerk et al. 2004). This hypothesis is of particular interest for ecosystems where bistability can lead to sudden shifts between ecosystem states under gradual external changes, implying potential ecological and economic losses (Rietkerk et al. 2004; Scheffer et al. 2001). Identifying patterns that only occur in the bistability area could lead to indicators of imminent shifts from a healthy to a degraded state in ecosystems.

Previous studies successfully explained vegetation patterns in arid ecosystems as a result of Turing-like instability in reaction–diffusion systems. Klausmeier (1999) proposed a model of vegetation patterns in semiarid regions based on plants and water dynamics. In this model, the main biological mechanism driving the formation of the patterns is that vegetation increases water infiltration in the soil. Klausmeier (1999) showed that vegetation stripes form on gentle slopes. Sherratt (2005) and Sherratt and Lord (2007) performed further analyses of the same model.

Several authors have proposed expanded versions of Klausmeier's (1999) model, which also allowed for pattern formation on flat ground. Okayasu and Aizawa (2001), HilleRisLambers et al. (2001), and Rietkerk et al. (2002) proposed models where the water budget was decomposed into soil and surface water. In these models, vegetation pattern formation is mainly due to the differential infiltration rate of water in bare and vegetated soil. Later, the model of Rietkerk et al. (2002) was modified to address the effects of seed dispersal (Thomson et al. 2008; Pueyo et al. 2008) as well as the seasonality and stochasticity of rainfall (Guttal and Jayaprakash 2007; Ursino and Contarini 2006) on vegetation pattern formation.

von Hardenberg et al. (2001) and Meron et al. (2004) developed a model based on groundwater and vegetation dynamics only (no surface water). In their model, the positive feedback between vegetation and water is based on two mechanisms: vegetation patches compete for water due to water uptake by root and vegetation reduces evaporation of water locally. More recently, the same group proposed a related model with two different variables for soil and surface water (Gilad et al. 2004), and they expanded this model to study mixed woody–herbaceous vegetation (Gilad et al. 2007a) and the contribution of both feedbacks to pattern formation under different environmental conditions (Gilad et al. 2007b).

An alternative modeling approach was developed by Lejeune and colleagues (Lefever and Lejeune 1997; Lejeune and Tlidi 1999; Lejeune et al. 1999, 2002). In this approach, the mechanism underlying pattern formation originates from plant–plant interactions only, meaning that there is no explicit water dynamics. Their framework is based on short-range activation and long-range inhibition between plants through shading by the canopy and competing for water by lateral roots, respectively.

The spatially explicit models described above contributed to previously developed theories of arid ecosystems, which were mainly developed with models that ignored spatial interactions (so-called mean field models). Interestingly, these mean field models already predicted the possibility of bistability in arid ecosystems (Rietkerk and van de Koppel 1997; Rietkerk et al. 1997; van de Koppel et al. 1997). The rationale behind this result is that plants may be able to have access to an extra source of resources when their density exceeds a critical threshold (Scheffer et al. 2005). Examples include deep-rooting plants that can access water from deeper aquifers or plants that increase infiltration rates and thereby enable infiltration of water that otherwise would be lost from the system due to runoff (Holmgren and Scheffer 2001; Scheffer et al. 2005).

Very few studies compared the behavior of mean field and spatial models (Lejeune et al. 2002; Meron et al. 2004; Rietkerk et al. 2002; von Hardenberg et al. 2001). None of these previous studies performed an analysis with the explicit aim to study the link between mean field bistability, spatial bistability, and regular spatial vegetation patterns. Thus, whether resource–concentration provides a direct link between the occurrence of bistability and regular spatial pattern formation remains unclear.

In particular, under which conditions (i.e., the range of parameter values referring to environmental factors and attributes of the system) does resource–concentration lead to bistability of the system? Is the history of the system (especially previous spatial organization of the system) playing a role in the size of the bistability area? Are there shapes of the regular patterns that only occur in the bistability area?

The aim of this study is to test the hypothesis that regular spatial patterns and spatial bistability are linked in systems driven by a resource–concentration mechanism, using a modeling approach. Our ultimate goal is to find out whether regular spatial patterns can be used as indicators of bistability and, consequently, of possible discontinuous transitions in ecosystems with resource–concentration. We used a modified version of the model of Rietkerk et al. (2002), which describes the vegetation dynamics in arid ecosystems based on one main driving mechanism: water infiltrates faster into vegetated ground than into bare soil, leading to a positive feedback between soil water availability and vegetation locally and the depletion of water from farther away from the plants. We studied a mean field version and a spatially explicit version of the model where vegetation and water diffuse through space. Both the nonspatial and the spatial models present bistability.

We derived the conditions leading to bistability and discontinuous transitions in the nonspatial model. We then analytically calculated the conditions under which regular patterns are expected to occur, using the so-called Turing method (HilleRisLambers et al. 2001), and we compared these analytical results with spatially explicit numerical simulations. We deduced the bistability area of the spatial model and, consequently, the link between spatial bistability and regular vegetation patterns. Finally, we checked how the history of the system, meaning the spatial organization of the system in the past, affects vegetation pattern formation and the stability of these patterns under environmental changes, including the size of the bistability area.

The model

In the model of Rietkerk et al. (2002), three state variables are considered: plant density P (in grams of carbon per square meter), soil water W (in millimeters), and surface water O (in millimeters). The model assumes that rainfall events in arid and semiarid ecosystems occur at an intensity exceeding the infiltration capacity of the soil. Hence, part of the rainwater infiltrates into the soil, while the remainder produces surface water and runoff routed to other spatial locations.

The model of Rietkerk et al. (2002) is a good compromise between ecological realism and mathematical tractability. In its original form, however, all the surface water that can possibly infiltrate in the soil will eventually infiltrate because surface water cannot be lost from the system in this model. Here, the model deviates from previously developed mean field models (Rietkerk and van de Koppel 1997; Rietkerk et al. 1997) where precipitation that does not infiltrate into the soil is indeed lost from the system. This loss can induce bistability because, at low vegetation density, too much water can be lost and the system might not be able to sustain the plant population (Scheffer et al. 2005). Because surface water cannot be lost from the original model of Rietkerk et al. (2002), there is no positive feedback and no bistability in the nonspatial version of this model. In fact, bistability and pattern formation simultaneously emerge when spatial interactions are included in the spatially explicit version of this model. Thus, to adequately study the link between resource–concentration, bistability, and regular vegetation patterns, we need to decouple bistability and pattern formation in the model of Rietkerk et al. (2002). This means that we need an operational positive feedback, also in the nonspatial model. More specifically, we need to include a surface water loss term from the system whose magnitude decreases with plant density (Rietkerk and van de Koppel 1997; Scheffer et al. 2005). We thus modified the model of Rietkerk et al. (2002) by adding a loss term for surface water. Note that this is a realistic assumption for arid ecosystems where precipitation is often concentrated in relatively short and intense events (e.g., Schwinning and Sala 2004). Therefore, not all water may infiltrate into the soil, especially in areas deprived of vegetation (Grayson et al. 2006). This excess of water may be lost from the systems through runoff (Holmgren and Scheffer 2001; Grayson et al. 2006; Scheffer et al. 2005). This additional loss term strengthens the positive feedback between plant growth and water infiltration. In the following, we reiterate the model formulation of Rietkerk et al. (2002), including a description of how the surface water budget was revised.

In arid ecosystems, vegetation cover is often a two-phase mosaic composed of densely vegetated patches and bare soil areas (d'Herbès et al. 2001). The two phases of the mosaic mainly differ in their infiltration capacity for water (Galle et al. 1999; Tongway and Ludwig 2001). Vegetation improves the structure of the soil because it stimulates the biological activity in the soil, its root system forms channels and aerates the soil, and its canopy intercepts raindrops and prevents crust formation (e.g., Holzapfel and Mahall 1999; Kelly and Walker 1976; Pugnaire et al. 1996). Thus, infiltration is higher under vegetation than in bare soil (Rietkerk et al. 2000). In this model, the infiltration rate is assumed to asymptotically approach a maximum with increasing plant density (Walker et al. 1981), an assumption supported by field observations reported in Rietkerk et al. (2000). Thus, after a rain event, water runs off in bare areas and mainly infiltrates in vegetated patches, which act as sinks of water.

As already mentioned, we added a loss term l _o O for surface water to the model of Rietkerk et al. (2000). With this new term, in the absence of terrain variations, the dynamics of surface water depth are modeled as:

$$ \frac{{\partial O}}{{\partial t}} = R - \alpha O\frac{{P + {k_2}{W_{\text{o}}}}}{{P + {k_2}}} - {l_{\text{o}}}O + {D_{\text{o}}}\Delta O $$

(1)

where R is the rainfall (in millimeters per day), α is the maximum infiltration rate (per day), k ₂ is the saturation constant of water infiltration (in grams per square meter), W _o is a measure of the infiltration contrast between vegetated and bare soil (dimensionless), D _o is the diffusion coefficient for surface water (in square meters per day), ∆ is the Laplace operator in x and y, and l _o is the surface water loss rate which is meant to reflect losses due to runoff (per day). Note that we did not model runoff explicitly because this would require more hydrologically oriented modeling of catchment characteristics, which is beyond the scope of the current paper. Instead, we made the simplifying assumption that runoff water was immediately lost from the system and that the per capita runoff rate was constant (following Eppinga et al. 2009). It is also noteworthy that our main results regarding the relation between spatial patterns and spatial bistability are not affected when l _o = 0 (in the mean field model, however, bistability is lost when l _o = 0).

The infiltrated soil water is lost due to plant uptake, to evaporation and drainage, and to lateral subsurface flow because of capillary forces. The soil water dynamic is modeled as follows:

$$ \frac{{\partial W}}{{\partial t}} = \alpha O\frac{{P + {k_2}{W_{\text{o}}}}}{{P + {k_2}}} - g\frac{W}{{W + {k_1}}}P - {r_{\text{w}}}W + {D_{\text{w}}}\Delta W $$

(2)

where g is the maximum specific water uptake (in millimeters per gram per square meter per day), k ₁ is the half-saturation constant of specific growth and water uptake (in millimeters), r _w is the specific soil water loss due to evaporation and drainage (per day), and D _w is the diffusion coefficient for soil water (in square meters per day).

The dynamics of plant biomass are modeled as:

$$ \frac{{\partial P}}{{\partial t}} = cg\frac{W}{{W + {k_1}}}P - {\text{d}}P + {D_{\text{p}}}\Delta P $$

(3)

where c is the conversion of water uptake by plants to plant growth (in grams per millimeter per square meter), d is the specific loss of plant density due to mortality (per day), and D _p is the plant dispersal (in square meters per day).

For all parameters apart from l _o, we used the parameter value of Rietkerk et al. (2000). l _o was calibrated so that the shape of the patterns exhibited by the model correspond to the patterns observed by Barbier et al. (2006) for corresponding rainfall values in West Africa.

Analyses

Nonspatial model

The nonspatial model was analyzed for equilibria and their stability using standard methods for ordinary differential equations (for analytical details, see Appendix 1). There is a critical rainfall level at which the desert steady state changes stability. If rainfall is below this critical level, the desert state is stable; otherwise, the desert state is unstable. The critical rainfall level depends on plant physiology, soil characteristics, and water infiltration characteristics (Appendix 1).

Turing analysis

Linear stability analysis was used to determine whether regular patterns can form. The principle of this analysis is to investigate the fate of a small initial heterogeneity in an otherwise uniform system (HilleRisLambers et al. 2001; Murray 2001; Turing 1952).

We call the patterns predicted by this method Turing patterns and the range of parameter value in which they occur Turing instability range. The Turing instability range can be analytically calculated, and the details of the analytical derivation are presented in Appendix 2. The analytical calculation of the Turing stability conditions involves a quasi-steady-state assumption. More specifically, we assumed that surface water dynamics are much faster than plant and soil water dynamics (HilleRisLambers et al. 2001; Appendix 2).

Spatially explicit simulations

Two-dimensional numerical simulations were also performed to complement the Turing analysis predictions. They correspond to Euler integration of the finite-difference equations resulting from the discretization of the diffusion operator.

Different grid sizes and boundary conditions were used depending on the type of simulations performed. When simulations were run on grids with changing parameter values along the x- and y-directions (mimicking gradients), the spatial mesh consisted of a grid of 500 × 500 cells and a “no flux of matter boundary condition” at the border cells was assumed. The use of “periodic boundary conditions” instead was not intuitively appealing because some parameter values change along the two directions. Otherwise, the spatial mesh consisted of a grid of 100 × 100 cells with periodic boundary conditions. In this latter case, the shape and size of the vegetation patterns was not affected by the choice of the boundary conditions (similar patterns were observed with no flux of matter boundary conditions).

Patterns were followed along transitions of decreasing rainfall. For each rainfall value, simulations were run until stationary patterns were reached (i.e., steady state). The initial conditions of the simulations were either “random” or “sticky.” In case of random initial conditions, simulations were started by randomly introducing vegetation in 1% of the cells (for each of those cells: P = 50 g C m⁻², W and O: spatially homogeneous equilibrium). In case of sticky initial conditions, the outcome of the model for a different rainfall value was used as an initial condition. The initial condition used is mentioned in the legends of the figures. The spatial domain is such that one cell is 2 × 2 m², and the integration time step is daily.

Results

We start with the nonspatial model. As a consequence of the loss term for surface water (l _o ≠ 0) and the resulting stronger positive feedback between vegetation and soil water availability, bistability can occur in the nonspatial model (Appendix 1), which was not the case in the model of Rietkerk et al. (2002). Bistability means that two states, one corresponding to a high vegetation density and the other to a desert deprived of vegetation, can coexist for the same set of rainfall conditions. Then, the transition from one state to the other occurs through a discontinuous transition. When the loss term of surface water is present in the model (l _o ≠ 0), the main mechanism driving the bistability of the nonspatial model is W _o, a measure of the infiltration contrast between vegetated and bare areas (Fig. 1; Figs. 4 and 5 in Appendix 1). When the infiltration rate of water in bare soil is high, water can infiltrate everywhere in the system, independently of the presence of vegetation, and the feedback between vegetation and water is weak. On the contrary, when the infiltration rate of water in bare soil is low, the amount of water that can infiltrate depends crucially on the presence of vegetation, which increases the infiltration rate. This creates a large difference in soil water availability between bare and vegetated areas, and the positive feedback between vegetation and water is strong. For high values of W _o, meaning a weak positive feedback, the transition to desert is continuous and no bistability occurs. When W _o decreases, the strength of the positive feedback increases and bistability arises (Fig. 1). Thus, a stronger positive feedback between vegetation and soil water leads to a larger bistability area in the nonspatial model (Fig. 1). The other parameters affect the size of the bistability area, but not its existence (Fig. 4 in Appendix 1).

Fig. 1

Effect of the infiltration rate of water into the soil (W _o) on the bifurcation diagrams of the nonspatial and of the spatial models. Black/grey lines vegetation density at equilibrium (vertical axis) as a function of rainfall, R (horizontal axis) in the nonspatial model. Black stable equilibria of the nonspatial model, grey unstable equilibria of the nonspatial model (which mark the limit between the basin of attraction of the desert equilibrium and the other stable branch corresponding to vegetation [black line]). Bistability occurs in the range of R values where the unstable equilibrium exists. For readability, the desert equilibrium was not plotted. For W _o = 0.9, the transition is continuous (no bistability), whereas it is discontinuous for W _o = 0.2. Using the parameter settings of HilleRisLambers et al. (2001) and Rietkerk et al. (2002) and l _o = 0.06, the nonspatial model exhibits bistability, and thus discontinuous transitions, for W _o < 0.70. Green lines vegetation density (average vegetation biomass per cell of the lattice) at steady state as a function of rainfall R for W _o = 0.2 with sticky initial conditions. Simulation grids were 100 × 100 cells, and boundary conditions were periodic. To obtain the dark green line, we started at a value of R where patterns are expected to occur, we run the simulation from random initial conditions (1% of the sites, randomly chosen in the lattice initially contain a density of vegetation of P = 50; W and O are at their value of the homogeneous equilibrium). We then followed the branch up (resp. down) by increasing (resp. decreasing) the value of R of 0.1, and running each of the simulations for 3,000 time steps. We used the final state of the former simulation as an initial state for the following. The desert equilibrium is stable from R = 0 to R = 2.5. Patterns go extinct at R = 1 and recover at R = 2.5. The dashed light green line (not numerically obtained) represents the unstable equilibrium between the two stable equilibria in the bistability area. Note that, for W _o = 0.9, the transition is the same in the spatial and in the nonspatial model because no spatial pattern occur along the transition. c = 10, g = 0.05, k ₁ = 5, α = 0.2, d = 0.25, l _o = 0.06, k ₂ = 5, r _w = 0.2

In the spatial model, the homogeneous desert is stable under the same range of conditions as in the nonspatial model, but vegetation can survive under harsher environmental conditions than predicted by the nonspatial model (Fig. 1). This is associated with pattern formation of the vegetation, which we explore in Fig. 2.

Fig. 2

Turing and non-Turing vegetation patterns, spatial and nonspatial bistability. a State diagram of the nonspatial model. Black homogeneous vegetation, grey bistability, white desert. Homogeneous vegetation survives until the hysteresis curve: $d = \left( {Rcg\alpha {W_{\text{o}}}} \right)/\left( {R\alpha {W_{\text{o}}} + {r_{\text{w}}}{k_1}\left( {\alpha {W_{\text{o}}} + {l_{\text{o}}}} \right)} \right)$ (Appendix 1). On the left of the white dashed vertical line (W _o = 0.70), transitions are discontinuous in the nonspatial model when bifurcation parameter, R, is changed, whereas they are continuous on its right. b Analytical solution for the Turing instability range (Appendix 2). Black homogeneous vegetation cover, light grey Turing patterns, white homogeneous bare soil. Dashed white line same interpretation as in a. c Spatially explicit simulations. Grid size = 500 × 500. Simulation time = 5,000 time steps. d Bistability area of the spatial model obtained by combining a and c: black homogeneous vegetation and no bistability (case 1), dark green homogeneous vegetation coexisting with homogeneous desert (case 2), light green Turing and non-Turing patterns coexisting with homogeneous desert (case 4). The bistability area of the spatial model is the total green area (dark and light green). Red Turing patterns outside the bistability area (case 3). e Areas of Turing and non-Turing patterns obtained by combining b and c. Black homogeneous vegetation, light grey Turing patterns predicted by the linear analysis, dark grey non-Turing patterns predicted by the spatially explicit simulation. The spatially explicit patterns from c can be seen in transparency behind the light and dark grey areas. The white dashed vertical line is the limit between discontinuous and continuous transitions in the nonspatial model (same as in a and b), and the gray dashed line is the same in the spatial model. c = 10, g = 0.05, k ₁ = 5, α = 0.2, d = 0.25, l _o = 0.06, k ₂ = 5, r _w = 0.2. Note that the different zones on d and e are obtained by comparison of the figures, not by direct numerical analysis

We now start the pattern analyses with a Turing analysis, which provides information about when regular patterns are expected to form, by investigating the fate of a small initial heterogeneity in a uniform vegetation cover (Levin and Segel 1985; Murray 2001; Turing 1952). Four possible model outcomes are possible: (1) homogeneous vegetation cover and no bistability; (2) homogeneous vegetation cover and bistability; (3) vegetation patterns and no bistability; (4) vegetation patterns and bistability. The Turing patterns, predicted by this analysis, occur always in a small zone under the harshest environmental conditions where the nonspatial model predicts vegetation maintenance (Fig. 2b; for the analytical calculations of the Turing instability range, see Appendix 2). Turing patterns occur both in conditions where vegetation coexists with homogeneous desert (Fig. 2b; case 4; on the left of the dashed white line) and in a range of conditions where Turing patterns are the only stable state of the system (Fig. 2b; case 3; on the right of the dashed white line). We now need to identify the bistability area of the spatial model, which we do by using spatially explicit simulations.

Spatially explicit simulations of the model show that vegetation patterns, more precisely labyrinths and spots, arise outside the Turing instability range (Fig. 2c). These non-Turing patterns arise due to the nonlinearities of the model, which are not taken into consideration when analytically calculating the Turing instability range, where the nonspatial model is linearized in the vicinity of the homogeneous vegetation equilibrium (Appendix 2). Non-Turing patterns prolong the Turing patterns in the desert area where the nonspatial model did not predict a possible maintenance of the vegetation (Fig. 2e). In the spatial model, vegetation can survive in the form of non-Turing patterns under harsher conditions than predicted by the nonspatial model. These non-Turing patterns are in a parameter zone where the homogeneous desert equilibrium is also stable (compare Fig. 2a, e). Numerical simulations confirm that the homogeneous desert state is stable under the same range of conditions in the spatial and in the nonspatial models (e.g., Fig. 1). Therefore, non-Turing patterns belong to the bistability area of the spatial model (case 4). The spatial bistability area is thus composed of the bistability area of the nonspatial model and, additionally, the area of the non-Turing patterns (Fig. 2d). In this case, the spatial bistability area is, therefore, always larger than the nonspatial bistability area due to the emergence of the non-Turing patterns in the spatial model. Note that spot-like and labyrinthine non-Turing patterns similar to those obtained here were recently observed and studied in Muratov (1997), Hastings et al. (1997), and Morozov and Petrovskii (2009).

As mentioned earlier, in our model, non-Turing patterns always coexist with stable homogeneous desert, but this is not necessarily the case for the Turing patterns (case 3 or 4; see Fig. 2b). We traced bifurcation graphs with rainfall, R, as a bifurcation parameter for different values of W _o, and we found that Turing patterns always prolong into non-Turing patterns under decreasing rainfall. This happens at the right of the white dashed line as well, although the range of conditions under which non-Turing patterns occur is then very small (Fig. 2e).

Our study has two important implications. The first implication is that, in this model, the occurrence of vegetation patterns indicates that the system is along a discontinuous transition. In particular, labyrinths and spots occur in the bistability area (case 4), meaning that they occur when the system is degrading, but once the system collapsed to the desert state, they would not occur if the system is recovering. So, the succession of patterns (in time) indicates whether the system is degrading or restoring, and the shape of the patterns indicates how far the system is from discontinuous transition to extinction. When rainfall decreases, the vegetation cover changes from homogeneous to gaps, to labyrinths, and to spots before becoming extinct (Fig. 2c). The spot patterns are thus the last patterns to occur along discontinuous transition to desertification. Therefore, they indicate proximity to desertification in this resource–concentration system, as already suggested by other model studies (Meron et al. 2004; Rietkerk et al. 2002, 2004; von Hardenberg et al. 2001). This result is not affected by the presence of the loss term of surface water, l _o O.

The second implication is that the parameter range under which discontinuous transitions occur in the spatial model is larger than in the nonspatial model (Fig. 2e). In particular, bistability occurs in the spatial model under higher values of W _o, meaning for weaker positive feedbacks between vegetation and water, compared with the nonspatial model. Moreover, in the spatial model, a stronger positive feedback leads to a larger bistability area (Fig. 2d).

Simulations show that there is another phenomenon that affects the area of regular vegetation pattern occurrence and, consequently, the area of spatial bistability: the type of initial condition (random or sticky; see the “Analyses” section for details), i.e., the history of the patterns (Fig. 3). Along a transition (e.g., R decreases), the usual sequence of patterns obtained from a random initial condition is homogeneous vegetation cover, gaps, labyrinths, spots, and desert when rainfall decreases (Fig. 3a–e; as commonly observed in other models of vegetation dynamics in arid ecosystems; Lejeune et al. 1999; Okayasu and Aizawa 2001; Rietkerk et al. 2002; von Hardenberg et al. 2001). However, with sticky initial conditions (i.e., using patterns obtained under other rainfall conditions as initial state), the patterns remain under harsher conditions (“stickiness of the patterns”; compare Fig. 3e, j, k). This means that the spatial organization of the vegetation makes the system more resistant to harshness of environmental conditions. Moreover, different pattern shapes emerge (new topology) with sticky initial conditions (Fig. 3f–j). Even band-like patterns are found on bare ground (Fig 3i; as also observed by Lejeune et al. 1999).

Fig. 3

Patterns occurring before extinction in discontinuous transitions. Spatially explicit simulations. Lattices of size = 100 × 100. Simulations time = 2,000 time steps. First row (a–e) patterns from random initial condition. From left to right, R = 1.8, R = 1.5, R = 1.3, R = 1.2, R = 1.00. Second row (f–j) patterns from sticky initial conditions. From left to right, R = 1.8, R = 1.5, R = 1.35, R = 1.25, R = 1.00. W _o = 0.2; k ₂ = 5; α = 0.2; r _w = 0.2; D _p = 0.1; D _o = 100; D _w = 0.1; d = 0.25; k ₁ = 5; g = 0.05; c = 10; l _o = 0.06. k Zoomed view of the low-rainfall extremity of the bifurcation graph, plotted with sticky initial conditions (as in Fig. 1) in solid line and with random initial conditions in dashed line. Note that the succession of patterns (a–e) can be seen as a cross-section of Fig. 2e for W _o = 0.2. Thus, a shows Turing patterns, while the rest (b–e) show non-Turing patterns

Discussion

We analyzed a model of vegetation dynamics in arid ecosystems, including a positive feedback between vegetation density and surface water infiltration, which leads to an increase in soil water availability locally and its consequent depletion from the plants' surroundings. The strength of this feedback depends on the infiltration contrast between vegetated and bare areas. Adding a loss term for surface water in the model of Rietkerk et al. (2002), we show that bistability can arise in the nonspatial and in the spatial models, leading to potential discontinuous transitions between vegetation and desert states. We provide the first complete analysis of the link between bistability and regular vegetation patterns in a model driven by resource concentration. Note that, in this model, it is very likely that different shapes of patterns coexist under certain ranges of conditions (Sherratt and Lord 2007). However, here, we are only interested in the bistability that occurs between a homogeneous desert state and a vegetated state, which is of relevance for understanding the behavior of systems approaching desertification.

In both the nonspatial and the spatial models, bistability arises if the positive feedback is strong enough. Bistability arises because a threshold vegetation biomass is necessary to activate the feedback between vegetation and water availability: if vegetation density is too low, vegetation is not able to concentrate enough water for its own survival. We show that the bistability area is always larger in the spatial model than in the nonspatial model, mainly because of the non-Turing patterns. Vegetation can live under harsher environmental conditions in the form of spatial patterns, but under these harsh conditions, the system is also bistable and thus likely to suddenly collapse to another stable state if perturbed (in agreement with Lejeune et al. 2002; Meron et al. 2004; Rietkerk et al. 2002; von Hardenberg et al. 2001).

Our model builds upon a considerable body of arid ecosystem models suggesting the possibility of regular pattern formation and discontinuous transitions (e.g., Gilad et al. 2004; HilleRisLambers et al. 2001; Klausmeier 1999; Lejeune and Tlidi 1999; Rietkerk and van de Koppel 1997; Rietkerk et al. 1997; Sherratt and Lord 2007; von Hardenberg et al. 2001). These previous models showed that bistability can occur in two different ways: plants may be able to harvest a resource that would otherwise be lost from the system (Rietkerk et al. 1997; Scheffer et al. 2005) or plants may be able to increase their access to resources by harvesting resources from their surrounding area (Rietkerk et al. 2002, 2004). The first type of bistability can be predicted by mean field models (Rietkerk and van de Koppel 1997; Rietkerk et al. 1997). The second type is bistability between non-Turing patterns and a desert state. For the first time, we examined the relation between these two types of bistability. More specifically, we tested (1) whether the Turing instability is always followed by the lower limit of the mean field bistability area when conditions get harsher and (2) whether the Turing patterns are always followed by the existence of non-Turing patterns when conditions get harsher. Both conditions are necessary requirements to verify the hypothesis that regular vegetation patterns indicate proximity to discontinuous transitions (Rietkerk et al. 2004). Our results revealed that both conditions are fulfilled in our model. In other words, if regular vegetation patterns are observed in our system, they indicate that the system is along a discontinuous transition to desertification and that the system would go extinct in a discontinuous manner if the external conditions would deteriorate. In particular, spots are the last patterns to occur before an extinction in a discontinuous manner, and thus they indicate proximity to a discontinuous shift in our model system with resource–concentration, mathematically formalizing the hypothesis of Rietkerk et al. (2004). Interestingly, however, our study also showed that even homogeneously covered vegetation may exhibit mean field bistability. Hence, the absence of regular vegetation patterns does not rule out the possibility of discontinuous transitions.

It seems quite well established that spots or irregular spot-like patterns are indicators of imminent desertification of arid ecosystems, in continuous and discontinuous ways (e.g., Kéfi et al. 2007, 2008; Meron et al. 2004; Rietkerk et al. 2002, 2004; von Hardenberg et al. 2001). However, their role as indicators specific to discontinuous transitions depends on the ecological mechanism generating them. Different ecological mechanisms can lead to apparently similar patterns. Gilad et al. (2004) developed a model of arid ecosystem dynamics quite similar to ours, but in which they included a positive feedback mechanism between aboveground and belowground biomass (so-called root-augmentation feedback) besides the infiltration feedback. In this model, they showed that the general sequence of biomass patterns (i.e., gaps, labyrinths, spots) is not affected, but their model allows for monostability of spot patterns: there is a precipitation range where the bare soil state is no longer stable but spots are still stable. Observing spots in nature is thus not enough to infer whether the system is close to a continuous or discontinuous transition to desertification. We need to know the underlying ecological mechanisms to be able to conclude about the value of spots as indicators of proximity from discontinuous shifts. If the resource–concentration mechanism operates, for example, soil water availability should be higher under vegetation patches than in bare areas.

It is ecologically relevant to know not only when a discontinuous transition will occur, but also how reversible such a transition may be. Restoring vegetation may be especially difficult under environmental conditions where the bare state is stable, i.e., when the infiltration rate of water in bare soil (W _o) is low in our model. Note that the degree of reversibility of a discontinuous transition is not indicated by the shape of the vegetation patterns since the same sequence of gaps–labyrinths–spots occurs for both low and higher values of W _o. However, in our model, the degree of reversibility is determined by the amount of water that is lost from the system when vegetation is gone (Rietkerk and van de Koppel 1997; Rietkerk et al. 1997; van de Koppel et al. 1997; Scheffer et al. 2005) and by the amount of water that can be harvested from surrounding areas (e.g., Rietkerk et al. 2002; Gilad et al. 2004, 2007a, b). Hence, a suitable indicator for the degree of reversibility of a proximate discontinuous transition should reflect both these characteristics. In line with previous studies by Shachak et al. (1998, 2008) and Gilad et al. (2004, 2007a), a possible indicator of reversibility could be an “infiltration ratio,” meaning the ratio between the infiltration rate in a vegetation patch and the calculated infiltration rate if the patch was bare (but with soil and surface water still being in the same equilibrium; see Fig. 6 in Appendix 3). This infiltration ratio could be relatively easily measured in the field.

Another interesting result is that the patterns are sticky, meaning that patterns are relatively similar to the patterns they originate from. Because of this stickiness, starting from a spatially organized system and gradually decreasing rainfall, vegetation can persist long into conditions where simulations starting from random initial conditions do not lead to vegetation maintenance. However, this only happens if the changes in rainfall are very small and slow; otherwise, the system collapses to a desert.

Both bistability and stickiness imply that the history of the system can affect the shape of the vegetation patterns observed under a given external condition. Although the patterns are static, looking at their shape at a snapshot in time could thus already provide information about the history of the system. Stickiness could also play an important role in the restoration of degraded systems where planting seeds or seedlings in spatial patterns might increase the recolonization success.

In a recent review about ecosystems exhibiting regular vegetation patterns, Rietkerk and van de Koppel (2008) highlighted the importance of scale-dependent feedback mechanisms between organisms and their environment in a variety of ecosystems, such as wetlands, savannas, mussel beds, coral reefs, ribbon forests, intertidal mudflats, and marsh tussocks. Besides showing regular patterns, these systems may exhibit bistability if scale-dependent feedback is indeed a main driver of their dynamics. Evidence of alternative stable states and hysteresis has actually already been reported for tidal flats (van de Koppel et al. 2001) and coral reefs (Mumby et al. 2007). Using our knowledge about the ecological mechanisms operating in these systems and, in particular, taking advantage of the scale-dependent feedback by using ecosystem engineers, could increase restoration success (Suding et al. 2004). Moreover, it might be possible to use the vegetation patterns themselves, by mimicking the patterns observed in healthy systems, to restore degraded systems.

Finding indicators specific to discontinuous transitions has attracted lots of attention in ecology lately, the main reason being the unexpected aspect of these transitions and their potential irreversibility. It is noteworthy though that transitions that are mathematically discontinuous are not always relevant for field ecology. In spatially organized systems, with consumer–resource positive feedback, appearance–disappearance of individuals may occur in patches of individuals, which are the spatial units of the system. This would lead to a mathematically discontinuous transition, but would not be perceived as such in the field. Moreover, even if the transition is discontinuous, a very small hysteresis loop is still very easy to reverse. In our model, we show that, for weak positive feedback, transitions to desert are discontinuous but the area of bistability, and therefore, the size of the hysteresis loop, is very small. Besides finding indicators of imminent shifts, a challenge for systems exhibiting catastrophic shifts is to be able to estimate the size of the hysteresis loops.

Appendices

Appendix 1: spatially homogeneous solutions

There is one trivial spatially homogeneous solution (i.e., lacking spatial structure): $\left( {{P^*} = 0,{W^*} = \frac{{\alpha {W_{\text{o}}}R}}{{\left( {\alpha {W_{\text{o}}} + {l_{\text{o}}}} \right){r_{\text{w}}}}},{O^*} = \frac{R}{{\alpha {W_{\text{o}}} + {l_{\text{o}}}}}} \right),\quad {\text{which}}\,{\text{is}}\,{\text{stable}}\,{\text{if:}}$

$$ d > \frac{{cg\alpha {W_{\text{o}}}R}}{{\alpha {W_{\text{o}}}R + {k_1}{r_{\text{w}}}\left( {\alpha {W_{\text{o}}} + {l_{\text{o}}}} \right)}}. $$

(4)

In the following, we call this equilibrium desert equilibrium. The stability condition of the desert equilibrium can be rewritten as follows:

$$ R\kern1.5pt<\kern1.5pt\frac{{{r_{\text{w}}}d{k_1}\left( {{l_{\text{o}}} + \alpha {W_{\text{o}}}} \right)}}{{\alpha {W_{\text{o}}}\left( {c{g} - d} \right)}}. $$

(5)

At this point, there is a transition towards a homogeneous vegetated state. This transition can be continuous or discontinuous. When the water infiltration feedback is strong enough, it is discontinuous, meaning that, for a range of rainfall rates, the model exhibits bistability between a uniform vegetated state and a bare state. The transition is discontinuous if:

$$ {W_{\text{o}}}\kern1.5pt<\kern1.5pt\frac{{{l_{\text{o}}}\left( {d{k_2} - c\left( {{g}{k_2} + {k_1}{r_{\text{w}}}} \right)} \right) + \sqrt {\vartheta } }}{{2\alpha {k_2}\left( {c{g} - d} \right)}} $$

(6)

in which:

$$ \vartheta = {l_{\text{o}}} \times \left( {4\alpha c} \right.\left( {c{g} - d} \right){k_1}{k_2}{r_{\text{w}}} + \left( {{k_2}\left( {c{g} - d} \right){{\left. { + c{k_1}{r_{\text{w}}}} \right)}^2} \times {l_{\text{o}}}} \right). $$

Nontrivial spatially homogeneous equilibria are solutions of the system:

$$ \left\{ {O = \frac{{\begin{array}{*{20}{c}} {W = \frac{{d{k_1}}}{{cg - d}}} \\ {R\left( {P + {k_2}} \right)} \\ \end{array} }}{{\begin{array}{*{20}{c}} {\alpha \left( {P + {k_2}{W_{\text{o}}}} \right) + {l_{\text{o}}}\left( {P + {k_2}} \right)} \\ {A{P^2} + BP + C = 0\left( {\text{Pol}} \right)} \\ \end{array} }}} \right. $$

with:

$$ A = \frac{d}{c}\left( {\alpha + {l_{\text{o}}}} \right), $$

$$ B = \frac{d}{c}\left( {\alpha {k_2}{W_{\text{o}}} + {l_{\text{o}}}{k_2}} \right) + \left( {\alpha + {l_{\text{o}}}} \right){r_{\text{w}}}\frac{{d{k_1}}}{{cg - d}} - R\alpha, $$

$$ C = {r_{\text{w}}}\frac{{d{k_1}}}{{cg - d}}\left( {\alpha {k_2}{W_{\text{o}}} + {l_{\text{o}}}{k_2}} \right) - R\alpha {k_2}{W_{\text{o}}}. $$

The polynomial Pol is convex. It has two solutions (one stable equilibrium and one unstable equilibrium) as soon as its min value $\left( {{P_{\min }} = - B/\left( {2A} \right)} \right)$ is negative. See Fig. 4 for the effect of l _o, d, and k ₁ on the bifurcation graph of the nonspatial model.

Fig. 4

Effect of the parameters (other than W _o) on the bifurcation diagrams of the nonspatial model. Black stable equilibria, grey unstable equilibria. c = 10, g = 0.05, W _o = 0.2, k ₁ = 5, α = 0.2, d = 0.25, l _o = 0.06, k ₂ = 5, r _w = 0.2

Fig. 5

State diagram of the nonspatial model: roles of l _o and W _o. Black homogeneous vegetation, grey bistability, white desert. Vertical axis l _o varies between 0 and 0.5. Horizontal axis W _o varies between 0 and 1. Other parameter values: c = 10, g = 0.05, k ₁ = 5, α = 0.2, d = 0.25, R = 2, k ₂ = 5, r _w = 0.2. For l _o = 0, the system shows no bistability independently of the value of W _o. As soon as l _o ≠ 0, bistability may occur in the mean field model under certain ranges of conditions (grey area on the figure), and W _o is the most important parameter determining this range of conditions (there is a critical value of W _o above which there is no bistability anymore)

Appendix 2: stability to spatially heterogeneous perturbations

We linearize the system of equation around the nontrivial equilibrium (P*, W*, O*):

$$ P = {P^*} + \varepsilon \left( {x,y,t} \right), $$

(7)

$$ W = {W^*} + \phi \left( {x,y,t} \right), $$

(8)

$$ O = {O^*} + \theta \left( {x,y,t} \right), $$

(9)

$$ \frac{{\partial P}}{{\partial t}} = \frac{{\partial \varepsilon }}{{\partial t}} = f\left( {P,W,O} \right) = \frac{{\partial f}}{{\partial P}}\left( {P - {P^*}} \right) + \frac{{\partial f}}{{\partial W}}\left( {W - {W^*}} \right) + \frac{{\partial f}}{{\partial O}}\left( {O - {O^*}} \right) $$

(10)

(the same for the two other equations).

We choose (HilleRisLambers et al. 2001):

$$ \varepsilon \left( {x,y,t} \right) = \varepsilon (t)\cos \left( {\frac{x}{{{l_1}}}} \right)\cos \left( {\frac{y}{{{l_2}}}} \right), $$

(11)

$$ \phi \left( {x,y,t} \right) = \phi (t)\cos \left( {\frac{x}{{{l_1}}}} \right)\cos \left( {\frac{y}{{{l_2}}}} \right), $$

(12)

$$ \varphi \left( {x,y,t} \right) = \varphi (t)\cos \left( {\frac{x}{{{l_1}}}} \right)\cos \left( {\frac{y}{{{l_2}}}} \right). $$

(13)

By replacing the partial derivatives (∂f/∂P, ∂f/∂W, ∂f/∂O) with their expressions, we obtain:

$$ \frac{{\partial \varepsilon }}{{\partial t}} = {P^*}\frac{{{{\left( {cg - {R_{\text{esp}}}} \right)}^2}}}{{cg{k_1}}}\phi - \varepsilon {Q^2}{D_{\text{p}}}, $$

(14)

$$ \frac{{\partial \varphi }}{{\partial t}} = \left( {\frac{{R{k_2}\left( {1 - {W_{\text{o}}}} \right)\alpha }}{{\left( {\alpha \left( {{P^*} + {k_2}{W_{\text{o}}}} \right) + {l_{\text{o}}}\left( {{P^*} + {k_2}} \right)} \right)\left( {{P^*} + {k_2}} \right)}} - \frac{d}{c}} \right)\varepsilon - \left( {{P^*}\frac{{{{\left( {cg - d} \right)}^2}}}{{{k_1}{c^2}g}} + {r_{\text{w}}}} \right)\phi + \alpha \frac{{{P^*} + {k_2}{W_{\text{o}}}}}{{{P^*} + {k_2}}}\theta - \phi {Q^2}{D_{\text{w}}}, $$

(15)

$$ \frac{{\partial \theta }}{{\partial t}} = - \frac{{\alpha R{k_2}\left( {1 - {W_{\text{o}}}} \right)}}{{\left( {\alpha \left( {{P^*} + {k_2}{W_{\text{o}}}} \right) + {l_{\text{o}}}\left( {{P^*} + {k_2}} \right)} \right)\left( {{P^*} + {k_2}} \right)}}\varepsilon - \left( {\alpha \frac{{{P^*} + {k_2}{W_{\text{o}}}}}{{{P^*} + {k_2}}} + {l_{\text{o}}}} \right)\theta - \theta {Q^2}{D_{\text{o}}} $$

(16)

with ${Q^2} = \left( {\frac{1}{{l_1^2}} + \frac{1}{{l_2^2}}} \right)$.

We use a quasi-steady-state approach, i.e., we assume that the dynamics of surface water acts on a much faster scale than changes in W and P (HilleRisLambers et al. 2001). θ is expressed as a function of ε and φ by solving ∂θ/∂t = 0. We replace θ by its expression in the two other equations. J, the Jacobian matrix at equilibrium, of the resulting system of two equations is given by:

$$ J = \left( {\begin{array}{*{20}{c}} { - {D_{\text{p}}}{Q^2}} & {{P^*}\frac{{{{\left( {cg - d} \right)}^2}}}{{cg{k_1}}}} \\ {E\left( {1 - \frac{{\alpha \left( {{P^*} + {k_2}W_{\text{o}}} \right)}}{{\left( {{D_{\text{o}}}{Q^2} + l} \right)\left( {{P^*} + {k_2}} \right) + \alpha \left( {{P^*} + {k_2}{W_{\text{o}}}} \right)}}} \right) - \frac{d}{c}} & { - {P^*}\frac{{{{\left( {cg - d} \right)}^2}}}{{{k_1}{c^2}g}} - {r_{\text{w}}} - {D_{\text{w}}}{Q^2}} \\ \end{array} } \right) $$

(17)

with $E = \frac{{R{k_2}\left( {1 - {W_{\text{o}}}} \right)\alpha }}{{\left( {{P^*} + {k_2}} \right)\left( {\alpha \left( {{P^*} + {k_2}{W_{\text{o}}}} \right) + {l_{\text{o}}}\left( {{P^*} + {k_2}} \right)} \right)}}$ whose Trace(J) < 0.

Pattern formation occurs when Det(J) < 0. Det(J) is given by:

$$ \begin{array}{*{20}{c}} {{\text{Det}}(J) = F\left( {{Q^2}} \right) = {D_{\text{p}}}{Q^2}\left( {\frac{{{P^*}{{\left( {cg - d} \right)}^2}}}{{{k_1}{c^2}g}} + {r_{\text{w}}} + {D_{\text{w}}}{Q^2}} \right)} \hfill \\ {\quad \quad \quad \quad \quad \quad \quad + \frac{{{P^*}{{\left( {cg - d} \right)}^2}}}{{{k_1}cg}}\left( {\frac{d}{c} + E\frac{{\alpha \left( {{P^*} + {k_2}{W_{\text{o}}}} \right)}}{{\left( {{P^*} + {k_2}} \right)\left( {{Q^2}{D_{\text{o}}} + {l_{\text{o}}}} \right) + \alpha \left( {{P^*} + {k_2}{W_{\text{o}}}} \right)}}} \right)} \hfill \\ {\quad \quad \quad \quad \quad \quad \quad - \frac{{{P^*}{{\left( {cg - d} \right)}^2}}}{{{k_1}cg}}E} \hfill \\ \end{array} . $$

(18)

Appendix 3: the infiltration ratio as an indicator of reversibility of discontinuous shifts

Fig. 6

Illustration of the “infiltration ratio” as a potential indicator of the reversibility of a discontinuous shift toward desertification. The graph shows the ratio between the actual infiltration rate and the infiltration rate if the patch was bare soil (but with the same water conditions). The result reflects that, at low infiltration rates of bare soil, it will be more difficult to restore vegetation