Fertility Island Formation and Evolution in Dryland Ecosystems

Vast dryland regions around the world are affected by the encroachment of woody vegetation, with important environmental and economical implications. Grassland-to-shrubland conversions are often triggered by disturbance of grassland vegetation, and the consequent formation of barren areas prone to erosion-induced nutrient losses. Inhibition of encroachment by erosion-induced depletion of soil nutrients contributes to the emergence of highly heterogeneous landscapes with shrub-dominated fertility islands surrounded by nutrient-poor bare soil. Here, we develop a process-based simplistic model thataccounts for the two competing processes of resource depletion and shrub encroachment by a non-linear diffusion mechanism. The proposed model is able to generate stable vegetation patterns with the same statistical properties as those observed in areas with well-developed fertility islands. We also show how a subsequent disturbance of shrubland vegetation can shift the dynamics toward states with smaller vegetation biomass. The process of land degradation may then occur through a number of irreversible intermediate transitions associated with losses in ecosystem function.


INTRODUCTION
The conversion of semiarid grasslands into shrublands (e.g., Scholes and Archer 1997) has been observed in many regions of the world, including North America (e.g., Archer 1989, van Auken 2000), South America (e.g., Adamoli et al. 1990), Africa (e.g., Scott 1966, Hudak andWessman 2001), and Australia (e.g., Sharp and Whittaker 2003).This process is often associated with the formation of heterogeneous landscapes with unfertile bare areas separated by nutrient-rich soil plots colonized by shrubs (Charley andWest 1975, Schlesinger et al. 1990).Known as "fertility islands" (Schlesinger et al. 1990, 1996, Schlesinger and Pilmanis 1998, van Breemen and Finzi 1998, Reynolds et al. 1999), these vegetated, fertile soil patches are the result of a landscape-scale redistribution of resources, with losses of nutrients and ecosystem function in the unvegetated areas, and accumulation of fertile soil particles beneath the shrub canopies (Schlesinger et al. 1990).In this paper, we focus on island-shaped vegetation patterns that originate from this specific process of soil redistribution; we do not address patterned landscapes caused by other mechanisms.
Both biotic and abiotic processes are commonly invoked to explain the autogenic fertilization of shrub-dominated areas, including plant uptake from the surrounding soils, nitrogen fixation by desert shrubs, animal activity, the funnelling of nutrientrich stemflow along shrub branches, and the ability of shrub canopies and soils to catch nutrient-rich soil particles transported by winds or water (e.g., Schlesinger and Pilmanis 1998).The erosion-driven degradation of intercanopy soils and the consequent emission and transport of nutrient-rich atmospheric dust are known to have important implications on regional and global climate (Ramanathan et al. 2001, Rosenfeld et al. 2001, Prospero and Lamb 2003), biogeochemical cycles (Duce and Tindale 1991, Swap et al. 1992, Okin et al. 2004), desertification (e.g., Schlesinger and Pilmanis 1998, Reynolds et al. 1999, Maestre et al. 2006), and human health (Dockery and Pope 1994).http://www.ecologyandsociety.org/vol13/iss1/art5/However, despite the recognized importance and widespread occurrence of grassland-to-shrubland conversions (van Auken 2000) and the growing interest in patterning mechanisms (e.g., Lefever and Lejeune 1997, Klausmeier 1999, von Hardenberg et al. 2001, HilleRisLambers et al. 2001, Barbier et al. 2006), to date no specific mathematical models exist that are capable of explaining the emergence of fertility islands.The lack of a process-based mathematical framework limits the ability to investigate possible changes in ecosystem structure and function under different climate and disturbance regimes, thereby preventing assessment of the stability and resilience of shrubland ecosystems.To this end, we propose a process-based simplistic model that accounts for the mechanisms typically associated with fertility island dynamics.Spotted vegetation can be obtained with a number of models invoking different mechanisms of pattern formation; these models are generally based on selforganized dynamics conducive to symmetrybreaking instability, including Turing-like instability (e.g., Rietkerk et al. 2002, Meron et al. 2004, Guttal and Jayaprakash 2007), short-range cooperative and long-range inhibitory interactions (e.g., Lefever and Lejeune 1997, Lejeune et al. 2002, D'Odorico et al. 2006), and instability induced by differential flow rates in an activator-inhibitor system (e.g., Klausmeier 1999, Okayasu and Aizawa 2001, Gilad et al. 2007).However, these mechanisms are different from the key processes characteristic of fertility island dynamics (e.g., Charley andWest 1975, Schlesinger et al. 1990).Here, we propose a new pattern formation model based on 1) facilitation between vegetation and the availability of soil resources, and 2) long-range spatial interactions accounted for by a long-range diffusion term.Field observations have shown that both these processes contribute to fertility island dynamics (Schlesinger et al. 1990, Reynolds et al. 1999).The proposed modeling framework is able to generate (non selforganized) vegetation patterns similar to those commonly observed in landscapes where the formation of fertility islands has been well documented (e.g., Okin et al. 2001).This framework can provide useful tools for the management of desert shrublands.

