Gradual Regime Shifts in Fairy Circles

TL;DR

This study combines empirical data and mathematical modeling to explore how fairy circles in Namibia undergo gradual regime shifts driven by water-vegetation interactions, illustrating self-organization and stable state transitions in pattern-forming ecosystems.

Contribution

It provides new evidence supporting fairy circles as a self-organized pattern caused by water-vegetation dynamics and highlights the gradual nature of regime shifts in such ecosystems.

Findings

Fairy circles are driven by water-vegetation interactions.

Transitions between stable states are spatially confined.

Regime shifts occur gradually through cascades of local transitions.

Abstract

Large responses of ecosystems to small changes in the conditions--regime shifts--are of great interest and importance. In spatially extended ecosystems, these shifts may be local or global. Using empirical data and mathematical modeling, we investigated the dynamics of the Namibian fairy circle ecosystem as a case study of regime shifts in a pattern-forming ecosystem. Our results provide new support, based on the dynamics of the ecosystem, for the view of fairy circles as a self-organization phenomenon driven by water-vegetation interactions. The study further suggests that fairy circle birth and death processes correspond to spatially confined transitions between alternative stable states. Cascades of such transitions, possible in various pattern-forming systems, result in gradual rather than abrupt regime shifts.

Tables

Table 1. Table 1: Relative change in the existence range of the hybrid states (snaking range), due to varied parameters compared to those used in the main text. In each case, one parameter is changed by either + 10 % percent 10 +10\% or − 10 % percent 10 -10\% , and the relative existence range is shown.

parameter	$+10\%$	$-10\%$
$\lambda$	$0.8913$	$1.1269$
$\eta$	$1.0179$	$0.9664$
$\nu$	$1.0604$	$0.9404$
$\rho$	$1.0183$	$0.9819$
${\delta }_w$	$1.0556$	$0.9413$

Equations

$${\partial }_{T}B$$

$${\partial }_{T}W$$

$${\psi }_i\left(S\right)=\sum _{j=0}^{S-1}{P}_{(i-j-L)}.$$

$$C_{\alpha }\left(S\right)=\frac{\sum _{i=1}^N({\psi }_i\left(S\right)-\overline{\psi \left(S\right)})\cdot ({\alpha }_i-\bar{\alpha })}{\sqrt{\sum _{i=1}^N({\psi }_i\left(S\right)-\overline{\psi \left(S\right)}{)}^2}\cdot \sqrt{\sum _{i=1}^N({\alpha }_i-\bar{\alpha }{)}^2}}.$$

$$B_T$$

$$W_T$$

$$H_T$$

$$G_B(\text{X},T)$$

$$G_W(\text{X},T)$$

$$G(\text{X},{\text{X}}^{\prime },T)$$

$$I$$

$$G_B(\text{X},T)$$

$$G_W(\text{X},T)$$

$$B_T$$

$$W_T$$

$$b=\frac{B}{K};w=\frac{W\Lambda }{K\Gamma };t=MT;x=X\sqrt{M/D_B}.$$

$${\partial }_{t}b=\lambda wb(1+\eta b{)}^2(1-b)-b+{\nabla }^{2}b,$$

$${\partial }_{t}w=p-\nu w(1-\rho b)-\lambda wb(1+\eta b{)}^2+{\delta }_w{\nabla }^{2}w.$$

$$\lambda =\frac{K\Gamma }{M};\eta =EK;p=\frac{\Lambda P}{K\Gamma M};\nu =\frac{N}{M};\rho =R;{\delta }_w=\frac{D_W}{D_B}.$$

Full text

