Vegetation Pattern Formation in Drylands

TL;DR

This chapter reviews mathematical models and observational studies of dryland vegetation pattern formation, emphasizing the integration of remote sensing data and the challenges of modeling across large spatial and temporal scales.

Contribution

It provides a comprehensive overview of conceptual models and observational evidence, highlighting future directions with remote sensing data integration.

Findings

Patterns are observable via satellite imagery.

Mathematical models help explain pattern formation.

Remote sensing data enhances model validation.

Abstract

This is a book chapter, written as a contribution to a new edition of "Dryland Ecohydrology", edited by Paolo D'Odorico, Amilcare Porporato, and Christiane Runyan, (to appear, Springer 2019). It aims to (1) describe some of the background to conceptual mathematical models of spontaneous pattern formation, in the context of dryland vegetation patterns, and (2) review some of the observational studies of the phenomenon. The chapter also highlights challenges and opportunities associated with the development of the models in light of increasing availability of remote sensing data. This includes both satellite imagery of the patterns and elevation data of the topography. The vast scales, in time and space, associated with the key processes further suggest avenues for improved mathematical modeling paradigms.

Tables

Table 1. Table 1: Contrasting the challenges of modeling pattern formation of vegetation patterns with those in the setting of a typical physical system where the theory was developed.

	Wetlands	Drylands
	(Pattern Formation in Fluid Dynamics)	(Pattern Formation in the Environment)
Equations	Navier-Stokes	models exist
Parameters	often excellent specs	some known, order of magnitude level
Time-Scales	$\ll$PhD thesis scale	decades-century
Space-Scales	table-top	meters to landscape-scale
Symmetries	excellent approximation	heterogeneities offer opportunities?
Mechanisms	well developed and validated	generic mechanisms invoked

Equations

$$\frac{\partial u}{\partial t}=ru-(k_0^2+{\nabla }^2{)}^{2}u+f(u),$$

$${\partial }_{T}B={\mathcal{F}}_G{\mathcal{F}}_I-{\mathcal{F}}_M$$

$${\mathcal{F}}_G(\mathbf{X},T)$$

$${\mathcal{F}}_I(\mathbf{X},T)$$

$${\mathcal{F}}_M(\mathbf{X},T)$$

$$w_{\alpha }(\mathbf{S})=\frac{1}{2\pi L_{\alpha }^2}{e}^{-|\mathbf{S}{|}^2/2L_{\alpha }^2}$$

$${\partial }_{t}b=(1-\mu )b+(\Lambda -1)b^2-b^3+\frac{1}{2}(L^2-b){\nabla }^{2}b-\frac{1}{8}b{\nabla }^{4}b$$

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

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

$${\partial }_{T}H$$

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

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

$$\mathcal{G}=gB\frac{W}{W+k},\mathcal{I}=\alpha H\frac{B+qf}{B+q}.$$

$$\nabla \cdot (W\nabla \zeta )\equiv \nabla \zeta \cdot \nabla W+({\nabla }^2\zeta )W,$$

Full text

1 Introduction

Some of the most beautiful examples of spontaneous pattern formation occur in drylands, where the spatial distribution of vegetation may appear as a remarkably self-organized patchiness on a scale spanning kilometers. Figure 1 shows examples of arced bands of vegetation, which alternate rhythmically with bare soil, in different drylands around the globe. Such examples of community-scale self-organization of water-limited ecosystems were first reported when aerial photographs accompanied hydrological survey work in the Horn of Africa, dating to the 1940s. The aerial surveys revealed the extent of these large-scale patterns, something that was not possible from the ground. Modern satellites have since provided a valuable tool for investigating these striking natural patterns.

A sense for “pattern” is well-honed in humans, and yet the term can lead to different interpretations or definitions. In this chapter, we restrict its scope to a mathematical one of possessing spatial symmetries. An idealization of vegetation bands, such as those shown in Figure 1, exhibit translation symmetry that is continuous along the bands. They also exhibit discrete translation symmetry, transverse to the bands; the pattern would look about the same if we translated everything by one stripe. In another type of vegetation pattern, a regular pattern of bare patches appears in an otherwise vegetated region. These types of patterns possess discrete rotation and translation symmetries, where rotating the pattern by a certain amount or translating it by a certain amount brings us to the same pattern. We focus on patterns that emerge in response to aridity stress, where patchiness may represent a strategy for harvesting the limited water from a region of suitable scale to sustain the patch. We would like to understand what sets the scale of the patches and the symmetries associated with their distribution. Even starting from knowledge of the biomass coverage fraction, these are not obviously determined. Instead, to explain the distribution of biomass, it’s necessary to understand the cooperative feedback between plants and nutrients, as well as the mechanisms for spatial redistribution of the limiting resources.

Modeling natural spatial patterns emerged as a challenging test bed for theories of pattern formation after a boom of nonlinear pattern formation research occurred in the 1980s based in laboratory sciences. This chapter highlights examples of both observational studies and mathematical modeling ones. Table 1 summarizes a number of the challenges that arise in modeling of the vegetation pattern formation process, including a lack of fundamental equations and a lack of controlled experiments. In addition, there is an enormous range of potentially relevant temporal and spatial scales. Timescales may vary from those associated with rain events (minutes, hours) to those associated with plant colonization and regeneration (decades, centuries). Likewise spatial scales range between those of individual plants (centimeters, meters), to those of the macroscopic community-scale patterns and landscape (decameters, kilometers). Often models are constructed by considering the most relevant processes, which may act on vastly different spatial and temporal scales, and incorporating them in the model on the time and length scale of the self-organized vegetation patches. These conceptually simple models are amenable to mathematical analysis. They are intrinsically limited, however, because they do not resolve processes on the appropriate scales, nor are they developed through explicitly averaging equations over finer scales to develop the simplified model. Because of this, it can be challenging to constrain parameters, or even to identify how the parameters correspond to measurable quantities. These issues are described through examples in Section 3.5.

The conceptual approach, in contrast to detailed mechanistic modeling, typically exploits our mathematical understanding of generic mechanisms for pattern formation. This perspective is at the heart of proposals for using pattern characteristics as early warning signs of ecosystem collapse, as described in Section 3.2. The expectation is that the pattern formation phenomena may rely on few details of the system beyond which spatial symmetries are present. This, afterall, explains why similar patterns are observed in such disparate physical, chemical and biological systems. In this instance, it can also explain why similar patterns appear in such far-flung regions of the globe, and in different formulations of the mathematical models of vegetation patterns. However, generic arguments often only apply in “a neighborhood” of an instability and if the patterns are of “sufficiently small” amplitude. Moreover, certain mathematical assumptions must be satisfied by the governing equations for these generic arguments to hold. These assumptions may be built into a conceptual mathematical model, but whether the assumptions hold cannot be checked against fundamental equations for the system (fundamental equations would be akin to Navier-Stokes for fluids), because such fundamental equations do not exist. This leads naturally to caveats on model predictions that are often unstated. These limitations on the predictions of models are important to make explicit in the context of predicting the fate of an ecosystem under changes in the degree of environmental stress, or in human land use.

This chapter focuses on observations and conceptual mathematical modeling of vegetation patterns that exhibit a dominant spatial scale and self-organization into a regular pattern. It omits other fascinating threads of dryland pattern formation research where there is no dominant scale, e.g. scale-free patch-size distributions that have been investigated in arid Mediterranean ecosystems (Kéfi et al., 2007) and, more globally, in the context of dryland ecosystem multifunctionality (Berdugo et al., 2017). Nor does it review literature on multiple spatial scale vegetation patterning, such as has been investigated in the context of so-called “fairy circles” (Bonachela et al., 2015; Getzin et al., 2016; Tarnita et al., 2017). The perspective is admittedly skewed towards the pattern formation one that emerged from investigating physical systems. Meron (2012, 2015, 2018) provides further background on the physics approach to pattern formation, specifically in the context of ecological modeling. Borgogno et al. (2009) is an excellent review article on the mathematical modeling frameworks for vegetation patterns. An overview of large scale regular pattern formation in a variety of ecosystems, including mussel beds, wetlands, drylands, ribbon forests and mudflats, is provided by Rietkerk and van de Koppel (2008).

This chapter is organized as follows. Section 2 reviews observational studies of vegetation patterns, starting with the early ground-based and aerial-survey ones from the Horn of Africa, and moving on to modern observational studies that are informed by mathematical models and often incorporate satellite data. Section 2 ends with a summary of the key ecohydrological processes that are implicated in the pattern formation process. Section 3 introduces some of the mathematical theory of pattern formation that informs studies based in model simulations. It describes mechanisms of vegetation pattern formation, emphasizing the way that our understanding of pattern formation in physical and chemical systems has influenced model development. It also presents a sequence of three conceptual mathematical models that frame the pattern formation problem, each of which adds a layer of complexity on its predecessor. Section 3 ends with a discussion of the challenges of constraining parameters in such parsimonious models. Section 4 turns to the modeling challenges that arise because many of the key processes act on vastly different time and space scales. It identifies some possible mathematical directions that might turn this challenge into an opportunity, both for improved modeling, by better constraining parameters, and for speeding up simulations by separating the “fast” from the “slow” scales. Section 4 ends by highlighting two specific directions as promising: modeling stochasticity of rain events explicitly, and exploiting landscape scale topographic heterogeneity to probe models. The chapter closes in Section 5 with our outlook, which emphasizes the importance of including humans in the equations.