THEORETICAL FRAMEWORK
We focus on the dynamics of fertility island formation that have been observed subsequent to the destabilization of arid grasslands such as those in the American Southwest (Schlesinger et al. 1990).The decline of these grasslands is often triggered by exogenous drivers (e.g., Schlesinger et al. 1990, Reynolds et al. 1999)-e.g., overgrazing or changes in fire management-whereas the formation of heterogeneous distributions of resources and vegetation has been argued to be induced by processes and feedbacks that are internal the system (e.g., Schlesinger andPilmanis 1998, Reynolds et al. 1999).Here, we deal with these processes without investigating in detail the causes of grassland destabilization.We just assume that, after the disturbance of the grass cover, a relatively uniform distribution of soil nutrients remains exposed to the erosive action of wind and water.We model the main processes capable of leading to the emergence of characteristic patterns known as "fertility islands," starting from these relatively uniform initial conditions.
The formation of fertility islands results from the interaction between two competing processes triggered by the disturbance of semiarid grasslands: 1) wind or water erosion in areas with thin vegetation canopy, and 2) the encroachment of shrub vegetation facilitated by the disappearance of grasses (i.e., of possible competitors), the consequent weakening of fire dynamics (van Auken 2000), and the autogenic fertilization of shrubdominated areas (Schlesinger andPilmanis 1998, Reynolds et al. 1999).Both nutrient loss from intercanopy soils and shrub auto-fertilization convert the relatively uniform initial configuration of grassland vegetation (Reynolds et al. 1999) into a highly heterogeneous distribution of resources, with a mosaic of nutrient-rich shrub-colonized soil patches (i.e., the fertility islands) bordered by nutrient-poor bare soil.This mechanism of resource redistribution acts through positive feedbacks with vegetation dynamics, as the loss (gain) of nutrients is facilitated by the decrease (increase) in grassy (woody) vegetation.
A broad range of scales is involved in the process of fertility island formation.In fact, soil erosion, and shrub establishment occur at the plot-to-field scale (i.e., ~1-10 m 2 ).Soil crusting, processes contributing to vegetation autogenic fertilization (e. g., nitrogen fixation, stemflow, and water and nutrient uptake by lateral roots), and the colonization of soil plots by neighboring shrubs are all examples of local dynamics, which involve short-range scales (i.e., within the range of a few meters).On the other hand, soil erosion caused by http://www.ecologyandsociety.org/vol13/iss1/art5/wind or surface wash, and mid-to-long range transport of seeds by wind or animals are examples of processes occurring at the field-to-landscape scale (i.e., >10 m 2 ).In the following, we model these complex interactions between resource redistribution and shrub encroachment through a system of two coupled differential equations, which express the dynamics of shrub biomass, V, and limiting resources, R. Both V and R are normalized variables ranging between 0 and 1.