The response of ecosystems to climate variability and anthropogenic disturbances is a fundamental aspect of ecology. Much attention has been devoted recently to large responses of ecosystems to small environmental changes or disturbances. Such responses, often termed “catastrophic regime shifts,” are conceived of as abrupt transitions between two alternative stable states that occur uniformly across the ecosystemScheffer2001nature . Spatially extended systems can respond in different ways to varying conditions, including pattern formationrietkerk2004science ; Deblauwe2008geb ; Meron2012ecomod and spatially confined transitions to alternative stable statesBel2012theo_ecol . When one of the alternative stable states is spatially patterned, a multitude of additional stable states can appear, each a mosaic of fixed domains that alternate between the uniform and the patterned states Lejeune2002pre ; Knobloch2008nonlinearity ; Meron2012ecomod . The existence of such hybrid states can strongly affect the dynamics of state transitions in fluctuating environmentsGandhi2015prl , and may lead to gradual rather than abrupt shiftsBel2012theo_ecol ; zelnik2013regime . This prediction of pattern formation theory has never been tested either in a real ecosystem or in a model describing the dynamics of a specific ecosystem.

Pattern formation is widespread in natural ecosystemsrietkerk2008tree ; Deblauwe2008geb , but good case models for studying gradual regime shifts are not abundant. An outstanding candidate is the Namibian fairy circle (NFC) ecosystem, which is fairly homogeneous and undisturbed, and recently has been the subject of intense research tlidi2008lnp ; picker2012ants ; tschinkel2012life ; juergens2013biological ; cramer2013namibian ; Fernandez-Oto2014philos_trans_A ; getzin2015adopting . The NFC ecosystem consists of a uniform matrix of perennial grass, punctured by circular gaps of sandy bare soil–the fairy circles–that on landscape scales, form nearly periodic patternspicker2012ants ; getzin2015adopting .

Various explanations for the formation of fairy circles have been suggested, including the release of poisonous gas and the feeding habits of ants and termitesbecker2000fairy ; juergens2013biological ; getzin2015adopting . These explanations, however, have not uncovered, as of yet, the small-scale feedbacks needed to account for the emergence of the large-scale orderrietkerk2004science ; rietkerk2008tree ; Meron2012ecomod . On the other hand, it is well established by model studies and confirmed by empirical observations that patch-scale biomass-water feedbacks can lead to regular landscape-scale vegetation patterns, and that periodic gap patterns, which highly resemble fairy circle patterns, can appear from uniform vegetation in response to water stressHardenberg2001prl ; Rietkerk2002an ; Gilad2004prl ; Borgogno2009gr ; Deblauwe2008geb . Indeed, recent detailed comparisons of fairy circle remote sensing data and model gap patterns show high similarity in static properties, such as hexagonal symmetry and spatial correlations of fairy circle size and distancegetzin2015adopting .

In this paper, we combine a theoretical analysis of a vegetation model, fitted to the biotic and abiotic conditions of the NFC ecosystem, with an empirical data analysis, to account for fairy circle dynamics in the NamibRand Nature Reserve from the years 2004 to 2013, thereby accomplishing two goals. The first goal is to substantiate the view of fairy circles as a pattern-formation phenomenon by complementing the current statistical evidencetlidi2008lnp ; cramer2013namibian ; getzin2015adopting with evidence based on the dynamics of the NFC ecosystem. The second goal is to demonstrate the feasibility of gradual regime shifts in the NFC ecosystem as cascades of unidirectional transitions across hybrid states.

I Model

The model we use is based on the Gilad et al. vegetation modelGilad2004prl ; Gilad2007jtb , which captures three different mechanisms of vegetation pattern formationKinast2014prl . However, applying the model to the NFC ecosystem (sandy soil, confined root zones) results in a simplified model describing the dynamics of aboveground biomass ($B$) and soil-water ($W$) areal densities. The simplified model captures a single pattern-forming feedback associated with the high rate of water uptake by the perennial grasses and the fast soil-water diffusion (relative to biomass expansion) in sandy soils. This mechanism leads to higher soil-water content in the bare fairy circles as compared with the vegetation matrixKinast2014prl , in agreement with reported observationsjuergens2013biological ; cramer2013namibian . The model equations read:

[TABLE]