2 Observational Studies

We look at some of the observational work on vegetation arcs, and other regular vegetation patterns, and summarize the key ecohydrological processes involved. The observational studies include aerial surveys and satellite image analysis, as well as some ground-based field investigations. Unlike early studies (circa 1950 to 1995), more recent investigations (1995 to present) are informed, in part, by models.

2.1 Early Investigations

Imagery of banded vegetation patterns dates to the first half of the 20th century. In the late 1940s, W.A. Macfadyen, a geologist working on a water supply investigation in what was then British Somaliland (MacFadyen, 1933), identified banded vegetation patterns on aerial photographs and reported them in a letter to Nature (Macfadyen, 1950a). He soon published a more extensive description in the Geographical Journal (Macfadyen, 1950b) where he reproduced aerial imagery of the bands taken in 1945 by the US Army Air Force (U.S.A.A.F.). Although there had been some aerial photographs of Somaliland dating to 1920 (Macfadyen, 1931) and an image of a single band appeared in a 1939 book focusing on soil erosion (Jacks et al., 1939), patterns had not before been reported. Macfadyen’s letter to Nature (Macfadyen, 1950a) generated additional reports of patterns in Sudan and in the Kalahari region.

Figure 2 shows one of the images reproduced in Macfadyen’s paper. The geographic locations of those images, after cropping any horizon lines, are represented by the blue boxes on the map. In addition to working with aerial imagery, Macfadyen surveyed a number of bands. He measured the local slope of the ground where arcs occurred ($0.2$ to $0.3$ % grades), the width of the vegetation bands ($10$ to $76$ m) and the wavelength (or spacing) of the bands ($70$ to $276$ m) (Macfadyen, 1950b). He noted that the arcs were oriented convex upslope in his region of study, thereby providing an indicator, from the aerial images, of the direction of a gentle prevailing grade. (One of his drawings of the arcing relative to the hillslope is reproduced in Figure11a.) In addition to characterizing the bands, he reported on the hydrology, linking it to the arcs. In a region near the bands, Macfadyen described rare, large rainstorms driving a film or sheet of water flowing overland (up to 25 mm deep). Macfadyen thought the relevant storm frequency might be as low as one or two such storms every few years. He speculated that overland water transport of organic matter, including seeds and animal excrement, by such rare storms might account for the origin of the arcs.

Macfadyen’s observations in the Horn of Africa were followed by those of Greenwood, in the same region (Greenwood, 1957). Greenwood focused on Royal Air Force (R.A.F.) aerial images dating to a 1951-52 aerial survey. The locations of images he reproduced in his paper are shown, in red, in Figure 2. The R.A.F. aerial survey images remain available via the Bodleian Library at the University of Oxford (Gowda et al., 2018). Making a site visit in 1955-56, Greenwood examined soil composition and made topographic observations. He linked water availability for vegetation to the inverse of soil clay content, noting that in Sudan, higher clay content soils require more rain to support trees than do less clayey soils. Water availability is currently an important component of modeling efforts, but the role of soil composition is hard to quantity in large scale studies of patterns due to dataset limitations (Bastiaansen et al., 2018).

Greenwood’s topographic observations extended Macfadyen’s report that arcs curve convex up-slope. Greenwood noted that they curve up-slope in shallow valleys, can be virtually straight on slight slopes, and can be convex down-slope when on low ridges. Arc curvature is discussed in Section 4.3 in the context of current mathematical modeling efforts. Greenwood also provided a mechanistic description of how bands might form. He identified vegetative litter mulch as creating a trio of positive feedbacks: it would increase penetration of water into the soil, reduce evaporation from the soil surface, and trap further organic litter from run-on. These are still believed to be important feedbacks in the system: feedback effects associated with increased infiltration and reduced evaporation have been incorporated into contemporary models (Rietkerk et al., 2002; Gilad et al., 2004) and a feedback related to the transport of surface organic matter is also sometimes included (Saco et al., 2007).

From Sudan, J.W. Wright reported grass patterns to Macfadyen in 1950 (Macfadyen, 1950b). In 1959, G.A. Worrall published accounts of those grass patterns (Worrall, 1959) and of tree patterns, in a region with sandier soil, in 1960 (Worrall, 1960b). In the grassy region, stripes were on high ground in between rivers, and arced convex downslope. Here the patterns sit on land with slopes between $0.2\%$ and $0.5\%$ with $8$-$12$ meters of vegetation separated by twice that width of bare ground. Worrall finds a substantial moisture differential between vegetated and bare ground. He observed rain led to moisture penetration of 10-15 cm in bare areas and 30-40 cm under grass, and that vegetation follows moisture, with lusher upslope growth, leading to upslope bands movement. Water transport was impacted by the formation of potholes at the upslope frontal edge of the stripes as well as by cracking of the soil within the stripe. Worrall identified feedbacks between water, wind, topography, soil, and vegetation as all contributing to the patchy biomass distribution associated with vegetation patterns (Worrall, 1960a).

After these initial works establishing the phenomena, observational efforts continued through the 1960s. Boaler and Hodge worked with patterns of grasses, shrubs, and Acacia in the Horn of Africa occurring on slopes of $0.2\%$ to $0.7\%$ (Boaler and Hodge, 1964), and also noted soil potholes impacting water transport at the front and in the body zone of the arcs. Observing a rainstorm, they found infiltration in the vegetated zone was five times deeper than in the bare region between arcs. Boaler and Hodge also speculated about an arc formation mechanism that starts with uniform vegetation under wet (or low grazing) conditions, with bare patches arising as the climate dries (or grazing increases). They imagine that plants adjacent to open spaces might be at a disadvantage, until there is a sufficient bare area for substantial runoff of water, creating a selection mechanism for the size of bare areas.

Hemming also thought about the arcs as being in balance with water availability via runoff and grazing pressure (Hemming, 1965). He thought increasing rainfall or decreasing grazing might widen individual arcs, while decreasing rainfall might reduce the number of arcs. In addition to thinking about how these systems stay in balance, Hemming contributed careful observational work. A map he constructed of an arc near Las Anod in northern Somalia is reproduced in Figure 3. He catalogued vegetation (Hemming, 1966), observed runoff, and measured soil saturation. For example, in one observation of sandy soil (8.64N, 47.58E) with 36 mm of gentle rain, the top 60 cm of soil saturated before runoff occurred. He also made an estimate of 15-30 cm/year migration speed, via upslope colonization, for his surveyed arc. This number is in the range of some contemporary measurements (Bastiaansen et al., 2018; Deblauwe et al., 2012; Gowda et al., 2018).

By 1970 striped patterns had been reported in the Somali Republic, the Republic of Sudan, Niger, Australia, Saudi Arabia and Iraq (Wickens and Collier, 1971). A distinction had been drawn between these stripes, that seem to form from a water-vegetation interaction, and those in Utah in the United States, Jordan, and Mauritania, where striping may be more related to wind-blown material (White, 1971). Theories of stripe formation via water-vegetation interaction were summarized by White (1971). He expanded upon the ideas of Boaler, Hodge, and Hemming by thinking about vegetation arcs across a precipitation gradient. He imagined that stripes might form via bare spots developing in uniform vegetation at the wetter end of the zone, and from colonization of largely bare surfaces at its arid margin. Greig-Smith’s 1979 presidential address to the British Ecological Society summarized the work of Boaler and Hodge (1964), Hemming (1965) and White (1971) on vegetation arcs and theories of their formation, while also providing a much broader context for non-randomness of vegetation distribution across many ecosystems (Greig-Smith, 1979).

By the 1970s, a list of factors that need to be present for vegetation stripes to occur was also well developed. Wickens and Collier (1971) offer a list of factors that is similar to a list offered by Boaler and Hodge (1964) a few years before. They name the following factors: patterns were known to occur in warm arid and semiarid regions on a variety of soils; patterns had been seen at rainfalls of between $100$ and $500$ mm per year, with rain often coming in heavy events; surface slopes ranged from $0.2$% to $2$%; shallow enough slopes that water runs off in sheets, rather than draining via channels; and across all of these locations, bare ground has a low infiltration rate (even with sandy soils, due to surface crust), while vegetated locations accumulate water.