Shrub Dynamics
The dynamics of shrub biomass, V, are expressed as the sum of two terms accounting for local processes and spatial interactions.The local dynamics are modeled as a logistic growth (e.g., Murray 1989), whereas the spatial interactions are expressed in the form of an integral term accounting for the effect of the surrounding vegetation on the encroachment rate.Thus, at any point x = {x,y} in the two-dimensional domain, Ω, the dynamics of shrub vegetation are modeled as: (1) where t is time, α is a parameter that modulates the growth rate, V cc is the ecosystem carrying capacity, i.e., the maximum amount of plant biomass sustainable with the available resources, D(R,V) is the diffusion coefficient, and k(|x' -x|) is a (positive) kernel function normalized in a way that at any point We express V cc as a function of R to account for the control exerted by resources on the rate of shrub growth.Because both unvegetated/nutrient-poor and shrub-dominated/fertile soil patches are stable in landscapes affected by fertility island formation, the dependence of the carrying capacity on the available resources needs to be capable of inducing bistability in the dynamics of V and R. To this end, the carrying capacity needs to be expressed as a nonlinear (quadratic) function of R (2) It will be shown that when the levels of available resources, R, are relatively low, vegetation growth is unsustainable and the system tends to the stable steady state V = R = 0. On the other hand, for relatively high values of R, both shrubs and resources tend to the stable steady state V = R = 1.We will return later to the properties of these bistable (local) dynamics (see "Stability Analysis").
The second term on the lefthand side of Eq. ( 1) accounts for the spatial dynamics of shrub encroachment at any point x in Ω.The rate of these dynamics depends both on (a) the shrub biomass existing in the surrounding area, and (b) the conditions (i.e., resources and vegetation biomass) existing at that point.The integral term in Eq. ( 1) models the dependence on the biomass of the surrounding vegetation.The kernel function, k(|x'-x|) weights the effect of the neighboring vegetation (i.e., shrub biomass existing at any nearby point, x') on shrub encroachment at point x, with k being a decreasing function of the distance between x and x'.The rate at which k decreases with the distance determines the range of the spatial interactions occurring in the system.
The integral formulation provides a more general representation of vegetation encroachment than the "classic" linear diffusion term (i.e., \nabla^2 V) used in a number of ecological models (e.g., Okubo 1980).In fact, diffusion leads to an unrealistic reduction of V in vegetated areas bordered by plots with lower vegetation density.To avoid this effect, the decrease in V is generally compensated by a local growth in the form of a Fisher equation (i.e., linear diffusion combined with logistic growth; Murray 1989).The integral form used in (1) does not induce any decrease in V within vegetated soil patches surrounded by areas with lower vegetation density.Moreover, it can be shown that, when the integral term is approximated by a series expansion in the neighborhood of x (i.e., for low values of |x'-x|), the low-order terms of this approximation include the diffusion operator (e.g., Murray 1989, p. 244), whereas higher-order terms account for longerrange interactions.Thus, the integral form in (1) is a generalization of a diffusion process, as it is capable of accounting both for diffusion (i.e., shortrange encroachment) and long-range interactions.In what follows, we will denote the spatial http://www.ecologyandsociety.org/vol13/iss1/art5/interactions modeled by (1) as "long-range" diffusion (Murray 1989).Because k is taken as a positive or null function, the effect of surrounding bushes is to enhance the rate of shrub encroachment (i.e., facilitation), whereas the effect of competition with neighboring plants is not considered in this study.In fact, the role of competition-facilitation dynamics in the process of vegetation pattern formation has already been investigated by other authors (Lefever andLejeune 1997, Lejeune et al. 2002) using a long-range diffusion model.Here, we use k ≥ 0 to show how fertility islands can emerge from non-linear interactions between shrub vegetation and resources without invoking plant competition.Therefore, the particular form of the kernel function can affect the shape and size of islands, but not the key mechanisms of pattern formation.
The dependence of the encroachment rate on the vegetation and resources existing at the point x is accounted for by expressing D as a function of both R and V. Due to this dependence, Eq. ( 1) becomes a non-linear long-range diffusion process.This nonlinearity plays a crucial role in the process of shrub encroachment in that it accounts for local limitations to shrub colonization due to resource availability.Because the encroachment rate at point x is likely to increase with an increase in the level of resources available at the same point, we take D as proportional to R.Moreover, the encroachment rate at point x is expected to decrease as V tends to the carrying capacity V cc .Thus, we express the diffusion coefficient as: (3) where D 0 is a parameter that accounts for the strength of the diffusion term in Eq. ( 1).
Non-local processes contributing to shrub recruitment (e.g., seed dispersal due to animals, wind, or surface runoff), are accounted for by modeling seedling dispersal as a Poisson process (e. g., Cox and Miller 1977) in space and time, with shrub seedlings occurring with probability, λ, per unit area and time.

Resource Dynamics
The temporal dynamics of R are modeled as a deposition-erosion process, with deposition rate linearly dependent on shrub biomass (due to autogenic fertilization) and erosion rate proportional to the existing resources and to the area, (1-V), left with no vegetation cover, (4) where β and γ are two (positive) coefficients expressing the relative importance of the deposition and erosion processes.
The term ε (<< 1) in Eq. ( 4) accounts for the net input of resources associated with atmospheric deposition.This process might significantly affect the resource budget by preventing a complete loss of resources from bare soils.However, deposition does not contribute to the mechanisms of pattern formation investigated by this paper, in that it does not induce any spatial heterogeneity in the distribution of vegetation and resources.The only way atmospheric deposition can favor fertility island formation is through its feedback with vegetation, i.e., by enhancing the accumulation of nutrient-rich soil particles beneath shrubs.However, this effect is already accounted for by the first term on the righthand side of Eq. ( 4).Thus, in what follows, we will simply take ε = 0.

Stability Analysis
The emergence of both fertility islands and vegetation patterns is investigated through the stability analysis of the homogeneous steady states of the system.To this end, we first rescale time as t' = αt and rewrite Eqs.Before investigating the effect of spatial interactions, we consider the local (i.e., temporal) dynamics expressed by Eqs. ( 5)-( 6) with D 0 = 0.In this case, the system has three steady homogeneous states: (V,R) 1 = (0,0); (V,R) 2 = (γ 2 /(β-γ) 2 , γ/(β-γ)); and (V,R) 3 = (1,1).Because both V and R range between 0 and 1 (i.e., are normalized with respect to their maximum values), (V,R) 2 is a possible state of the system only if β > 2γ.The stability of (V,R) 2 and (V,R) 3 can be assessed through classical linear stability analysis (Cross and Hohenberg 1993, Drazin 2002): it is found that (V,R) 2 is always unstable, whereas (V,R) 3 is stable.Because the eigenvalues of (V,R) 1 have vanishing real parts, this state is a degenerate (i.e., non-hyperbolic) point, and its stability needs to be investigated numerically, in that it is completely determined by the non-linear terms of the dynamics (Hartman-Grobman theorem (Guckenheimer and Holmes 1983)).It is found that this steady state is always stable.Thus, the system is bistable and converges to either one of its two stable states, (V,R) 1 and (V,R) 3 , depending on the initial conditions.The dotted line in Fig. 1a marks the boundary between the attraction basins of these two stable points.With initial conditions corresponding to points located below (above) that line, the system tends to the stable state (V,R) 1 (or (V,R) 3 ).
To investigate the effect of the spatial term in Eq. ( 5), we repeat the previous analysis with values of D 0 ≠ 0. In this case, the system has the same three (homogeneous) steady states with the same stability conditions as before.However, as shown in Fig. 1a (dashed line), the boundary (or "potential barrier") of the attraction basins of the two stable states, (V, R) 1 and (V,R) 3 , changes with respect to the local dynamics (dotted line (D = 0)).When V is smaller than the carrying capacity (solid line in Fig. 1a) the size of the basin of attraction of ((V,R) 3 ) increases at the expense of that of ((V,R) 1 ), suggesting that spatial interactions favor the establishment of shrub vegetation and the attainment of stable vegetated conditions.In fact, because the diffusion coefficient is positive for V<V cc = R 2 (Eq.( 3)), in these conditions the spatial term in Eq. ( 5) favors the process of shrub encroachment.
Figure 1a suggests that, on nutrient-rich soils, the establishment and growth of isolate shrub seedlings can be observed (Fig. 1a, zone A).As resources are lost from unvegetated (or sparsely vegetated) areas, the system moves toward the boundary between the two basins of attraction in Fig. 1a: in these conditions, shrub establishment can still occur, although the growth rate decreases as R decreases.
As the point (V,R) representative of the state of the system crosses the potential barrier for the case D = 0 (dotted line in Fig. 1a), the local resources are unable to sustain the establishment and growth of isolated shrubs (Fig. 1a, zone B).Thus, the positive feedback between vegetation dynamics and resources leads the system to a stable unvegetated and unfertile state.However, in the presence of surrounding vegetation, the spatial dynamics can still allow for colonization of the bare soil by shrubs, because of the favorable environment (e.g., higher moisture and nutrient contents) provided by the nearby vegetation.These spatial dynamics favor the growth of shrub clumps even in environments that have undergone significant resource depletion.On the other hand, a further decrease in soil resources (Fig. 1a, zone C) shifts the system into the domain of attraction of (V,R) 1 for the case D≠ 0, indicating that, despite the existence of spatial interactions, shrub establishment is inhibited by resource limitations.At this point, fertility islands are unable to expand further and vegetation patterns remain "frozen."The landscapes resulting from these dynamics exhibit vegetated, nutrient-rich soil patches surrounded by unfertile bare soil.Other processes (e.g., fire dynamics, pest outbreaks, grazing) may undoubtedly contribute to the formation of this patchy landscape.However, this study aims to develop a minimalistic framework for modeling the key mechanisms that determine the formation of fertility islands, without invoking the role of complex interactions with the disturbance regime.
Thus, this analysis shows that the dynamics expressed by Eqs. ( 5)-( 6) exhibit two key features: 1) bistability induced by the non-linear dependence of the carrying capacity, V cc , on the resources, and 2) the formation of patterns of vegetation and resources induced by the mechanism of long-range non-linear diffusion.Unlike the short-range diffusion (i.e., <nabla> 2 V) these mechanisms are http://www.ecologyandsociety.org/vol13/iss1/art5/It should be noted that several existing models of pattern formation are able to generate landscapes with spotted vegetation resembling vegetation patterns associated with the formation of fertility islands.Recent studies were either based on the socalled Turing's (1952) instability (Klausmeier 1999, von Hardenberg et al. 2001, HilleRisLambers et al. 2001) or on neural models (Lefever andLejeune 1997, Lejeune et al. 2002).In the first case, patterns emerge as a result of the instability induced by at least two competing processes of linear diffusion having very different diffusivities (e.g., Murray 1989).Models based on Turing instability (Klausmeier et al. 1999, von Hardenberg et al. 2001) require at least two state variables to generate vegetation patterns, whereas in neural models, patterns emerge from the joint effect of short-range activation and long-range inhibition.Lefever and Lejeune (1997) used a neural model to explain the formation of vegetation patterns: in this model, short-range activation represented the effect of positive intra-specific interactions (facilitation) acting within a shorter range (e.g., subcanopy autogenic fertilization) than the inhibition induced by plant competition.In our study, we propose a http://www.ecologyandsociety.org/vol13/iss1/art5/model (Eqs.( 1)-( 2)) that is conceptually different from both frameworks.In fact, in our case, only one state variable diffuses and there are no explicit negative intra-specific interactions.Moreover, unlike Turing instability models, non-linear longrange diffusion is the key mechanism here for pattern formation.In fact, because shrub encroachment takes place only when sufficient soil resources are locally present, the diffusion rate of vegetation depends (non-linearly) on the soil resource, and diffusion is inhibited for lower values of available resource.