In the biomass equation, $\Lambda$ is the biomass growth rate coefficient, $K$ is the maximal standing biomass, $E$ is a measure for the root-to-shoot ratio, $M$ is the mortality rate, and $D_{B}$ represents the seed-dispersal or clonal growth rate. In the soil-water equation, $P$ is the precipitation rate, $N$ is the evaporation rate, $R$ is a dimensionless factor representing a reduction of the evaporation rate due to shading, $\Gamma$ is the water-uptake rate coefficient, and $D_{W}$ is the effective soil-water diffusivity in the lateral ($X,Y$) directions. The model parameters used in this study are different than those used by Getzin et al. getzin2015adopting , and were estimated from published data. We refer the reader to the SI Appendix for more details on the derivation of the model and to the Methods and Materials section for details on the parameter estimations and values.

II Results

We first study the stationary solutions, in one spatial dimension (1d), along the precipitation axis. Fig.1A shows uniform solutions representing bare soil and uniform vegetation and periodic solutions representing gap patterns. Within the bistability range of uniform vegetation and gap patterns, there is a subrange with many stable hybrid states, consisting of confined domains of uniform vegetation in an otherwise periodic vegetation pattern (and vice versa), (inset B). The figure further shows 1d spatial profiles of a periodic gap pattern (panel C) and of a hybrid state that describes a periodic pattern with one missing gap (panel D). The corresponding 2d patterns are shown in panels E and F, and similar fairy circle patterns (obtained from satellite images) are shown in panels G and H. Similarly to other models of dryland vegetation Borgogno2009gr ; zelnik2013regime , the model also exhibits spotted and striped patterns (see the SI Appendix), for lower values of precipitation.

Next, we use the model, along with high-resolution satellite images of the NamibRand area from 2004 to 2013, to suggest that fairy circle “birth” and “death” processes (appearance and disappearance of gaps)tschinkel2012life correspond to transitions between hybrid states. Figure 2A shows the birth of a fairy circle following a drought in 2007, while Fig. 2C shows the death of a fairy circle following a spate (a period of excess precipitation) in 2008. Figure 2B and 2D show the same behaviors in model simulations, starting with initial conditions that mimic the initial fairy circle patterns in 2004 and applying precipitation downshift and upshift to mimic the drought and spate. The precipitation shifts were chosen to be strong enough to take the system out of the subrange of hybrid states (in opposite directions), thus inducing hybrid-state transitions.

In addition to the appearance and disappearance of gaps, Fig. 2 also shows changes in gap sizes; a moderate drought expands the gaps, while a spate contracts the gaps. To gain a deeper understanding of these dynamics, we correlated fairy circle data that were extracted from satellite images with rainfall data from several meteorological stations in the NamibRand region, and compared the observed correlations with those predicted by the model. Figure 3 shows the cross-correlations between the accumulated rainfall over different time periods (integration time) and the average fairy circle size (panel A) and the fairy circle number (panel B). Analogous cross-correlations using model simulations (panels C and D) show a qualitative agreement. The simulation results are based on averaging over 100 regions and long time series, and therefore, they show smooth curves. Panels E and F show the same quantities for the model simulations but averaged over fewer regions and using shorter time series. The curves in these panels resemble the curves for the real data (panels A and B). The stronger correlations in the simulations (compared with those in the field) between the gap size and the precipitation is due to the fact that the simulated system is homogeneous, while the field is heterogeneous and also because the simulated precipitation is not intermittent as it is in the field. The results highlight the different time scales for changes in the size and in the number of fairy circles. For fairy circle size, the correlation magnitude is maximal for a relatively short integration time (of a few years), while for fairy circle number, the correlation magnitude does not decrease with time throughout the time range of available data. This suggests that the timescale for changes in the fairy circle size is much shorter than the time scale for changes in the fairy circle number.

We turn now to the question of gradual regime shifts and ask whether climate variability in the form of repeated short droughts can induce a gradual regime shift. To this end, we simulated the model equations, starting with a uniform vegetation state that is locally disturbed to form a few gaps, and using precipitation downshifts that take the system outside the hybrid-state subrange periodically in time. As Fig. 4A shows, such a precipitation regime drives the system towards a nearly periodic patterned state (rightmost panel) by successive shifts to hybrid states of larger gap numbers. The opposite transition is shown in Fig. 4B, in which a series of short spates drives the system from a nearly periodic gap pattern, through a multitude of hybrid states, to a uniform vegetation state. Systems initially in a uniform vegetation state (see the SI Appendix) or a perfectly periodic patterned state (Fig. 4C), exposed to the same series of droughts or spates, show only minor fluctuations in the biomass density and remain practically unchanged.