In the 1980s and 1990s, case studies were reported from Mexico, extensive research was ongoing in Australia, and a large-scale study of the Sahel was undertaken in Niger. The work of Cornet and Montaña in the Chihuahuan desert of Mexico (Cornet et al., 1988; Montaña et al., 1990; Montaña, 1992) adds observations of soil moisture, plant migration, and plant community evolution drawn from nearly a decade of study. Their observations at the transition zone where striped vegetation gives way also drove a hypothesis about landscape evolution. They posited that the water drainage system gave way from gully-flow to sheet flow, leading to the development of vegetation stripes (Cornet et al., 1992). In Australia, work by Ludwig and Tongway (1995); Tongway and Ludwig (1990) followed up on the early research of Slatyer (1961) on Acacia aneura grove/intergrove patterns. They identified microtopography (Tongway and Ludwig, 1990), developed monitoring methods (Tongway, 1995), and also made observations of patch size distribution and its sensitivity to rainfall (Ludwig et al., 1999).

In Niger, the HAPEX-Sahel study (Goutorbe et al., 1997) included a “Tiger bush” (banded vegetation) study site. The February 1997 issue of the Journal of Hydrology was a special issue reporting on the study. The issue included detailed rainfall information (Lebel et al., 1997), evaporation measurements (Wallace and Holwill, 1997; Kabat et al., 1997) and carbon dioxide flux (Levy et al., 1997) in the banded vegetation patterns, as well as evapotranspiration measurements at a savannah site nearby (Tuzet et al., 1997). A September 1999 special issue of the journal Catena includes additional case studies from Niger (Chappell et al., 1999; Galle et al., 1999; Valentin et al., 1999), Mexico (Janeau et al., 1999), and Australia (Dunkerley and Brown, 1999; Macdonald et al., 1999) addressing water transport, water harvesting, band migration, and soil cations distribution. In addition, a more recent paper documents an ecohydrological study of a patterned site in the Stockton Plateau of West Texas in the United States (McDonald et al., 2009).

2.2 Recent observational work

Due to advances in model development (see Section 3.4) that occurred in the late 1990s, more recent observational work has often been influenced by modeling efforts. Barbier et al. (2014) have reviewed a number of recent observational results in the context of model predictions. Starting with the work by Barbier et al. (2006), observations are often presented in the context of model predictions, and make use of the multidecadal aerial images to identify changes in vegetation (e.g., in Niger over a forty year interval). The use of multiple images over time was pioneered in work by Wu et al. (2000), where they took advantage of three decades of aerial photographs in Niger to examine fragmentation. Barbier et al. (2006) found that extended drought in their study region in Niger was accompanied by a shift from homogeneous vegetation to spotted patterns, and that the shift was lessened in areas protected from wood cutting and grazing. They connected their results to modeling efforts by noting that the observed spatial wavelength of the patterning appeared to be endogenous. They saw this as offering support for models that rely on a pattern formation mechanism to generate patterns (see Section 3.3.)

Recent observational work has often made use of remotely sensed data, and sometimes of contemporary statistical learning methods. Some new pattern locations were identified by Deblauwe et al. (2008) after studying correlates of patterning using a MaxEnt machine learning algorithm developed for identifying factors of ecological niche segregation (Elith et al., 2011). Figure 4 shows some of the pattern locations in North America, northern Africa, and Australia, where they are in black over a map based on the Köppen-Geiger climate classification system, indicating that patterns are mainly observed in warm arid (BWh) or warm semiarid (BSh) regions.

There has also been some work oriented towards directly constraining model parameters (Barbier et al., 2008), however, many recent studies focus on measuring quantities that are common outputs of a model, such as pattern morphology or pattern wavelength, as a function of precipitation or other environmental factors (Deblauwe et al., 2012, 2011; Penny et al., 2013). Pattern morphologies include homogeneous vegetation cover, vegetation with bare patches (gaps), labyrinthine patterns, isolated patches or spots of vegetation, and bare ground. In models, different morphologies are associated with different levels of mean annual precipitation, and if precipitation levels were to change with time, transitions between pattern states is often predicted. See Section 3.2 for more about pattern morphology and critical transitions.

In a study of Sudan by Deblauwe et al. (2011), they analyze satellite imagery from 1967, 1988, and 2001-2002 using Fourier-based methods to identify pattern morphology as a function of the rainfall gradient across the study region (approx. 350mm/yr to 500 mm/yr). Figure 5 in Section 3.2 is reproduced from their work, and indicates model predictions of pattern morphology in an aridity-slope parameter space. Their observations qualitatively matched the predictions of how patterns would be distributed along a rainfall gradient. In addition, the time interval in their study included decades of low precipitation. The change in precipitation coincided with regions that transition from gapped patterns to labyrinths as well as regions that transition from labyrinths to spotted patterns. This type of morphological change with reduced precipitation aligns with model predictions. In addition to comparing pattern morphology, features such as the wavelength and orientation of bands relative to the local slope (Bastiaansen et al., 2018; Penny et al., 2013), the width of bands (Gowda et al., 2018), the local biomass (Bastiaansen et al., 2018), and arc curvature (Gandhi et al., 2018a), have been identified using satellite or aerial imagery and compared to model predictions. Other work has focused on vegetation patch size distribution (Kéfi et al., 2007) as an indicator of desertification (Berdugo et al., 2017; Génin et al., 2018; Moreno-de las Heras et al., 2011).

2.3 Summary of Key Ecohydrological Processes

Following the basic pattern formation concepts, often attributed to Alan Turing, vegetation patterning is thought to be the result of an interplay between short-range cooperative and long-range competitive plant interactions (Borgogno et al., 2009; Lejeune et al., 2004; Rietkerk and van de Koppel, 2008). Here we summarize some of these feedbacks, many of which have been incorporated in mathematical models of vegetation pattern formation.

Positive feedbacks between plants and water include the ability of plants to concentrate resources (e.g., runoff interception, shadow, nutrient-rich litter) providing more favorable conditions for plants to establish. In particular, vegetation cover intercepts runoff thus favouring infiltration and increasing local soil water availability (Meron, 2015; Paschalis et al., 2016; Saco et al., 2007; Valentin et al., 1999; Ursino, 2005), and protects against soil erosion (Ludwig et al., 2005; Saco et al., 2007; Valentin et al., 1999). In addition to this runoff-runon mechanism, infiltration rates under vegetated patches are enhanced as a consequence of the increased soil macroporosity related to biological activity (e.g., termites and earthworms) and vegetation roots (Saco et al., 2007). In contrast, bare areas are maintained by formation of soil crusts that reduce water infiltration, preventing vegetation from establishing and growing (Assouline et al., 2015; Borgogno et al., 2009; Meron, 2015). In addition to this positive infiltration feedback, vegetation cover may decrease temperature fluctuations, exposure to solar radiation, and wind desiccation, thus reducing water losses and further favouring plant growth. Another positive feedback mechanism is the protection from fires and grazing (Borgogno et al., 2009). The concentration of resources (i.e., water and nutrients) in vegetated and fertile soil patches – contrasted to unfertile bare areas – has led to the concept of “fertility islands” (Ridolfi et al., 2008; Schlesinger et al., 1990). Species that are capable of modifying the abiotic environment by resource redistribution, and thus facilitating the growth of other species, are known as “ecosystem engineers” (Gilad et al., 2004; Jones et al., 1994). Conversely, competition for resources among plants is exerted thorough the root system which extends beyond the canopy area and extracts water and nutrients from the bare patches (Meron, 2016; Rietkerk and van de Koppel, 2008). This negative feedback takes place over a range that is larger than that of facilitative interactions.

In addition to these positive and negative feedbacks, a number of exogenous climatological and geomorphological factors can further affect vegetation patterning. In fact, spatial heterogeneities in soil properties and topographic features are likely to alter water runoff, infiltration, solar radiation and erosion processes (Paschalis et al., 2016), ultimately changing the spatial redistribution of resources and thus the pattern formation mechanism. A number of studies empirically analyzed the relation between pattern geometry and “external” environmental conditions, such as rainfall, temperature, ground slope, soil type and other topographic variables (Augustine, 2003; Deblauwe et al., 2008; Penny et al., 2013). Analysis of the location of vegetation patterns (i.e., their “ecological niche”) provided evidence for a correlation between vegetation patterns and environmental constraints (namely, humidity index, mean temperature and its seasonality, rainfall amount and seasonality, soil wilting point and ground slope) at global scale (Deblauwe et al., 2008). Pattern morphology also seems to be highly related to large scale properties of the landscape such as slope orientation, hillslope gradient, and soil type, as well as smaller scale features (e.g., ridges, streams, roads, and soil peculiarities), which can produce deviations from the global trends (Penny et al., 2013).

3 Nonlinear Dynamics Perspective: Critical Transitions and Pattern Formation