RESULTS AND DISCUSSION
The spatial and temporal dynamics of shrub vegetation and available resources expressed by Eqs. ( 5) and ( 6) generate vegetation patterns (Fig. 1) similar to those observed in regions where the existence of fertility islands has been welldocumented (e.g., Tiedermann and Klemmedson 1986).Figure 2(a, c) shows two model-simulated vegetation patterns that are almost indistinguishable from two images of shrublands from the American Southwest (Fig. 1b, d).The irregularity of the size and shape of vegetation patterns arises from the coalescence of adjacent vegetated areas and from the randomness both of initial conditions and of seedling occurrences.The elongation and anisotropic growth induced in real fertility islands by preferential directions in resource redistribution (Okin andGillette 2001, Schade andHobbie, 2005) can be accounted for by using a (possibly random) anisotropic diffusion in Eq. ( 5) for the encroachment process, e.g., adopting an asymmetric kernel function, k.
To quantitatively compare model results with data, we characterize the spatial configuration of shrub vegetation through the cluster size distribution of vegetated areas.Each cluster is defined as a set of connected vegetated sites according to a fourneighbor topology (e.g., Gould and Tobochnik 1996).For the two examples shown in Fig. 2, the model results compare favorably with the data (see Fig. 3), suggesting that Eqs. ( 5)-( 6) are able to generate landscapes of fertility islands with the same statistical properties as those observed in the field data.These patterns persist even once the system approaches its steady state.In fact, the nonlinearity embedded in the long-range diffusion term of Eq. ( 5) halts shrub encroachment when the resources available in the barren soils surrounding shrub-dominated areas become insufficient.Thus, the dependence of shrub diffusion on local resources prevents the complete encroachment of the landscape by shrubs, and stabilizes vegetation in configurations exhibiting well-defined fertility islands (Fig. 2).Notice that, in this model, vegetation patterns do not emerge as a result of prescribed spatially correlated heterogeneities associated with soils or micro-relief, but from the mutual inhibition between the two processes of resource depletion and shrub encroachment.Figure 4 shows the temporal dynamics of the spatial means (µ V and µ R ) and standard deviations (σ V and σ R ) of the fields of V and R, until the attainment of the steady-state conditions.The net increase in the standard deviation of R (Fig. 4) is in agreement with field observations reporting that the process of resource redistribution accompanying shrub encroachment generates highly heterogeneous spatial configurations of shrub biomass and resources, consistent with the notion of fertility island formation (Schlesinger et al. 1990, Bhark andSmall 2003).The process of land degradation leading to the loss of grassland vegetation and the formation of nutrient-poor barren areas is partially limited by the encroachment of shrubs, which locally prevent the complete depletion of soil resources and increase nutrient availability through mechanisms of autogenic fertilization.
The stability and resilience of this new state of the system (i.e., shrubland) is here investigated to assess its susceptibility to shifts to other states in response to changes in the disturbance regime (Scheffer et al. 2001, D'Odorico et al. 2005).Field experiments (e. g., Schlesinger andPilmanis 1998, Tiedermann andKlemmedson 1986) suggest that the removal of shrubs from fertility islands is generally followed by erosion-induced depletion of soil nutrients and consequent loss of ecosystem function.We simulate the effect of shrub uprooting by setting to zero (Fig. 3) the value of V at time t' = 250 (i.e., after the system has reached a steady state).The subsequent increase in erosion rates is found to deplete the resources accumulated in the fertility islands (Fig. 4b), although the regeneration of shrubland vegetation by seedling dispersal (from external sources) and growth may allow for a partial recovery of the system (Fig. 4), suggesting that the autogenic fertilization induced by regenerated shrubs is able to compensate for only part of the erosion-induced loss of soil resources.In fact, the system reaches a state with lower average vegetation density, as shown in Fig. 4a, b.These results, which are  Fig. 4. Model-simulated temporal evolution of the spatial means (µ V , µ R ) and standard deviations (σ V , σ R ) of the fields of V and R using the same parameters as in Fig. 2a.The removal of shrub vegetation at time t' = 250 is followed by a period of resource depletion concurrent to the establishment, growth, and encroachment of new shrub seedlings, and by a period of resource accumulation due to autogenic mechanisms induced by shrub vegetation.At a steady state, both shrub biomass and resources attain asymptotic levels lower (by about 43% in the case of V) than those existing before time t' = 250, whereas the asymptotic values of the coefficients of variation, σ V /V and σ R /R, increase.Similar analyses (not shown) obtained with the same parameters as in Fig. 2c lead to even lower post-disturbance steadystate values of V. http://www.ecologyandsociety.org/vol13/iss1/art5/obtained with a constant rate, λ, of seedling occurrence, do not account for possible positive feedbacks between λ and the shrub biomass.A decrease of λ in the post-disturbance dynamics would set the new asymptotic state to even lower values of V and R, including-in some casesdesert conditions.
To show the broad variety of configurations that can be generated by the model (Eqs.( 5)-( 6)) we investigate how the (spatial) standard deviation, σ V , of V (in the asymptotic state reached by the system) varies with the parameter D 0 (Fig. 5).For low values of D 0 only a sparse shrub canopy develops because shrub seedlings are unable to encroach and colonize their surroundings.As D 0 increases, the spatial term in the dynamics of V becomes stronger, thereby leading to the formation of shrub clumps, which increase in size as D 0 increases and eventually coalesce.Higher values of D 0 lead to the complete encroachment of the landscape by shrub vegetation and to a consequent decrease in the heterogeneity (i.e., σ V ) of the system (Fig. 5).Thus, depending on the parameters of the system, these dynamics can lead to different (stable) landscapes, including homogeneous unvegetated/vegetated soil, and patchy landscapes with fertility islands and bare soil.A qualitatively similar variety of patterns can be obtained either by increasing the parameters β and λ-i.e., by increasing the rate of resource deposition and of seedling dispersal-or by decreasing γ (i.e., the resource erosion).The environmental conditions, soil properties, and vegetation characteristics of a specific landscape will determine which of these mechanisms plays the dominant role in the formation of fertility islands.

CONCLUSIONS
The non-linear diffusion model developed in this study is capable of generating vegetation patterns through two competing mechanisms-namely, resource depletion and shrub encroachmentacting in conjunction with positive feedbacks between resource accumulation and shrub biomass.Field observations indicate that these patterns are the key mechanisms contributing to fertility island formation.The patterns generated by the proposed model are steady and have statistical properties comparable to those observed in real patterns from arid and semiarid shrublands.The same framework suggests that the shrubland state may have only a limited resilience and be prone to shifts to lower-resource/lower-biomass states, a process also known as "desertification."The stability of these subsequent states suggests the existence of a high degree of irreversibility in the process of land degradation.These possible scenarios could have a crucial impact on the management of desert shrublands.In this sense, field investigations are needed to better document the spatio-temporal evolutions of shrub-dominated landscapes, and to test the validity of the modeling results presented in this paper.

Fig. 2 .
Fig. 2. a) Example of vegetation patterns generated by the model in a numerical simulation on a square 500 x 500 grid.The parameters are β' = 0.9, γ' = 0.155, with random, uniformly distributed initial conditions for R (in [0.42, 0.57]) and V initially equal to 0. Seedling dispersal occurs at rate λ = 1.7 • 10 -3 pixel -1 per unit time, with each seedling having initial biomass, V 0 = 0.01.The kernel function, k, is different from zero (and positive) only in the eight nearest neighbors of each pixel, with k being a linearly decreasing function of the distance, and the sum of k calculated for each pixel and its nearest neighbors being equal to 1.The spatial coordinates are rescaled as x' = x √(α /D 0 ) and y' = y √(α /D 0 ).b) Aerial view (taken from Okin and Gillette (2001), Fig. 3c) of vegetation from the Jornada site, in the Chihuahuan Desert (New Mexico).c) Same as a) but with λ = 7 • 10 -3 pixel -1 per unit time, and random initial conditions for R taken as uniformly distributed in [0.4,0.5].d) Satellite view of vegetation patterns from a shrubland, 30 miles southwest of Las Cruces, New Mexico (32° 12' 28'' Lat N; 106° 51' 15'' Long W).Images in Fig. 2b, d were converted from the original images to a grayscale raster format, and transformed into a binary matrix (using a threshold filter).

Fig. 3 .
Fig. 3. (a) Cumulative distribution of cluster areas, A, determined for the shrub vegetation shown in Fig. 2a, b.(b) Distributions of cluster areas calculated for the vegetated areas shown in Fig. 2c, d.The distributions obtained for model-simulated vegetation patterns (Fig. 2a,c) closely match those observed in two examples of shrubland ecosystems (Fig. 2b, d).