III Discussion

We regard the model results described above and their agreement with field observations as compelling evidence for the view of fairy circles as a pattern-formation phenomenon and for the existence of hybrid states in the NFC ecosystem. This evidence is based on the following grounds: Spatial properties – the model, adjusted to the particular context of the NFC ecosystem, predicts gap patterns that resemble the spatial distribution of fairy circles in the field, as Fig. 1 shows. Moreover, the model predictions agree with recent field observations of higher soil-water content within the fairy circles, relative to vegetated domainsjuergens2013biological ; cramer2013namibian . Dynamics – fairy-circle birth and death processes can be reproduced with model simulations by temporal escapes outside the existence range of hybrid states (Fig. 2), supporting the interpretation of these processes as hybrid-state transitions. Correlations with precipitation – significant correlations between the rainfall and the number and size of fairy circles indicate the importance of biomass-water interactions in the dynamics of the NFC ecosystem. These correlations and the time scales associated with them, short for size variations and long for number variations, are in agreement with the model predictions (Fig. 3). The time-scale difference can be understood using pattern formation theory; gap-size changes are related to the convergence to a steady state, while gap-number changes are related to the dynamics of fronts that separate domains of uniform vegetation and gap patternsGandhi2015prl , structural and wavelength changes of the patterned state, and other nonlinear processes.

This evidence does not rule out alternative hypothesesVlieghe2015ecol_ento ; juergens2015ecography , but these hypotheses cannot be considered as satisfactory explanations of the fairy circle phenomenon unless they uncover the local feedbacks that are responsible for the emergence of large-scale order and account for the correlation of fairy circle dynamics with rainfall variability. The model suggested here requires additional evidence in order to fully uncover the origin and dynamics of the NFC. First, measurements of the soil-water diffusivity and the dependence of the soil-water extraction rate on the vegetation density are required in order to support the mechanism suggested in our model. Second, the precipitation fluctuations in the field are very large and, according to the simplified model, are expected to drive the system towards one of the uniform states. We believe that accounting for the spatial heterogeneity of the ecosystem, i. e., topography, soil heterogeneity and other factors, would result in more robust patterns in the model. Third, the exact biological origin of hybrid states in the NFC system remains unclear. Controlled field experiments of the FCs dynamics, combined with the mathematical framework suggested here, will be able to provide conclusive evidence for the hypothesis proposed here and reveal the relevant biological mechanisms. In particular, experiments in areas with different soil textures, topography, and plant species, together with manipulations of the precipitation, are expected to better clarify the relative importance of specific plant physiology and ecosystem characteristics in shaping the observed NFC dynamics.

The results presented here are not limited to the NFC ecosystem and may be observed in other pattern-forming systems. The appearance of hybrid states in systems exhibiting a bistability of uniform and patterned states implies that transitions, from uniform vegetation to a gap pattern and vice versa, may take place through successive hybrid-state transitions driven by periodic escapes outside the existence range of hybrid states. Dynamics of this type allow for spatially extended ecosystems to exhibit gradual rather than abrupt regime shifts.

IV Model Analysis

Vegetation pattern formation in this simplified model results from a finite wavenumber (Turing) instability of the uniform vegetation state. This instability requires strong water uptake, quantified by the parameters $\Gamma$ and $E$, and fast soil-water transport, characterized by $D_{W}$, relative to biomass dispersal, quantified by $D_{B}$Kinast2014prl . The high hydraulic conductivity, typical of sandy soils, acts in favor of this condition, but field observations point towards the possible role of an additional factor–the large underground termite populations found in the NFC ecosystemjuergens2013biological . The underground network of channels that the termites form may increase soil-moisture transport in the lateral directions, leading to higher effective values of $D_{W}$. Studying this possible termite effect is not only significant for understanding fairy circle formation, but also for understanding regime shifts, as the presence of termites may affect the existence and stability ranges of the uniform and patterned vegetation states, as well as those of the hybrid states.