The classic study by Pearson (1993) of pattern formation in the Gray–Scott model (Gray and Scott, 1983) demonstrates the rich and often unexpected pattern forming behavior by even the simplest of reaction-diffusion systems. It is therefore not surprising that models of this type can reproduce observed spatiotemporal patterns. Section 3.1 discusses the development of the mathematical framework within nonlinear dynamics that connects the emergence of patterns across a broad range of contexts. The key idea is that as a parameter changes a uniform (vegetated) state loses stability and a patterned state emerges. The usefulness of this approach lies in the fact that the predictions are often insensitive to details of a given model, at least for parameters near where the patterns form. This is described further in the context of early warning sign proposals in Section 3.2. Section 3.3 provides more details of the mathematical mechanisms behind “pattern forming instabilities” in vegetation models, with a chronology of their development in Section 3.4. We also include some details of three particularly influential phenomenological models in Section 3.5, and conclude in Section 3.6 by illustrating the challenges faced when trying to constrain parameters in these models.

3.1 Pattern Formation

The study of spontaneous spatial pattern formation took off as a sub-field of nonlinear dynamics research in the 1980s, although it has a much longer history within fluid dynamics research. See Cross and Hohenberg (1993) for a general review from physics, and Crawford and Knobloch (1991) for a review in the fluid dynamics setting. The nonlinear dynamics perspective emphasizes “generic mechanisms”. For example, those of bifurcation theory, which is a mathematical framework that helps elucidate the origin of bistability, hysteresis, and spontaneous oscillations in nonlinear systems. Studies in nonlinear dynamics helped to identify fundamental routes to chaos, and also the transitions to patterns in spatially extended nonlinear systems. With a focus on generic mechanisms, connections are made between strikingly similar emergent behavior in chemical, biological and physical systems.

At the heart of pattern formation research is a basic notion of “spontaneous symmetry-breaking”, which is associated with the onset of an instability. Consequently, genericity arguments apply for special parameter sets that are close to those where an instability sets in. These arguments do not require much detailed knowledge about the system; just knowledge about the type of instability (e.g. whether it is oscillatory or steady/monotonic), and of the associated spatial symmetries of the system, and then one can predict what kinds of patterns are possible (Golubitsky et al., 1988; Hoyle, 2006). While the stability of these patterns depends on details of the model nonlinearities, the characteristic spatial and temporal scales are often determined, approximately, by the quantitative linear computations that are required to determine the modes of instability. However, the results are technically only valid for “small amplitude patterns” near the onset of instability in parameter space. Some care must be taken when extending the results to the states actually observed in the environment, which go between well-developed community-scale vegetation patches and bare soil with zero biomass, and hence are not “small amplitude” patterns.

Much of the research on pattern formation in the 1980’s and 90’s involved a healthy back and forth between laboratory experiments and theory. An interesting counter point to this, which is relevant to vegetation pattern modeling, involves the so-called “Turing mechanism”. This mechanism is attributed to a mathematical argument of Alan Turing, which appeared in his 1952 paper “The Chemical Basis of Morphogenesis” (Turing, 1952). In it, he shows that diffusing and reacting chemical species, which he calls “morphogens”, may undergo an instability that leads to a spatial pattern, i.e., the diffusing chemicals do not evolve to have spatially uniform concentrations. An essential ingredient for this pattern formation process is that the two species diffuse at drastically different rates: an activator chemical diffuses slowly, acting only at short-range, while an inhibitor chemical diffuses rapidly, and acts at long-range. The formation of a pattern with an emergent characteristic lengthscale is, on the surface, quite surprising because it relies on diffusion, which we generally associate with destroying any sharp boundaries, such as those that appear in the patterns. The emergent pattern lengthscale is tied to the competition between two drastically different diffusive lengthscales and this scale cannot be simply identified with some obvious physical scale of the problem. Interestingly, Turing’s theoretical insight, which suggested that “dappled patterns” might emerge when the reacting chemical species diffuse at different rates, was not actually realized in a laboratory for decades following his proposal. The difficulty was to have the chemical species diffuse at different rates; this was finally achieved in the laboratory in the early 1990’s (Castets et al., 1990; Ouyang and Swinney, 1991), when a chemical reactor was made of a gel which only inhibited motion of the (larger) “activator” molecules.

The Swift–Hohenberg equation (Swift and Hohenberg, 1977) is a particularly simple model that demonstrates how patterns form via a symmetry–breaking instability of a spatially uniform state. The model is a partial differential equation (PDE) for a scalar field $u(x,y,t)$,

[TABLE]

where $\nabla^{2}\equiv\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$ is the Laplacian, and $f(u)$ is a suitable nonlinear function of $u$. The PDE was originally introduced by Swift and Hohenberg, with the nonlinearity $f(u)=-u^{3}$, as a simplified model of Rayleigh-Bérnard convection. In that setting, the scalar field $u$ describes the mid-plane vertical velocity of a thin layer of fluid, heated from below, and the driving from an imposed temperature gradient is captured by the control parameter $r$. The linear stability of $u=0$ to spatially periodic perturbations $e^{ikx}$ is determined by examining their growth rate $\sigma_{k}=r-(k_{0}^{2}-k^{2})^{2}$, which is a simple quartic polynomial in $k$ with a maximum at $k_{0}$. It follows that for $r<0$, the heat conduction state $u=0$ is linearly stable, and for $r>0$, the conduction state is unstable and gives way to a convection pattern. In this case, the characteristic wavelength of the pattern is $\sim 2\pi/k_{0}$ near onset, and in the case of Rayleigh-Bénard convection this scale is set by the thickness of the fluid layer.

As illustrated in Section 3.5, many models of vegetation pattern formation have features that are qualitatively similar to the Swift–Hohenberg equation (including in some cases the quartic growth functions $\sigma_{k}$). These similarities are only expected to be valid in a neighborhood of parameter space near the pattern forming instability, e.g. near $r=0$ for $|u|\ll 1$. Validating models of vegetation pattern formation is challenging when they are intended to go beyond a Swift-Hohenberg model. Not least of these is that, unlike Rayleigh-Bénard convection, controlled laboratory experiments do not play a role, due to the long timescales and large spatial scales of community-scale ecosystems. Nonetheless, the power of a generic mechanism, which typically operates for parameters near some stability threshold, resides in its being largely insensitive to the details. This approach is described in section 3.2, where thresholds are often associated with these special parameter sets near critical transitions.

3.2 Critical Transitions and the Dream of Early Warning Signs

One framing of the potential significance of dryland vegetation pattern formation is in the setting of early warning signs of desertification and ecosystem collapse. For instance, in many mathematical models, as the precipitation parameter decreases, there is a well-defined sequence of transitions in the pattern morphology prior to complete collapse to a barren desert state (Deblauwe et al., 2011; Gowda et al., 2016; Meron et al., 2004; Rietkerk et al., 2004; von Hardenberg et al., 2001). This is illustrated in Figure 5 in the particular case of the “flat terrain patterns”, where, with increasing aridity uniform vegetation gives way to a matrix of bare patches called “gaps”. With further increases in aridity stress, the gaps merge to a labyrinthine pattern that eventually disintegrates to a new matrix, now of vegetation “spots”, before its ultimate collapse to bare soil. The pattern sequence can exhibit significant hysteresis, due to the nonlinear feedbacks that favor bistability of states. Thus it is not easy to reverse a transition from one type of state to another, and in this case, the spot pattern, as precursor to complete collapse, would serve as an “early warning sign” of desertification. A perspective piece by Rietkerk et al. (2004) identifies other ecosystems, including peatlands and savannas that may show the same pattern sequence when subjected to increasing environmental stress.

The similarity between gap/spot patterns and patterns of hexagonal Rayleigh-Bénard convection cells suggests that the same generic mathematical pattern formation framework might be applied to both. In this analogy, convection roll patterns would be supplanted by the labyrinthine patterns. (See Figure 6.) Exploiting this analogy, Gowda et al. (2014) determined the minimal mathematical ingredients for the pattern sequence gaps$\to$labyrinths$\to$spots. It required a change of sign of a quadratic coefficient in a bifurcation problem, as the precipitation parameter decreased. The ideas were firmly grounded in a small amplitude, weakly nonlinear theory of pattern formation for which retaining only the leading quadratic and cubic terms in a Taylor expansion can capture nonlinear essentials. As such, though, it can make no claims about large amplitude, “fully nonlinear” patterns. In a follow-up paper to this framing, Gowda et al. (2016) tested it in the setting of the Rietkerk model, described in Section 3.5, by performing both the weakly nonlinear analysis, and extensive numerical simulations outside of the small amplitude regime. In the setting of this particular model, they were able to test the validity of their theoretical quadratic coefficient “proxy” for the pattern sequence, and trace the ecological origin of the sign change to key feedbacks incorporated in the model. Thus, at least for the Rietkerk model, they showed that the gaps$\to$labyrinths$\to$spots sequence is robust, i.e. insensitive to changes in other (fixed) model parameters.

While a number of different mathematical models show a gaps$\to$labyrinths$\to$spots sequence as the precipitation parameter is decreased (holding other parameters fixed), it is a challenge to observe this sequence of transitions, at a given location over time, in nature. Perhaps the best observational evidence for this sequence is by Deblauwe et al. (2011). In their study they classified vegetation patterns in a 22,255 km2 region that has a strong precipitation gradient, north to south, as it transits out of the Sahel in Sudan. This pattern sequence relies strongly on the terrain being sufficiently flat, as indicated in Figure 5. In many regions, including the patterned region of Sudan, the flat terrain patterns turn to banded ones transverse to the grade (Deblauwe et al., 2011) when there is even a very modest slope. Thus the possibility of a spot pattern as a warning sign of ecosystem vulnerability may vanish with this change in terrain. For possible early warning signs in the context of banded patterns, as precipitation decreases, simulations and analysis of mathematical models have suggested that band spacing may increase (Siteur et al., 2014b), the width of vegetation within the bands may decrease (Yizhaq et al., 2005), and that the speed of fronts, recovering from disruption, might slow (Zelnik and Meron, 2018).

The early warning signs framework relies on identifying the appropriate generic mechanism associated with the critical transitions. Often this exists as an asymptotic (long time) state rather than the transient one associated with, say, a $\sim$decadal timescale. Simulations of mathematical models investigate the response to stress on the system, usually by varying the precipitation parameter with everything else fixed, and trusting it is not too sensitive to heterogeneity in the environment. One expects, though, that a gradual change in aridity will change the composition of the plant community. Interestingly, when one looks at long term change of individual vegetation bands in the Horn of Africa, much is strikingly unchanged, and often the changes that are noted do not align with a clean story from models. For example, a recent study by Gowda et al. (2018) measured the change of thousands of individual vegetation bands, between their 1950’s aerial photographs and their modern satellite images. They found that the bands were surprisingly unchanged in regions of Somalia where there were no signs of human impacts. Moreover, in the human impacted areas, the most striking change to the bands that remain intact is that they thickened. Some of this band widening is evident in Figure 7, which is a false color overlay of a 1950’s aerial image, geo-referenced with a modern satellite one. The upslope migration of the bands, at a rate predicated by Hemming (1965), is also evident in this overlay.

We note that there is an alternative approach to early warning signs, not based in pattern, but instead based on a statistical physics analysis of the scaling distribution of vegetation patch size, and its deviation from a power-law distribution (Kéfi et al., 2007). Spatial heterogeneity has been noted in this case to impact the performance of early warning indicators (Génin et al., 2018). For a more comprehensive review of early warning signs in the context of spatial patterns in ecological systems, see (Kefi et al., 2014).

3.3 Mathematical mechanisms underlying pattern formation in vegetation models

The work of Gierer and Meinhardt (1972) provides a conceptual framework for interpreting Turing’s mathematical results on diffusion driven instability in terms of short-range facilitation and long-range inhibition. Many models of spatial pattern formation in ecology, and more broadly, now commonly rely on this as a general principle for generating pattern forming instabilities in which a uniform state loses stability to perturbations of some finite wavelength (Levin and Segel, 1985). This interpretation leads to a generalization of the Turing mechanism described in Section 3.1. Some heuristic insight can be gained by considering a small local perturbation of the slowly diffusing activator and rapidly diffusing inhibitor. The local increase in the activator leads to local increase in both the activator and inhibitor. The inhibitor will spread more quickly and depress the concentration of both the activator and inhibitor surrounding the local disturbance. If the diffusion rates are sufficiently disparate, the inhibitor can leave the local perturbation without suppressing it and the perturbation can grow because an excess of the activator remains. This type of instability is often invoked in vegetation models under an assumption that water “diffuses” much more quickly than biomass on flat terrain. Biomass, due to cooperative affects of increased infiltration, acts as the “activator”. The “lack-of-water” (as a positive quantity) would then be considered an inhibitor (Meron, 2018).111For example, a transformation from the water field $W$ to a “lack-of-water” field $W_{L}=W_{0}-W$ in the Klausmeier model (Klausmeier, 1999), for some constant $W_{0}$ such that $W_{L}>0$, can lead to an activator-inhibitor system for $(B,W_{L})$, where $B$ is the biomass field.

A generalization of the Turing instability can appear in kernel-based models where the temporal dynamics at a location are governed by integral terms that capture influence from the neighborhood surrounding a particular location. For example, this is often used to capture plant growth for cases when root systems harvest water over an extended area surrounding a given plant. This provides long range inhibition while short range facilitation can be obtained through moisture trapping over smaller areas. From the perspective of symmetry-breaking, the resulting bifurcations are indistinguishable from Turing bifurcations. The dynamics of both types of systems can therefore be described by the same class of reduced models in appropriate limiting cases.

A variation of this general framework for pattern formation is further exemplified by models of banded vegetation patterns on hillslopes. In this case the water undergoes fast advection instead (i.e. downslope flow) of fast diffusion, leading to a so-called Turing-Hopf instability (Rovinsky and Menzinger, 1992). An immediate consequence of spatial asymmetry, between uphill and downhill, is that patterns are expected to travel and might then be captured in models as traveling wave patterns. Vegetation models such as Klausmeier’s (Klausmeier, 1999) predict the travel direction to be uphill when advection is introduced, consistent with some of the observational evidence at a number of locations (Deblauwe et al., 2012). In a Turing system where a pattern is generated by diffusion, an asymmetry can also be introduced, e.g. by adding an advective term. Asymmetry can also be incorporated into kernel-based non-local models, typically by introducing asymmetry into the integration kernel. Both of these cases generate similar behavior of traveling wave patterns that can model banded vegetation with an upslope colonization.

The perspective discussed in this section focuses on the “formation” of the pattern at a particular point in parameter space through an instability of a uniform state. There are a number of possible scenarios by which an ecosystem can evolve in time to a particular patterned state. Environmental changes could bring a once stable uniform state to an unstable region of parameter space and the pattern could grow from perturbations or fluctuations. For instance, a decrease in water availability could lead to the emergence of patterns from small-scale variation in an initially uniform vegetation state. Another possibility is that the vegetation patterns colonize a region of bare soil via an invasion front. See Samuelson et al. (2018) for a bare bones model that illustrates the latter scenario. Studies of front dynamics in vegetation models have suggested restrictions on what scenarios are possible (Zelnik et al., 2017). The prediction in this case is that the colonization of bare soil by vegetation patterns may not be gradually reversed by a receding front, and the patterns are instead lost all at once through a uniform collapse to bare soil. Other studies have suggested that the pattern selected as a result of colonization may be more consistent with empirical data than the patterns formed from small perturbations of a uniform state (Sherratt, 2015).

3.4 Development of Mathematical Models for Vegetation Pattern Formation

The first models of vegetation pattern formation seem to appear around the mid 1990’s (Mauchamp et al., 1994; Thiery et al., 1995), and are discrete-time and discrete-space cellular automata models. The state of each location on a spatial grid, at each time, is determined, via some set of defined rules, by the state of itself and its neighbors at the previous time. The state often captures characteristics about both biomass and water at a given point in space and time. The models can be extremely mechanistic in the way they capture the various processes involved, leading to a large number of parameters that must be determined from data. This modeling approach tends, therefore, to be very closely connected to observation and often includes comparison to field studies.

The first efforts to incorporate ideas from pattern formation come from Lefever and Lejeune (1997). On a historical note, Lefever made important contributions to the study of pattern formation in nonlinear optics (Lugiato and Lefever, 1987) and in chemical systems. Together with Prigogine, he formulated the Brusselator model (Prigogine and Lefever, 1968). The Lefever and Lejeune model is a conceptual model that captures dynamics of a continuous biomass field using a integro-differential equation accounting for plant-to-plant interactions mediated by the environment. Other work from this branch of modeling has focused on spatially localized patches of vegetation as examples of dissipative structures (Lejeune et al., 2002).

Explicit modeling of environmental dynamics by coupling a continuous biomass field to a continuous water field in a pair of reaction-advection-diffusion equations is first introduced by Klausmeier (1999) and later analyzed in detail by Sherratt (2005); Sherratt and Lord (2007). Klausmeier’s original model is formulated for banded vegetation patterns on a uniform hillslope. Recent work has generalized the model so that it can capture flat-terrain patterns as well (van der Stelt et al., 2013). Using one model to generate both flat-terrain and uniform hillslope patterns, a goal first introduced by von Hardenberg et al. (2001), may sometimes be appropriate. An interesting case where they appear in proximity is in the region of Sudan studied in (Deblauwe et al., 2011). It is not the case that the two types of patterns always occur in association, e.g. flat terrain patterns have not been described in the Horn of Africa, where banded patterns are prevalent (Gowda et al., 2018). It is an open question whether the mechanisms are fundamentally different between formation of flat-terrain and of banded hillslope patterns.