The bifurcation diagram in Fig. 1 A,B was calculated using a numerical continuation method (AUTO software). The dynamics of the system was simulated using the pseudo-spectral method with periodic boundary conditions in Matlab, in either one or two dimensions. The simulations in Fig. 2 were made using initial conditions of steady states of the system, each derived from a region of the 2004 satellite image. This was done by taking the relevant region in the original image, segmenting it into vegetation and bare soil areas, and setting the values of the biomass and soil-water areal densities to those corresponding to the uniform solutions (uniform vegetation and bare soil, respectively). The system was then integrated forward in time, with a constant value of precipitation rate, $P$, until a steady state was reached.

V Parameter Estimations

In all model simulations, we used the following parameter values: $E=7[m^{2}/kg]$, $K=0.4[kg/m^{2}]$, $M=10.5[1/yr]$, $N=15[1/yr]$, $\Lambda=0.9[(m^{2}/kg)/yr]$, $\Gamma=12[(m^{2}/kg)/yr]$, $R=0.7$, $D_{B}=1.2[m^{2}/yr]$, $D_{W}=150[m^{2}/yr]$. These values were either taken from data published in the literature or estimated using such data, as described below. The precipitation rate, $P$, was varied within the range $50<P<150[mm/yr]$, which is similar to the range of annual rainfall in the NFC ecosystem.

The maximum standing biomass, $K$, was estimated based on maximal biomass densities reported for semi-arid savannah grassesnepolo2012short . The mortality rate, $M$, was estimated based on the mortality data of Stipagrostis Uniplumis as a representative example of an abundant perennial grass in the NFC ecosystemzimmermann2010grass . The reported mortality rate was measured under normal conditions (i.e., including all the growth factors). The mortality rate in the model represents the natural mortality alone, and therefore, we used a higher mortality rate. The biomass diffusion parameter, $D_{B}$, was estimated based on published clonal growth and seed dispersal datacain1997clonal . The value for the shading factor, $R$, was inferred from data on reduced evaporation by tree canopieswallace1999modelling . No relevant literature was found for the root-to-shoot parameter, $E$, for the perennial grass; the chosen value is similar to that used in an earlier studyGilad2007jtb for woody dryland species. The evaporation rate, $N$, was estimated using the relation $N=P/W$ for a uniform bare soil state, where the soil-water density, $W$, was determined from water-content sampling data inside a fairy circlejuergens2013biological . The water-uptake rate coefficient, $\Gamma$, was inferred from differences in the levels of soil water within and outside the fairy circles, and from biomass measurementswesuls2010perceptions . The biomass growth parameter, $\Lambda$, was determined from estimates of the ratio of transpiration to dry biomass for tropical grassErickson2012Water and the previously estimated value of $\Gamma$.

Finally, the value of the soil-water diffusivity, $D_{W}$, was estimated using the Van Genuchten model for water retention curves vanGenuchten1980ClosedForm . Using parameters corresponding to sandy soilVGparams and the values of soil-water content that correspond to those found in the fieldjuergens2013biological , we estimate $D_{W}$ to be in the range $0.1-100[m^{2}/year]$. We used a value that slightly exceeds the upper limit of this range to account for the possible enhancement of lateral water transport due to the network of underground channels formed by termite populationsjuergens2013biological . We emphasize that the model results reported here are robust and can be obtained with many different sets of parameters (see the SI Appendix (Sensitivity Analysis)).

VI Earth Data Analysis

We analyzed high resolution satellite images taken over the period of 2004-2013 in the NambiRand Nature Reserve. The snapshots shown in Figs. 1 and 2 were processed from four images, captured in 2004, 2008, 2010 and 2013, whose spatial resolution is as high as $0.5[m]$. The rgb data shown were derived using Pan-Sharpening with ERDAS software when rgb was available (2004, 2010 and 2013), and false-color was used for the 2008 image.