Contemporary reaction-advection-diffusion models are moving towards capturing more mechanistic details of various processes thought to play a key role in pattern formation. For example, the model due to HilleRisLambers et al. (2001); Rietkerk et al. (2002), which we describe below as ”the Rietkerk model,” subdivides the water into a surface water and soil moisture to better capture the biomass feedback on the infiltration of water from the surface into the soil. Gilad et al. (2004) also incorporate this three-field (two water and one biomass) modeling approach. Nonlocal seed dispersal has been explored (Bennett and Sherratt, 2018; Pueyo et al., 2008; Thompson et al., 2008). Some focus on the effects of seasonality or stochasticity of rainfall (Kletter et al., 2009; Konings et al., 2011; Siteur et al., 2014a). Others couple the hydrological processes to landform evolution via sediment transport dynamics (Saco et al., 2007). There is also now a push to model the interaction of multiple species of vegetation within these patterned ecosystems (Gilad et al., 2007a). The influence of grazers (van de Koppel et al., 2002; Siero, 2018) and interactions with termites (Pringle and Tarnita, 2017) have been considered as well. Even hybrid models that capture plants with an agent-based approach and couple them via a continuous water field have been proposed as a way to incorporate individual plant characteristics (Vincenot et al., 2016).

3.5 Examples

We now describe three examples from the literature that serve as foundations for many contemporary reaction-advection-diffusion models of vegetation patterns. We assume a two-dimensional spatial domain throughout this discussion, though many model studies restrict to one spatial dimension for simplicity. The Lefever–Lejeune model (Lefever and Lejeune, 1997), perhaps the first to take a pattern formation perspective, captures self-interactions of biomass through nonlocal integral kernels that exhibit short-range facilitation and long-range inhibition. The Klausmeier (1999) model explicitly models water dynamics coupled to biomass evolution on a uniform hillslope in a pair of equations reminiscent of the Gray–Scott chemical reaction model (Gray and Scott, 1983). Pattern formation in this case occurs as a result of a fast downhill advection of water and slow diffusion of biomass, and does not produce patterns in absence of a hillslope. A model by Rietkerk and collaborators (HilleRisLambers et al., 2001) subdivides water into a fast-moving surface flow and a slowly-diffusing soil moisture that are coupled to a slowly-diffusing biomass field in order to better capture biomass feedback on infiltration of water from the surface into the soil. The increased complexity of the Rietkerk model relative to the previous two examples illustrates the trade-offs between fidelity and parsimony in models of this type. See Zelnik et al. (2013) for a more detailed mathematical comparison of versions of these models to a version of a model by Gilad et al. (2007b).

Lefever–Lejeune Model

The Lefever–Lejeune model (Lefever and Lejeune, 1997) has a very simple structure, with the dynamics of biomass field $B(\mathbf{X},T)$ depending on a growth term $\mathcal{F}_{G}$ multiplied by an interaction term $\mathcal{F}_{I}$ minus a mortality term $\mathcal{F}_{M}$.

[TABLE]

Each of these terms represents an integral that captures nonlocal interaction between the plants:

[TABLE]

Here the integrals are over a two-dimensional spatial domain, i.e. both $\mathbf{X}$ and $\mathbf{S}$ represent two-dimensional vectors. The “replication constant” $\lambda$ captures the growth rate of the plants, the “cooperativity constant” $\Omega$ characterizes plant facilitation, the “biomass packing limit” $K$ plays the role of a carrying capacity, and the “mortality constant” $\eta$ characterizes biomass loss. The integration kernels $w_{\alpha}$ with $\alpha=G,I,M$ are typically assumed Gaussian, and on a two-dimensional spatial domain, take the form

[TABLE]

The lengthscales responsible for the pattern-forming instability are contained in a growth kernel $w_{G}$ in the growth term $\mathcal{F}_{G}$, with a characteristic “dissemination lengthscale” $L_{G}$ for seed dispersal, and a competition kernel $w_{I}$ in $\mathcal{F}_{I}$, with a characteristic “inhibition length” $L_{I}$. There is also a “toxicity lengthscale” $L_{M}$ associated with the mortality kernel $w_{M}$, but it is not expected to play a role in the mechanism for pattern formation, and is often taken to be zero (i.e. $\mathcal{F}_{M}=\eta B$).

With $L_{M}=0$ and in the limit of low biomass densities and small biomass gradients222The key role of these assumptions is to ensure that substituting a Taylor expansion of the biomass into the integrals provides good approximations., the equation can be approximated by the fourth-order PDE

[TABLE]

where time has been rescaled by birth rate $t=\lambda T$, space has been rescaled by inhibition length $\mathbf{x}=\mathbf{X}/L_{I}$ and biomass has been rescaled by carrying capacity $b=B/K$. The dimensionless parameter $\Lambda=K\Omega$ measures availability of resources and $\mu=\lambda/\eta$ is a birth-to-death ratio. This reduced equation shares many qualitative characteristics with the Swift–Hohenberg equation (1); indeed a variant of the model (2) reduces exactly to an equation of Swift–Hohenberg type in an appropriate limit (Lefever et al., 2009).

The Lefever–Lejeune model does not explicitly track the dynamic evolution of water. It does, through the choice of parameters and functional form of the kernels, implicitly make assumptions about the environment. Because the model does not attempt to capture details of hydrology or environment, the assumptions about these are minimal. On the other hand, it can be difficult to use this model to study the influence of the environment since it is often unclear how any given process is encoded. For example, rainfall is not explicitly included as a parameter, but $\mu$ indirectly captures the influence of aridity under the assumption that plants lose biomass at a higher rate in more arid environments. However one might expect that $\Lambda$, which characterizes resource availability, is effected by aridity as well. As another example, the effect of a uniform hillslope can be captured in the model by introducing anisotropy in a weighting function that breaks symmetry along the uphill-downhill direction.

Klausmeier Model

The Klausmeier model (Klausmeier, 1999) couples the dynamics of the biomass field, $B$, to the dynamics of a water field, $W$, on a uniform hillslope where uphill is in the $+X$ direction:

[TABLE]

Water is input into the system with precipitation rate $P$, assumed to be the constant annual mean. The advection rate, $V$, depends on the slope of terrain and represents the overland water flow; $V>0$ means the hill slopes upwards in the $+X$-direction and the water is transported in the opposite downslope direction. The soil evaporation rate is given by $L$. The transpiration, $RWB^{2}$, drives biomass growth with a water-use efficiency parameter $J$. The quadratic dependence on biomass in the growth captures, heuristically, the positive biomass feedback on infiltration of water into the soil. The mortality rate $M$ captures biomass loss due to both natural death, but also due to grazing and to human pressure. Seed dispersal is modeled through linear diffusion with rate $D$.

Figure 8 shows a typical bifurcation diagram associated with a pattern-forming vegetation model such as (5), with shading to indicate the range of mean annual precipitation parameter $P$ where patterns exist. Two important examples of patterned states are shown; one bifurcates at the Turing-Hopf point that determines the right boundary, and the other is a single “pulse” or “oasis” state that marks the left boundary of existence. Also, shown in this figure are examples of growth rate functions $\sigma_{k}$ associated with small periodic perturbations $e^{ikx}$ of uniform solutions.

One advantage of Klausmeier’s model over the Lefever–Lejeune equation is that each term in the $W$ equation can be associated with a specific ecohydrological process. For example, the topography influences water transport and so it is clear which term is capturing the hillslope’s effect. However, because the $W$ field captures a conglomeration of hydrological processes, across timescales, and both on the surface and in the soil, it is not consistent with the observed hydrology in many ways. For example, the model predicts that the water field peaks in bare soil regions as can be seen in Figure 8(c), yet soil moisture is observed to be higher in the vegetated zones after rain events, as described in Section 2.1.

The reaction terms of the Klausmeier model are a reparametrization of the Gray–Scott model, a classic activator-substrate system. The generalized Klausmeier/Gray–Scott model (van der Stelt et al., 2013), which additionally includes a water diffusion term, can be transformed into the standard form of Gray–Scott model studied by  Pearson (1993) in the case that the water transport is purely diffusive.333However, the parameters used in the context of vegetation patterns correspond to unphysical parameter choices in the original chemical kinetics setting. This purely diffusive version of the model has also been shown to reduce to an equation of Swift–Hohenberg type, i.e. Equation (1), for the pattern amplitude in an appropriate limit (Gandhi et al., 2018b).

Rietkerk Model

The model proposed by Rietkerk et al. (2002); HilleRisLambers et al. (2001) takes an additional step to remedy the poor hydrological predictions of the Klausmeier model by subdividing the water into a surface water field and soil moisture field. An ordinary differential equation precursor to this model, aimed at studying the effects of grazing in water-limited and nutrient-limited ecosystems (van de Koppel et al., 1997), was published in the same year as Lefever and Lejeune’s model. The Rietkerk model for vegetation patterns is given by

[TABLE]

where water is input to the surface by precipitation with constant rate $P$ and quickly diffuses along the surface with rate $D_{h}$. After infiltrating into the soil, water diffuses with rate $D_{w}$, is lost by evaporation or leakage deeper into the soil at rate $r$, and is lost by plant uptake characterized by the nonlinear function $\mathcal{G}$. Seed dispersal is modeled by linear diffusion with rate $D_{b}$, plant mortality from grazing or otherwise occurs at a rate $m$ and growth is given by $c\mathcal{G}$. The nonlinearities of the model arise through the growth and infiltration functions, given by

[TABLE]

The parameter $g$ sets the maximum growth rate and the parameter $k$ sets the soil moisture required for half that growth rate to occur. The maximum infiltration rate is set by $\alpha$ and the dimensionless feedback parameter $f\in[0,1]$ controls the ratio of the bare soil to vegetated soil infiltration, e.g. it represents the infiltration reduction introduced by a soil crust. Specifically, bare soil infiltrates at reduced rate $\alpha f$, while when the biomass $B\gg q$, the infiltration approaches its maximum rate $\alpha$.

Separating surface water $H$ from soil moisture $W$ allows for more detailed modeling of the biomass feedback on infiltration of water from the surface into the soil, a process long thought to play an important role in pattern formation of these ecosystems, as described in Section 2.1. While this improves predictions regarding soil moisture distribution, the surface water is still problematic to interpret. In order to account for the assumption that water is input at a constant rate year-round, the infiltration rate is reduced by an extra factor of the fraction of time it rains during the year. Even with this caveat, the Rietkerk model provides a sound framework for exploring predictions from extensions that capture additional processes in more detail. However, more detailed modeling comes at the cost of increased complexity. The scaled version of the Rietkerk model has six dimensionless parameters, as opposed to three for the dimensionless versions of the Klausmeier and Lefever–Lejeune models.

3.6 Challenge: Constraining Parameters in Parsimonious Models

Models appear along a spectrum in terms of level of details about the underlying processes, but can be broadly classified as either “conceptual” or “mechanistic.” Mechanistic models attempt to capture the actual processes with high fidelity, but this often leads to increases in mathematical complexity and number of parameters. Despite the parsimony associated with conceptual models, they are able to reproduce the observed patterns and are often amenable to mathematical analysis. However, the connection between the model parameters and the actual processes being modeled can be tenuous. Constraining model parameters based on data can therefore be challenging and the values obtained can be difficult to interpret. Moreover, the parameter values associated with a given process can vary over orders of magnitude between models. The typical approach to estimating parameters in the more conceptual models involves a combination of constraints based directly on observation or measurement, and inferred constraints that rely heavily on model assumptions. The inferred constraints can be based on characteristics of the bare soil state, uniform vegetation state or patterned state.

We now illustrate the process for choosing parameter values within the most parsimonious of models within the reaction-advection-diffusion framework, following the approach of Klausmeier (1999). A similar approach based on observed pattern characteristics and linear predictions has been used to constrain two key dimensionless quantities in the Lefever–Lejeune model (Lejeune et al., 1999). More detailed studies have also attempted to fit the plant interaction kernel using data from measurements of the root system (Barbier et al., 2008).

In Klausmeier (1999), the constant rainfall parameter $A$ was set to 300 mm/year, a typical mean annual precipitation for semiarid regions. Values for plant water use $J=0.002$ kg Dry Mass/kg H2O for trees ($J=0.003$ kg Dry Mass/kg H2O for grasses) and mortality $M=0.18$ year*-1* for trees ($M=1.8$ year*-1* for grasses) were taken from estimates in a previous modeling study (Mauchamp et al., 1994), which were based on data about maximum water use efficiency and carbon costs required to maintain leaves. The evaporation rate $L=4$ yr*-1* is set by typical values of water from observation and predictions about the bare soil state. The transpiration rate $R=1.5$ mm H2O/yr/(kg dry mass)2 for trees ($R=100$ mm H2O/yr/(kg Dry Mass)2 for grasses) is set by typical biomass values from observation and predictions about the uniform vegetation state. This leaves the two transport parameters to estimate. One way to do this would be to constrain the water advection parameter $V$ and the biomass diffusion parameter $D$ so that the the wavelength $\lambda$ and migration speed $c$ of the patterns agree with typical observed values. Based on observations by Deblauwe et al. (2012), we use $\lambda=10-200$ m and $c<60$ cm/yr and assume onset of patterns occurs with MAP of 120-350 mm/yr. Figure 9 shows that this constrains both parameter to approximate ranges of $V\sim 5-3500$ m/year and $D\sim 1-70$ m2/year for trees. The ranges for grasses are more constrained: $D\sim 1-15$ m2/year and $V\sim 5-15$ m/year. For comparison, Klausmeier (1999) took $D=1$ m2/year and $V=365$ m/year and typical values for the analogous transport parameters used in the Rietkerk modelare surface water flow speed of $V=3700$ m/year and biomass diffusion rate of $D_{b}=37$ m2/year. We note that typical water flow speeds are actually on the order of meters per second (i.e. $10^{7}-10^{8}$ m/year), so these values, based on a model with constant precipitation or “drizzle”, do not align with those following a rain event that might lead to the observed sheet flow as was described in Section 2.1.

The above estimates of the transport parameters rely on predictions about wavelengths and migration speeds of spatially periodic vegetation patterns at the onset of small amplitude patterns, i.e. at the Turing-Hopf bifurcation point where the spatially uniform state loses stability. However, the observed patterns correspond to large amplitude states in the model and appear in parameter regimes where a range of wavelengths are stable. This predicted multistability can be summarized by the so-called “Busse balloon”, originally developed in the context of Rayleigh-Bénard convection (Busse, 1978). These balloons, for example in the $(P,k)$ parameter plane, are based on nonlinear calculations that characterize the boundaries of stability for patterns as a function of the precipitation level $A$ and the pattern wavelength ($2\pi/k$). Such calculations, which are a nontrivial undertaking, have been carried out for the generalized Klausmeier/Gray–Scott model (van der Stelt et al., 2013). Recent progress has been made by generating observational Busse balloons to compare with predictions from this model (Bastiaansen et al., 2018). For a given site in Somalia (i.e. with fixed aridity) the empirical “Busse balloons” of band migration speeds, and of biomass level, are estimated as a function of the local wavelength of the pattern and the local elevation gradient. This allows for comparison between qualitative trends in the data and predictions from a model. For the generalized Klausmeier/Gray–Scott model, the observations are consistent with predicted trends in biomass and migration speed as a function of the pattern wavelength but do not match predictions about trends with local elevation grade.

4 Opportunities for Modeling: Time-scale Separation and Heterogeneity

Section 3.6 highlighted some of the challenges associated with modeling surface water flow using annually averaged quantities. Modeling these processes on the appropriate fast hydrological timescale means that model parameters correspond more closely to measurable quantities. Capturing fast processes presents a different set of challenges in terms of model complexity and numerical computation time. For example, a study that uses a process based ecohydrological model can make detailed predictions about water dynamics on short timescales (Paschalis et al., 2016). However, the simulations are limited by computational complexity to a few decades using static vegetation patterns on an idealized hillslope. Section 4.1 summarizes some mathematical techniques that exploit separation of scales in models to simplify analysis and/or generate reduced descriptions. Section 4.2 describes the important challenge of capturing rainfall as a stochastic driver in models of dryland vegetation patterns. Section 4.3 describes ways that spatial heterogeneity, due to large-scale topographic variations, might provide an opportunity for better connecting observations with models.

4.1 Scales in Time and Space

Many of the reaction-advection-diffusion models described in Section 3 are formulated on the timescale of biomass growth and death and capture emergent dynamics on a century scale. In this section we highlight progress and potential for capturing processes on appropriate scales within multiscale models and then reducing them using formal mathematical techniques. The advantage, when such reduction is possible, is that it provides a framework for capturing the (sometimes unintuitive) interactions between the fast time and short lengthscale processess on the long time and large scale dynamics.

Because water and biomass evolve on different temporal and spatial scales in dryland ecosystems, mathematical methods such as “geometric singular perturbation theory” that rely on separation of scales have been applied in conceptual models by separating the fast dynamics associated with water transport from slow dynamics of biomass dispersal (Bastiaansen and Doelman, 2018; Carter and Doelman, 2018; Sewalt and Doelman, 2017; Siero et al., 2015). In contrast to Turing-type analysis that can identify existence and stability of small amplitude (weakly nonlinear) patterns near onset of an instability, these methods apply for large amplitude (strongly nonlinear) patterns far from a pattern forming instability. The approach is to separately compute the fast (or small-scale) dynamics and the slow (or large-scale) dynamics in reduced systems and then connect these solutions together. The key to developing reduced descriptions in this approach rests on the result that the solution constructed in the limit that the scales do not interact persists when the interaction between scales is small but nonzero.

Multiple scales methods are already used for dealing with the timescale separation between water transport and biomass dispersal which are transport terms in the model. Such methods could also be used to capture the timescale separation that occurs within the interaction terms, between the slow biomass growth and death rates and the fast water infiltration rate. In the context of geometric singular perturbation theory, this would suggest using a pair of reduced models: one for during and shortly after rain events in which biomass is fixed and the water fields evolve on a fast time, and one for after the surface water completely infiltrates into the soil in which soil moisture and biomass evolve on a slow timescale. The complete dynamics could then be approximated by “gluing” the solutions from these two systems together.