The data used for Fig. 3 were taken from monthly precipitation data and from measurements of fairy circles features using satellite images. The precipitation data were derived from interpolations of both the readings of meteorological stations within the NambiRand Nature Reserve (for the years 2000-2012) and of the NOAA data center (for the years 1990-1999) for two locations in Namibia. The numbers and sizes of the fairy circles were estimated from eight satellite images, the four previously described (from 2004, 2008, 2010 and 2013), and four others (from 2005, 2009, 2011 and 2012) taken from Google Earth data. From each of the eight images, four regions of $180[m]x180[m]$ were chosen and analyzed. A proprietary segmentation algorithm, involving morphological operations and thresholding with a CLAHE algorithm, was used to identify each fairy circle. The identification results were then reviewed manually to minimize errors caused by the automatic process, and in total, about 800 fairy circles were identified and measured for each time point.

The cross-correlation measures shown in Fig. 3 were calculated by comparing the accumulated precipitation, $\psi$, with either the number of fairy circles, $\alpha$, or the average size of a fairy circle, $\beta$. The accumulated precipitation is defined as:

[TABLE]

Here, $P_{k}$ is the precipitation at time $k$, $L$ is the lag time and $S$ is the integration time, with all times given in number of years. In order to fit the rainy season in the field, we considered the annual precipitation between November 1 and the following October 31. The lag time, $L$, was set to zero (i.e., considering the accumulated rain during a period $S$ ending at the end of the previous year) for the fairy circle size and $L=1$ for the fairy circle number; these values yield the most significant correlations between the precipitation and the fairy circle dynamics. The cross-correlation for the fairy circle number $C_{\alpha}$ is defined as:

[TABLE]

$\bar{A}$ stands for the average of $A$ ($\bar{A}\equiv\frac{1}{N}\sum\limits_{i=1}^{N}A_{i}$), and $N$ is the number of time points, with the value of $8$ for the observation data, and either $8$ or $100$ for the simulation data. The cross-correlation for the fairy circle size, $C_{\beta}\left(S\right)$, was calculated in the same manner.

Acknowledgements.

We wish to thank Nils Odendaal, the Chief Executive Officer of the NamibRand Nature Reserve, for providing us the precipitation measurements of various meteorological stations within the reserve, and we also thank Elad Eizner and Noa Levi-Ohana for their help with the segmentation and analysis of the satellite data. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant number [293825], and from the Israel Science Foundation under grant number 305/13.

Appendix A Derivation of the simplified model

We present the derivation of the simplified model that we use and analyze in the main text. Our starting point is the vegetation model of Gilad et al. Gilad2007jtb that includes three pattern-forming feedbacks. Using the characteristics of the fairy circles ecosystem, we explain the assumptions and approximations that allow us to replace the integrals of the root-augmentation feedback with simpler algebraic terms and to decouple the aboveground water dynamics from the soil-water and vegetation dynamics.

The dryland vegetation model of Gilad et al. Gilad2007jtb describes the coupled dynamics of the areal densities of vegetation biomass (B), soil water (W) and surface water (H), all having the dimension of mass per unit area. The dynamics is described with a temporal resolution that is smaller than the typical time scale for changes in the biomass (dictated by the biomass growth rate, mortality and dispersion rate) but is large enough to allow averaging the intermittent nature of the precipitation. The spatial scale resolution is larger than the scale of a single plant (in order to allow the description of a continuous vegetation density) but smaller than the typical patch size. Restricting our interest to flat terrains, the model reads:

[TABLE]

where

[TABLE]

The terms $G_{B}$ and $G_{W}$ involve integration over a root kernel ($G$) that represents the spatial extent of the root zone in the lateral directions. Applying the model to perennial grasses, e.g. Stipagrostis Ciliata (common in the fairy circle ecosystem), which form patterns with a typical length scale of  $10m$ picker2012ants and have a root girth of approximately $0.5m$ Midgley1999Book , allows us to assume that the kernel function is much narrower than the biomass and soil-water distributions. Under this condition, we can approximate the kernel by a Dirac delta function. Formally, this is done by taking the limit $S_{0}\rightarrow 0$ in Eq. 3c, where $S_{0}$ represents the lateral root-zone size of a seedling. Following this approximation, we may replace the integrals of Eqs. (3c) by the algebraic forms:

[TABLE]

The fairy circle ecosystem consists of sandy soil. This soil type is characterized by a high rate of surface water infiltration that is comparable to the infiltration rate in vegetated soil. To account for the absence of a significant infiltration contrast between bare and vegetated soil, we set $f=1$. The infiltration rate then becomes a constant, $I=A$, independent of the biomass $B$, and the equation for the surface water variable $H$ decouples from those for $B$ and $W$. This equation has a single stationary uniform solution, $H_{0}=P/I$, which is always linearly stable. Since $H$ is the fastest variable, we can assume that, on the much slower time scales over which $B$ and $W$ significantly change, it has already equilibrated at $H_{0}$. Inserting the solution $H=H_{0}$ into the equation for $W$, we obtain the two-variable model:

[TABLE]

In the biomass equation, $\Lambda$ is the biomass growth rate coefficient, $K$ is the maximal standing biomass, $E$ is a measure for the root-to-shoot ratio, $M$ is the mortality rate, and $D_{B}$ represents the seed dispersal or clonal growth rate. In the soil-water equation, $P$ is the precipitation rate, $N$ is the evaporation rate, $R$ is a dimensionless factor representing a reduction of the evaporation rate due to shading, $\Gamma$ is the water-uptake rate coefficient, and $D_{W}$ is the effective soil-water diffusivity in the lateral ($X,Y$) directions, assumed to be a constant, independent of the state variables, that represents linear diffusionGilad2007jtb . We refer the reader to an earlier publicationGilad2007jtb for additional information about the original model. The model considered here (Eqs. (5)) predicts, in addition to the states shown and discussed in the main paper, the appearance of spotted and patterned states under low precipitation rates as shown in Fig. 5.

Appendix B Sensitivity Analysis

In order to test the sensitivity of the results (presented in the main text) to the set of parameters, we first write the model in a dimensionless form. Thereby, we reduce the number of parameters and simplify the analysis.

Equations (5) describe the dynamics of the aboveground biomass ($B$) and the soil water ($W$) in a dimensional form. Translation of the model equations to a dimensionless form is achieved by rescaling the state variables $B,W$ and the space and time coordinates as follows:

[TABLE]

In terms of these dimensionless quantities, the model reads:

[TABLE]

[TABLE]

The dimensionless parameters here are related to their dimensional counterparts by the following relations:

[TABLE]

We are thus left with only six parameters, for which we use $p$ as the main bifurcation parameter, and as such, its value changes throughout our analysis. The basis for the main results is the existence of stable hybrid states; therefore, we tested the conditions under which they occur and their existence range (range of precipitation rate, $p$). We changed the other five parameters, one at a time, within $\pm 10\%$ of their value. We found that in all cases, the hybrid states exist, with their existence range slightly changing. The relative changes in their existence range are shown in Table I.

As can be seen, the parameter $\lambda$ has the strongest effect on the snaking range. The parameters $\eta$ and $\delta_{w}$ control the pattern-forming feedback in our model Kinast2014prl . Therefore, we present a more thorough analysis of their effect on the snaking range. In Fig. 6, we show the snaking range as a function of either $D_{W}$ (left panel) or $E$ (right panel). Fig. 7 shows two examples of the bifurcation diagram for a relatively big change in one of these parameters.

In addition to the sensitivity of the snaking range to the parameters the size of the gaps and the distance between neighboring gaps are also affected by the parameters. The histogram of the FC sizes, in the regions for which the satellite images were analyzed, is shown in Fig. 8.

In Fig. 9 we show the average size of the FCs and the average distance between them as predicted by our model for constant precipitation rate of $P=102[mm/yr]$ and for a range of the feedback control parameters $D_{w}$ and $E$. It is shown that there is a range of values of these parameters that correspond to the FCs size observed in the field. The range of possible values for these parameters, which have not been measured directly in the field, may be estimated from these graphs.

Appendix C Predicted response of uniform vegetation to a series of droughts

In the main paper, we showed that the response of vegetation to a series of droughts or spates depends on the initial condition. To complement the information we present in Fig. 10 the response of a uniform vegetation to a series of droughts.