4.2 Stochastic Drivers

Noy-Meir (1973) defined arid and semiarid ecosystems as “water-controlled ecosystems with infrequent, discrete, and largely unpredictable water inputs”. Indeed, the timing and magnitude of rainfall events in these areas has a large random component (Nicholson, 1980; Noy-Meir, 1973; Rodríguez-Iturbe and Porporato, 2007). While energy inputs (e.g., radiation and temperature) exhibit a fairly continuous variation over the year, rainfall usually comes in discrete “pulses” of very short duration (Huxman et al., 2004; Noy-Meir, 1973; Rodríguez-Iturbe and Porporato, 2007; Schwinning and Sala, 2004; Schwinning et al., 2004) (see Figure 10.) The intermittent and unpredictable character of rainfall results in fluctuations in soil moisture availability, which controls vegetation and ecosystem functioning (e.g., resistance to drought, productivity, germination). The response of the system depends on the magnitude of each pulse as well as the time interval between rainfall events, namely rainfall depth and frequency (Rodríguez-Iturbe and Porporato, 2007).

The role of stochastic rainfall variation in vegetation pattern formation is still poorly understood, thus opening up a number of possibilities for the development of modeling frameworks that explicitly account for the role of environmental stochasticity (Borgogno et al., 2009; Rohani et al., 1997). Moreover, predicted changes in rainfall frequency and intensity under climate change (IPCC, 2013; Pascale et al., 2016; Seager et al., 2007; Solomon et al., 2009) may dramatically impact dryland ecosystems. Understanding how these environmental fluctuations impact vegetation self-organization will help us assess the potential impact of a changing climate on these vulnerable ecosystems, and their possible desertification  (Rietkerk et al., 2002; Kéfi et al., 2007; Paschalis et al., 2016).

In some cases, vegetation models of the type described in Section 3 have been modified to include random environmental drivers (Kletter et al., 2009; Siteur et al., 2014a) or augmented by coupling to climactic conditions (Konings et al., 2011) as a way of studying the influence of stochastic rain ittermitancy and drought on deterministically formed patterns.

Noise-induced pattern formation (see Sagués et al. (2007) for a review from a physics perspective) has also been posited alternate mechanism for vegetation patterning. D’Odorico et al. (2006) show that the random switch between wet and dry conditions associated intermittent rainfall may trigger pattern formation even in the absence of a pattern forming bifurcation. In fire-prone savannas, the emergence of vegetation patterns is explained by purely noise-induced mechanisms linked to fire randomness and fire-vegetation feedbacks (D’Odorico et al., 2007).

4.3 Topographic Heterogeneity

Environmental heterogeneity occurs on many different scales and its influence can often be observed through characteristics of the vegetation patterns. For instance, topographic variation on the scale of an individual plant, microtopography, has been investigated in connection to a transition from banded vegetation patterns to vegetation following irregular drainage patterns by McGrath et al. (2012). Their model predicts that microtopography can disrupt the formation of vegetation bands and lead to an irregular drainage-type pattern of vegetation aligned along the slope. Other model studies have additionally sought to capture the influence of vegetation patterns in the formation of terracing of the underlying terrain along the hillslope through a coupling with erosion and soil transport dynamics (Baartman et al., 2018; Saco et al., 2007).

The pattern formation models of Section 3 are posed for ecosystems that evolve over time and space deterministically. However, because the timescales for the biomass birth and death processes are so long, it is challenging to test the model predictions through observations via time-lapsed images. An alternative to testing predictions against time-series is provided by the natural environment with its topographic heterogeneity. As described in Section 2, topography has long been known to play an important role for vegetation patterns.

Topographic data at 30-90m resolution (Farr et al., 2007) is available worldwide, which is a comparable spatial scale to some of the individual vegetation bands. Can we exploit this natural environmental driver, with all its varied realizations over space, to test and improve the models? Arcing of vegetation bands, described below, provides an example of how this might proceed.

The earliest observational study from the Horn of Africa (Macfadyen, 1950b), noted that the vegetation bands tend to align transverse to the prevailing elevation gradient and arc convex upslope (Figure 11a,b). Later studies (Boaler and Hodge, 1964; Greenwood, 1957) point to examples where the vegetation bands were arced convex downslope and suggest that the arcing direction depends on the location of the patterns relative to cross-slope curvature in the terrain: the bands tend to follow elevation contour lines (Figure 11c). Modern observations (Penny et al., 2013) support this claim that the bands tend to follow contours, but find that small-scale disruptions in topography or water flow can lead to significant deviation. The first modeling efforts that rely on ideas from pattern formation (Lefever and Lejeune, 1997), suggested that arced vegetation segments can form spontaneously on a uniformly sloped hill (with straight elevation contour line); similar emergent arcing on a uniform slope was found in other models as well (Meron et al., 2007; Siero et al., 2015).

Within the context of conceptual reaction-advection-diffusion pattern formation models, Gandhi et al. (2018a) made predictions about arcing direction by introduce a general topography $Z=\zeta(X,Y)$. They replaced the transport term $V\partial_{X}W$ in (4) with one based on a topography-dependent flux,

[TABLE]

in order to investigate the impact of a cross-slope curvature on patterns. From (9), we see that the transport leads to both an advective contribution ($\nabla\zeta\cdot\nabla W$), and a change in the (local) water loss rate proportional to the terrain curvature ($(\nabla^{2}\zeta)W$). On “ridge lines”, where $\nabla^{2}\zeta<0$, loss is increased due to water diversion, while in “valley lines”, where $\nabla^{2}\zeta>0$, loss is decreased due to water accumulation. The prediction from this model is that, as aridity increases, patterns shift from being predominantly on higher elevation ridge lines, convex downslope, to being spatially confined to lower elevation regions, convex upslope (as we see in Figure 11b). These model predictions find some traction when comparing arcing direction of vegetation bands in different regions that have different aridity levels (Gandhi et al., 2018a). This modeling framework makes another, as yet untested, prediction: all bands migrate uphill, regardless of their curvature, and low-lying valley patterns with the same band spacing, are predicted to migrate faster than those on higher ground, all else being equal.

5 Outlook

People study patterns because they are beautiful, but also because the semiarid regions where patterns occur sustain human life. There may be information about the fragility of the system embedded in the patterns. Studying the patterning and vegetation patchiness that arises in these areas may also offer insights into the efficient use of resources in these water-limited places. Changes in vegetation cover with changing climate also have the potential to feedback into the climate system as a forcing Rietkerk et al. (2011); Boussetta et al. (2015). Understanding land use, land use change, and vegetation evolution in patterned regions thus has implications for overall modeling of the climate system.

Climate change is predicted to pose a challenge in many drylands (Huang et al., 2017), expanding their extent and altering patterns of precipitation. This has implications for how well drylands can continue to support populations that depend upon them. Climate change projections for the Sahel and the Horn of Africa contain uncertainty, with models predicting different signs to rainfall change in some locations (Rowell et al., 2016).

Human interactions with the ecosystem include, but are not limited to, direct harvesting, to impacting grazing pressure, and to land use change. Other impacts (such as disrupting soil crust) have a number of human-driven causes (Belnap and Eldridge, 2001) and may be more difficult to parameterize in a model. An example of this disruption can be seen in the extensive tracking shown in Figure 12a. This figure, taken from Gowda et al. (2018), compares bands from 1952 and 2011. In the 2011 image, tracks are likely from vehicles. This soil disruption has spread vegetation out on the landscape and interrupted stripes. Human actions can also generate regions of fertility in these systems. For instance, creating patches has been used for restoration in a semiarid area (Ludwig and Tongway, 1996).

Land use choices also impact the extent of land that may exhibit patterns. For example, in Figure 12b, near Fort Stockton, Texas in the United States, rain-fed vegetation patterns are interrupted by irrigation driven agriculture. The circular regions that are more highly vegetated are due to irrigation, while the banded patterns nearby are the local vegetation cover. There has also been an expansion of cropland in the Sahel over the satellite era (Hiernaux et al., 2016). Changing land use, perhaps evolving with population pressure, sits outside the scope of current modeling efforts for vegetation patterns.

There is a disconnect between the reality of human-impacted systems, and the types of idealized, ecosystem-only models, that are straightforward to analyze. Human pressures on biomass, on soil, and on land use, are difficult to parameterize in models. In addition, human impacts have downstream effects on human livelihoods, so there is likely a feedback between human pressure and biomass productivity that is not captured in modeling work that neglects humans. In the future, data driven modeling approaches combined with high-throughput remote sensing data have the potential to change how we frame this problem. Observation-based work intrinsically includes human interactions with the system. The growing back-and-forth between observations and modeling makes this a rich area of study.

Acknowledgements: We are grateful to Karna Gowda for insightful discussions and, along with Jake Ramthun, for help with figures. This work was supported in part by the National Science Foundation grant DMS-1440386 to the Mathematical Biosciences Institute (PG) and National Science Foundation grant DMS-1517416 (MS).