Research articleLocal properties of patterned vegetation: quantifying endogenous and exogenous effects
Abstract
Dryland ecosystems commonly exhibit periodic bands of vegetation, thought to form due to competition between individual plants for heterogeneously distributed water. In this paper, we develop a Fourier method for locally identifying the pattern wavenumber and orientation, and apply it to aerial images from a region of vegetation patterning near Fort Stockton, TX, USA. We find that the local pattern wavelength and orientation are typically coherent, but exhibit both rapid and gradual variation driven by changes in hillslope gradient and orientation, the potential for water accumulation, or soil type. Endogenous pattern dynamics, when simulated for spatially homogeneous topographic and vegetation conditions, predict pattern properties that are much less variable than the orientation and wavelength observed in natural systems. Our local pattern analysis, combined with ancillary datasets describing soil and topographic variation, highlights a largely unexplored correlation between soil depth, pattern coherence, vegetation cover and pattern wavelength. It also, surprisingly, suggests that downslope accumulation of water may play a role in changing vegetation pattern properties.
1. Introduction
Pattern formation occurs in numerous ecological and biological systems, where it has been linked to reaction–diffusion (Turing)-type instabilities [1], hydrodynamic instabilities [2] and potential (variational) dynamics that maximize or at least increase ecological productivity [3,4]. Patterns have been observed in dryland vegetation [5], in bogs and wetlands [6,7], in mussel beds [8], termite mounds [4] and other systems [9]. In many cases, ecological patterns form an intermediate realization between environmental states in which the entire landscape is either colonized or bare. As such they often indicate the presence of bistable states, which are characterized by the potential for critical and locally irreversible transitions [10]. Both observations and theoretical work suggest that transitions between vegetated and desertified (unvegetated) conditions in patterned systems are preceded by striking changes in the morphology of the vegetation patterns, which thus act as an early-warning sign of deteriorating ecosystem health [11]. Understanding the controls on the morphology of vegetation patterns therefore has practical interest in terms of ensuring that observed changes are interpreted correctly.
Vegetation pattern formation in dryland ecosystems is a global phenomenon, ranging from random distributions of bare soil and vegetation canopies [12] to highly organized spatial distributions with identifiable length scales and orientations [1,5,13]. Intermediate cases such as power-law (scale-free) clustering [14–16], and dendritic structures in which vegetation concentrates along drainage lines [17,18] are also observed. All of these morphologies can be related to the presence, strength and directionality of positive feedbacks that concentrate resources (such as nutrients, soil carbon and water) that sustain plant life in a localized region near the plants [19–24]. These feedbacks have led to the moniker ‘ecosystem engineers’ being applied to perennial plants in dryland ecosystems: they create the conditions necessary for their own survival [25,26].
The formation of periodic vegetation patterns is strongly linked to the coupling of (i) water redistribution to the vicinity of plants and (ii) competition between individual plants for access to this water resource [3,5,27,28]. The dynamics of patterned vegetation systems has formed a focus of ecohydrological and nonlinear dynamics research, motivated by the increasingly well-demonstrated connection between pattern morphology and desertification risk [1,11,29,30]; and by the inherently interesting nonlinear dynamics of these systems [3,31]. Field and theoretical work has identified the important roles of vegetation dynamics [11], seed dispersal processes [32,33], root morphology [3,29], surface flow dynamics [34] and climate feedbacks [35] in modulating pattern morphology.
While our understanding of vegetation pattern dynamics is improving, there remain undeniable differences between simulated vegetation patterns and the natural ones observed in the field. Natural patterns are characterized by a degree of disorder and heterogeneity on multiple scales that is not reproduced in idealized models. Disorder can arise intrinsically through the existence of a range of stable wavelengths and orientations, with transitions in space and time occurring through local pattern defects [36]. In addition, the pattern wavelength and orientation can be strongly influenced by either wavelength-scale inhomogeneities [37] or larger gradients [38] in the underlying system, as well as the presence of boundaries [39] or noise [40].
In natural systems, this means that disorder observed in vegetation patterns could either reflect intrinsic features of the pattern-forming process (endogenous effects), or could reflect spatial changes in soil structure and local topography (exogenous effects). For instance, Thompson et al. [32] explored whether the unrealistically smooth nature of many models of vegetation biomass distribution was an artefact of representing seed dispersal as a diffusive process. Greater fidelity between modelled and observed (disordered) vegetation patterns was achieved by representing plant population migration with a seed dispersal kernel. More simply, theoretical treatments of vegetation patterning usually neglect variations in soils and topography, or impose periodic boundary conditions that remove the differences in water availability between the top and bottom of a hillslope. Such simplifications will inevitably lead to idealized representations of the vegetation response and obscure interactions between the intrinsic pattern structure and the spatial structure of the landscape. For example, McGrath et al. [17] demonstrated that the orientation and wavelength of vegetation patterns modelled with realistic boundary conditions changed between the top and bottom of the hillslope. In control simulations with periodic boundary conditions, vegetation bands formed with a single wavelength and were orientated at 90° to the direction of water flow, in agreement with other modelling studies. When realistic boundary conditions were imposed, however, the bands near the top of the hillslope curved to lie perpendicular to the ridgeline. Similarly, the effects of spatial heterogeneity in soil properties on pattern formation have not been widely explored (although see [32]). A number of theoretical studies indicate that local biomass, band properties and band coherence should vary with changing minimum and maximum infiltration capacities and soil properties [41,42], and there are some tantalizing hints that subsurface features, such as the ironstone underlying tiger bush in Niger, calcrete hardpans underlying banded patterns in Texas, and silcrete hardpans underlying mulga bands in Australia [43–45], may have an association with patterned vegetation.
Predictions about the interaction of vegetation patterns with changing soil or topography can be investigated using remotely sensed datasets, an approach with a long and growing history. Vegetation patterns in Africa were first observed from light aircraft flights [46,47]. Initial analyses of the patterns demonstrated their spatial regularity on the basis of the two-dimensional Fourier power spectrum [48]. More recently, Deblauwe et al. [30] undertook large-scale analyses of morphological trends in vegetation patterns in the Sudan by linking several remotely sensed datasets: surface imagery from the System for Earth Observation (SPOT) satellite, topographic data from the Shuttle Radar Topography Mission (STRM), and rainfall data from the Tropical Rainfall Measuring Mission (TRMM). This allowed analysis of pattern morphology at 400 m resolution over multiple square kilometres. Although higher resolution datasets are available, they have generally either been analysed only over relatively small spatial scales (≈1 km2) [29] or investigated from the perspective of identifying morphological change over time [49].
A key open question is therefore to characterize small-scale irregularities in vegetation patterns, and to classify them based on whether they arise due to either endogenous dynamics or variation in exogenous features imposed by the landscape. Understanding the implications of such variation on the resilience and stability of ecosystems has important consequences for identifying and preventing potential desertification.
In this paper, we report on the development of methods suitable for quantifying local patterns within high-resolution aerial photography, and for relating those features to ancillary datasets describing soil and topographic variation. Our site, located near Fort Stockton, TX, USA [43], was selected because it exhibits vegetation bands and also offers several characteristics that facilitate studies of covariation. First, the aerial photography covers a large region at a resolution (0.5 and 1 m) comparable to the diameter of perennial vegetation canopies, allowing for the observation of individual plants. Second, high-resolution (10 m spatial, 15 cm vertical) elevation datasets are available from the US National Elevation Dataset (NED), allowing changes in the orientation and gradient of the hillslope to be mapped on scales that are much smaller than the pattern wavelength. Finally, the study area has been mapped as part of the Soil Survey Geographic (SSURGO) national soils database, and contains considerable variation in soil type. The availability of these datasets provides a valuable opportunity to investigate correlations between vegetation patterning and soil characteristics over tens of square kilometres.
Our analysis depends on these spatial datasets, and is subject to their inherent limitations. These limitations include the resolution, spatial artefacts and the risk of spurious correlation, given that vegetation features are often considered when mapping soil boundaries. To minimize the effects of these issues, we focus on two large-scale hypotheses: (1) the wavelength and orientation of the vegetation pattern are locally coherent but exhibit both rapid and gradual variation and (2) the variability in vegetation pattern features will correlate with soil and elevation features in predictable ways. We note that by focusing on soil and topographic features, we inherently assume that spatial variations in pattern morphology arise due to spatial heterogeneity in landscape properties, rather than due to the nonlinear dynamics of the pattern-forming process itself. In patterns far from threshold, features such as defects, dislocations, grain boundaries and boundary conditions can result in spatial heterogeneity in pattern properties, even when all other fields are homogeneous [50]. To control for this possibility, we also briefly address the following null hypothesis: (3) large-scale variations in pattern properties can be explained by the nonlinear interactions associated with vegetation pattern formation.
We test the major hypotheses through the development of a localized Fourier technique to identify pattern wavenumber and direction. To address Hypothesis 1, this technique was applied to high (0.5 m) resolution aerial photography. To address Hypothesis 2, the technique was applied to a 188 km2 area, allowing local pattern properties to be correlated with the site soil and topographic properties. Hypothesis 3 was tested by identifying a subset of the 0.5 m resolution image containing banded patterns that were close to ideal (i.e. reproducible by a model). We applied the Fourier analysis technique to both observed and simulated patterns, and compared variability in the wavenumber and direction fields in the modelled and observed patterns.
We find that our two major hypotheses are satisfied. Local patterns are oriented in the same direction as the topographic slope and the pattern wavelength decreases for steeper gradients. Deviations from these trends are associated with the presence of ridges, stream channels, anthropogenic features or changes in soil type. Different soil types within the study area determine the pattern boundaries and the pattern morphology: shallow soils are associated with highly coherent, shorter wavelength patterns, and deep soils with patterns that become incoherent and increase in wavelength near stream channels. Finally, the modelling test validated our decision to focus on landscape and soil features. Even in the most uniform region of the observed patterns, the variability in the real pattern properties exceeded that which could be retained in the steady-state solution of a physical model that closely reproduced the mean pattern properties.
2. Material and methods
(a) Study site and data
The study site is a 188 km2 area located approximately 30 km northwest of Fort Stockton, TX, USA (coordinates: 31°05′ N, 103°03′ W) [51]. The climate is hot and dry, with 370 mm annual rainfall on average, mean summer maximum temperatures of approximately 35°C and winter minima near freezing. The site is part of a large cattle ranch and is used for grazing. Dominant vegetation species include tarbush (Flourensia cenua), bunch and sod grasses (Aristida purpurea, Bouteloua curtipendula and Scleropogon brevifolius), and mixed mesquite (Prosopsis glandulosa) and juniper (Juniperus pinchotti) brush [43]. The site contains a striking spatial pattern consisting of bands of continuous vegetation cover lying over bare soil [43] (figure 1).
Figure 1. (a) Binary image of the Fort Stockton region, covering an area of approximately 188 km2. Vegetation is shown as black and bare soil is shown in white. (b) Detail of size 260×260 m2, showing the original detailed image. (c) Binary image detail. (Online version in colour.)
High-resolution imagery (0.5 and 1 m pixels) were obtained from Digital Globe, and from the National Agricultural Imaging Project (NAIP) [52]. We used the highest resolution images from Digital Globe for fine-scale analysis and a comparison between modelled and observed pattern morphology. We analysed the NAIP images, which cover the whole area, to relate local pattern properties to soil and topographical features.
We classified the image pixels as ‘vegetation’ or ‘no vegetation’ using a supervised classification based on total brightness [53]. We used brightness because the perennial vegetation was not actively photosynthesizing when the photographs were taken, meaning that standard vegetation indices could not discriminate the vegetated locations. An example of the resulting binary image is shown in figure 1. Darker colours represent vegetation cover and the lighter colours bare soil. The insets show an original and classified image over a 260×260 m2 window. We obtained elevation data from the NED provided by the US Geological Survey and soil maps from the SSURGO Database.
(b) Fourier windowing method
To quantitatively test our hypotheses, we developed a quasi-local technique to obtain spatial maps of pattern wavelength λ and pattern orientation Θ for the binary images. The technique provides information about the local wavevector 
, similar to that provided by short-time Fourier transforms or wavelet-based approaches used in timeseries analysis [54,55]. Local wavevectors are useful for classifying convection patterns [56] and for identifying pattern defects [57,58]. Here, we applied a two-dimensional Fourier transform to obtain the power spectrum within a square, moving window. Local wavelength and pattern orientation were identified for each window, and the results averaged for all windows that spanned a given location. This straightforward technique is suitable for identifying the dominant pattern properties in noisy images with irregular patterning. More elegant and rigorously local techniques, based on the ratios of the spatial derivative of the banding pattern [57], could not be applied because the vegetation bands deviate substantially from a sinusoidal pattern in space. The drawbacks of the short-time Fourier transform—namely the tendency to truncate long wavelengths due to the finite window size, and poor localization of short-wavelength components [59]—do not pose significant difficulties in the current application, provided that the window is appropriately sized with respect to the local wavelength of the pattern.
For each window of size L×L, we obtained the two-dimensional fast Fourier transform 
of the pattern f(x,y). For each window’s 
, we calculated the power spectrum 
. As L increases, the likelihood of identifying a single wavelength and orientation decreases. Conversely, as L decreases, the k-resolution in Fourier space is reduced. To optimize both the localization of measurements and the resolution, we chose L to be at least 5λ. The window was applied to overlapping regions of the pattern at 20 m intervals.
The power spectrum measures the power contributed to the pattern by each wavevector k. To separate the local wavelength from its orientation, we decomposed each k into its magnitude (wavenumber) 
and its orientation θ. Owing to symmetry, the pattern orientation is defined only between 0 and π. To identify the dominant k in each window, we binned S(k) into annular rings of width k=6π/L. To deconvolve the natural 1/k scaling of the image [60,61], we computed the total power within each annular ring, S(k). The location peak of this total power (rather than the mean power) is used to define the most energetic wavenumber, k1. To compensate for the large bin width, we computed the location of the weighted average k1≡〈kS(k)〉/〈S(k)〉 over all rings that formed part of the peak and contained more than 75% of the peak power. We discarded windows where k1 corresponded to the limits of the window function or pixel resolution, and, to avoid discarding sites where L truncated λ, visually confirmed that these windows corresponded to unpatterned areas. Finally, we discarded windows where the mean power of the pattern-forming wavenumber was less than that of noise. To determine the dominant pattern orientation, we binned S(k) into segments of width π/8 and computed the average power, S(θ). The dominant orientation was located via the weighted average θ1≡〈θS(θ)〉/〈S(θ)〉. The most energetic values (k1,θ1) for each window were used to generate spatial maps λ(x,y) and Θ(x,y) representing the local pattern wavelength and orientation. For each x,y location, the assigned value of λ and Θ comes from the average 2π/k and θ of all windows which include that x,y.
Although the procedure so far assumes a two-dimensional pattern with a single wavelength and orientation within each window, the power spectra regularly contained additional peaks. We evaluated the uniqueness of the local pattern properties by comparing the dominant peak to its most distant energetic peak. Energetic peaks were defined as those containing more than 75% of the power of the dominant peak. The most distant peak was then defined as the energetic peak with wavenumber k2 and orientation θ2 that maximized |k1−k2| and |θ1−θ2|. Uniqueness metrics were then defined for both the wavenumber and the orientation as
Matlab code that performs the Fourier windowing method outlined in this section is provided as the electronic supplementary material.
(c) Construction of the spatial fields used for analysis
The local Fourier analysis was run twice, using windows with L=260 m and 400 m, applying the window to overlapping regions on 20 m intervals, and generating smoothly varying maps of the local pattern properties. The two window sizes ensured that the largest wavelengths could be analysed (L=400 m), and exploited the greatest feasible localization given the median pattern wavelengths (L=260 m). The L=400 m method identified large λ in regions where the L=260 m method identified no pattern; while the L=260 m method allowed the pattern properties to be identified in regions with dimensions less than 400 m. To combine the outcome, we averaged the results of the analyses for overlapping cells (resulting in a minimum change in the identified Θ and λ), and retained the uniquely identified Θ and λ in cells where either mechanism failed. Since the pattern properties in windows located near the edges of the pattern are influenced by both patterned and unpatterned regions, we trimmed the resulting λ(x,y) and Θ(x,y) fields by 130 m to discard the most strongly affected regions. The results of the complete analysis are shown in figure 2.
Figure 2. Maps of (a) local wavelength λ(x,y), with colour scale indicating wavelength in metres and (b) local orientation Θ(x,y), with colour scale indicating orientation in radians. Areas in white do not contain patterns as recognized by the Fourier windowing method. (Online version in colour.)
We developed a map of the soil type on the 20 m grid by interpolating the SSURGO data. We smoothed the NED elevation data to remove high-resolution artefacts [62], and computed the topographic gradient (steepness) and aspect (orientation) over the same 20 m grid. We made direct comparisons between λ(x,y), Θ(x,y) and the soil and topographic features. To cope with the large, noisy dataset, we grouped the data into bins of equal size and related the median and interquartile range of the local pattern characteristics to the predictor variables across the bins. The behaviour of these summary statistics and the dataset was analysed using least-squares regression. We analysed the frequency of occurrence of deviations between the pattern and hillslope orientations ΔΘ and the frequency of occurrence of non-unique pattern orientation and wavelength, QΘ and Qλ, as a function of location (e.g. near topographic minima (drainage lines), maxima (ridges), anthropogenic features (roads and trails) and the edges of the pattern).
(d) Comparison to models
To quantify the degree to which local pattern variations can be explained as non-ideal pattern features that can arise from the nonlinear dynamics of the pattern-forming mechanisms, we performed a control analysis using simulated data. Our test region is a 500×800 m2 region of the 0.5 m resolution pattern. Analysed using a 200×200 m2 window, the pattern in this region contained a single high-energy wavenumber and direction peak, indicating that it could be reasonably modelled with homogeneous model parameters. We used the pattern in this region as the initial condition for a physically based vegetation band-forming model [1,33]. We used this model because it (i) represents the surface runoff feedback mechanisms that were observed to occur at this site by McDonald et al. [43] and (ii) simulates vegetation bands that retain curvature, variability and other non-ideal behaviours (compared to more idealized models [63]) and thus has the potential to generate spatially variable pattern properties endogenously. We selected model parameters by stepping through reasonable parameter combinations and selecting the combination that minimized the variation between the observed vegetation pattern and the model prediction after 5000 timesteps [32]. This method allows us to identify a parameter set for the model that approximates the observed patterns as a stable steady-state solution. Using both the model output and the original binary images, we apply the same Fourier windowing method and compute key spatial statistics (mean, variance and autocorrelation length, taken as the lag at which the autocorrelation halved) for the resulting λ(x,y) and Θ(x,y).
3. Results
(a) Local properties of vegetation patterns
Using the Fourier windowing method, we find that 44% of the 188 km2 area contains patterned vegetation with a clearly identifiable wavelength and orientation. The Fourier windowing method fails to detect some areas where vegetation patterns are visually identifiable. These include regions where vegetation bands are confined to fingers of a particular soil type that are 
m in extent, meaning that the pattern cannot be identified as the dominant Fourier mode in any given window. There are also some regions where a pattern can be identified by eye, but is so disordered that it falls below the noise threshold. Approximately 10% of the study area consists of isolated patches of patterning which offer little opportunity to explore spatial variations. Instead, we confine our analysis to the 34% (≈64 km2) of the image that consists of spatially connected vegetation patterns. Within these regions, the pattern has a mean wavelength of 
m (reported variations are standard deviation unless otherwise specified) and an autocorrelation length of approximately 800 m (i.e. ≈12λ). The probability distribution function (PDF) of Θ(x,y) values is peaked in the north–south direction, but spans a full 0≤Θ<π range. Table 1 provides summary statistics describing the average pattern characteristics and topography.
| mean | s.d. | corr. length (m) | |
|---|---|---|---|
pattern wavelength,  (m) |
63 | 13.8 | 800 |
mean pattern orientation,  (rad) |
1.4 | 0.76 | 600 |
| hillslope orientation (rad) | 1.2 | 0.75 | 1400 |
| hillslope gradient | 0.007 | 0.002 | 1900 |
While seven distinct soil types occur in the study area, two of these, the Delnorte and Reakor associations, contain 97% of the vegetation patterning. These soils are distinguished by a relatively low clay and high silt content. The Delnorte association is characterized by a very shallow soil depth (23 cm) due to the presence of calcium-carbonate-based hardpan or petrocalcic horizon. By contrast, Reakor association soils are at least 2 m deep. The full range of soil properties is shown in table 2.
| Reakor | Delnorte | |
|---|---|---|
| soil depth (cm) | 200 | 23 |
| hydraulic conductivity (Ksat, μm s−1) | 2.1 | 59.2 |
| clay content (%) | 31.5 | 8.7 |
| sand content (%) | 6.8 | 42.8 |
| silt content (%) | 61.7 | 21.1 |
We find that three factors broadly determine the local pattern properties: the orientation of the slope, the steepness of the slope and the soil type. The orientation Θ of the pattern is almost completely determined by the underlying slope orientation. As illustrated in figure 3a, approximately 81% of the pattern is oriented within ±π/8 rad of the hillslope. Deviations between the pattern and hillslope orientations, denoted ΔΘ, are discussed further in §3c. Second, we observe a significant relationship between the local pattern wavelength λ and the hillslope gradient. Figure 3b illustrates the moderately strong, significant decline in the median wavelength λ of equally sized data bins with increasing hillslope gradient (r2=0.63,p=6×10−3). In addition, we found that the pattern wavelength is influenced by the soil type. Figure 3c shows PDFs of λ for each soil type. While Delnorte soils have a unimodal PDF (mode λ=53 m), the Reakor association soils have a bimodal distribution. One mode corresponds to that of the Deltnorte soils (λ=52 m), but the other modal wavelength is much longer, with λ=71 m. Further analysis of this effect is provided in §3d.
Figure 3. (a) Frequency distribution of ΔΘ, the deviation between the local pattern orientation and the local topographic orientation. The slope orientations deviate by less than π/8 for 80% of the pattern. (b) Comparison of the local pattern wavelength λ to the local topographic slope. Data were collected into 10 bins of equal size (≈11 000 data points per bin). The box and whisker plots indicate the median (central lines); the 25th and 75th quartiles (box limits) within the bin. The whiskers extend two inter-quartile ranges from the median. The thick line indicates a linear fit to the medians. (c) PDFs of the pattern wavelengths associated with the two major soil classes, indicating the bimodality in λ on Reakor association soils. (Online version in colour.)
(b) Comparison to model
The results are shown in table 3. We find that while the mean properties of the model and observed pattern are similar (a consequence of calibrating the model to these means), Θ and λ vary three to five times more in the observed pattern than in the modelled pattern. The underlying topography is still more variable, and presumably causes the variability of the observed pattern. Thus, even in the region of the study site with the greatest uniformity in the pattern properties, the spatial heterogeneity observed in the real patterns exceeded the pattern variability that could be produced by a physical model. These results justify our attribution of the remaining variability in the pattern properties across the site to environmental variation rather than to defects or initial condition effects arising from the nonlinear dynamics of the system.
| mean | s.d. | corr. length (m) | |
|---|---|---|---|
natural pattern wavelength  (m) |
53 | 4.1 | 350 |
model pattern wavelength  (m) |
54 | 1.7 | 360 |
natural pattern orientation  (rad) |
1.7 | 0.39 | 420 |
model pattern orientation  (rad) |
1.6 | 0.07 | 390 |
| hillslope orientation (rad) | 1.7 | 0.5 | 120 |
(c) Non-uniqueness in pattern attributes
While the local pattern wavelength and orientation are on average set by the local hillslope gradient and orientation, we nonetheless observe regions with significant deviation from the overall trend. As shown in figure 3a, the deviation between pattern and hillslope orientation, ΔΘ, is more than π/8 over approximately 20% of the pattern. There are even regions in which ΔΘ≈π/2, where the vegetated bands run parallel rather than perpendicular to the local hillslope gradient. Figure 4 shows the spatial distribution of the orientation uniqueness metric QΘ across the patterned region, highlighting regions of non-uniqueness. About half of the windows with large ΔΘ also contain more than one pattern direction (non-uniqueness in Θ(x,y)). As illustrated in figure 5, non-unique pattern orientations are clustered near streams (50% increase in frequency of QΘ<0.75 relative to the remainder of the pattern), ridges (50% increase in frequency), roads (30% increase in frequency) and the pattern edge. We also observe regions with large deviations (ΔΘ>π/8) but nonetheless unique pattern orientations (QΘ≥0.75). Such regions also occur near ridges and streams more frequently than in the remainder of the pattern (≈30% more frequent in each case).
Figure 4. Spatial distribution of the uniqueness metric QΘ for the pattern orientation, and magnified cases illustrating different examples of non-uniqueness in the pattern orientation. (a) Non-uniqueness in orientation associated with a road, which adds north–south oriented structure to the local northeast–southwest oriented banding due to upslope ponding of runoff. (b) A region of changing pattern orientation and wavelength associated with the interleaving of two contrasting soil types. (c) The short length scales on which the pattern can change orientation as it bends around a sharp ridge. A proportion of the region where ΔΘ≈π/4, as indicated by the different director arrows, has QΘ>0.75. (d) The deviation between pattern and slope orientation around a stream channel; again there are regions where ΔΘ≈π/2 and QΘ>0.75. (Online version in colour.) Figure 5. Relative increase in frequency of (a,b) non-unique pattern orientations and (c) wavelengths, for windows within 250 m of local elevation minima (streams), maxima (ridges), roads, and pattern boundaries. All frequencies are measured relative to the frequency in the remainder of the pattern. Non-unique orientations, for example, occurred 50–60% more often near streamlines and ridges than in the bulk of the pattern. (Online version in colour.)
The insets in figure 4 illustrate regions of complexity in Θ(x,y). Figure 4a,b shows regions of non-unique pattern directions. In figure 4a, we show a perturbation of the local pattern by a road. As illustrated, the upslope edge of roads is globally associated with increased vegetation, while the downslope edges are typically bare. Roads thus generate linear features that can create non-uniqueness in the local pattern orientation. Figure 4a also shows a second anthropogenic feature, a storm drainage outlet, which further decreases QΘ west of the road. Figure 4b illustrates how changes in soil type can lead to low QΘ. In this example, fingers of the Reakor association soils are interleaved with the Delnorte association soils, causing rapid changes in pattern wavelength and orientation, and consequently multiple energetic values of θ and k within the 260 m windows. Figure 4c illustrates a region where the pattern orientation changes rapidly, turning through approximately π rad within a single 260 m window. Such rapid change inevitably leads to low QΘ because there are multiple pattern orientations located in a single window. There are also, however, regions near the ridge crest where QΘ≈1 and the pattern is locally oriented perpendicular to the slope. Figure 4d is representative of the complex vegetation transitions that occur near stream channels. Here, too, the pattern lies perpendicular to the local hillslope orientation, and rather than curving into the streamline, remains broadly aligned with the band orientation away from the stream.
(d) Soil-type effects
While the vegetation patterns on the Delnorte soils have a unimodal wavelength distribution, the Reakor soils exhibit a bimodal distribution (figure 3). The pattern associated with the second peak of λ on the Reakor soil consists of bands of bare soil within a matrix of vegetation cover, inverting the distribution in the remainder of the pattern. Visually, this ‘anti-stripe’ pattern is more disordered than the remainder of the pattern seen on the Reakor association soils. In the power spectra, these anti-stripe peaks contain less than half the power of the short wavelength, coherent peaks. Thus, on the Reakor association soils the pattern varies through space from short-to-long wavelength, ordered-to-disordered patterns, and lower-to-higher biomass. Because increased biomass in drylands implies an increased access to water, we examine how λ varies as a function of the distance to the nearest stream (as a proxy for availability of water). The results are shown in figure 6. While patterns on the Delnorte soils did not show a distance-dependent λ (r2=0.08 and p=0.4), the patterns on Reakor soils show a trend of increasing λ with decreasing distance to the nearest stream (r2=0.96 and p<5×10−7).
Figure 6. Distribution of the soils within the pattern-forming areas with respect to landscape position, referenced to the nearest stream channel or local elevation minimum for (a) Reakor and (b) Delnorte association soils. Although both soil types occur across the range of landscape positions, the Reakor association soils are more strongly associated with riparian areas. Distribution of pattern wavelengths on the (c) Reakor association and (d) Delnorte association soils with respect to the nearest stream, again based on equally sized bins of the dataset. There is a strong (r2=0.96) and significant (p<5×10−7) decline in the wavelength on Reakor soils when moving from the riparian areas to the uplands. There is no correlation between landscape position and pattern properties on the Delnorte association (r2=0.08 and p=0.4). (Online version in colour.)
4. Discussion
The local Fourier metrics confirm Hypothesis 1 by showing that the wavelength and orientation of the vegetation patterns are typically coherent on scales of 600–800 m, but can change on scales of ≈20 m when the slope orientation changes over the same scale. We were unable to reproduce the observed variation with a model even within regions of relatively uniform patterning, confirming Hypothesis 3 and further suggesting that exogenous factors rather than the pattern-forming mechanisms drive the spatial variability in λ and Θ. In the light of these findings, and to address Hypothesis 2, our subsequent analysis focuses on the relationship between exogenous factors and the local pattern properties.
Our analysis indicates a hierarchy of controls on the morphology of the vegetation patterning: at a global scale the pattern morphology is determined by the slope orientation and the hillslope gradient, while the soil type imposes the template for the pattern-forming region. At smaller scales, complexity in the form of multiple local pattern orientations and deviations from the global trends are associated with ridges, streams, roads and changes in the soil type. Organization of the pattern characteristics on small scales arises due to a combination of soil type and topographic context: homogeneous, well-defined patterns with a unimodal λ on Delnorte soils, regardless of the topographic context. On Reakor association soils, the pattern λ increases in locations closer to streams, while the pattern simultaneously becomes more disordered and vegetation cover increases.
Several of these relationships support previous observations. The average hillslope gradient of 0.7% at the Fort Stockton site conforms to observations and predictions of a minimum hillslope gradient being needed for band formation [41], and exceeds the minimum slope gradient observed in other settings (e.g. 0.25% in the Sudan [30]). In addition, the inverse relationship between pattern wavelength and the hillslope gradient (see figure 3) has been identified elsewhere [30,64–66]. However, the one model [67] that has been used to study wavelength–gradient correlations predicts the opposite trend: steepness causes larger wavelengths [41,68]. We explored this trend with a different model [1] and find that it also produces the opposite trend to that observed in natural systems. While both models predict that multiple pattern wavelengths are stable for a given hillslope gradient [33,69], the selection of particular longer or shorter wavelengths within that range is clearly at odds with observations.
The close correspondence between hillslope and pattern orientation at the Fort Stockton site is a near-universal feature of vegetation patterning. The observed deviations generally arise due to local effects that disrupt the pattern, rapidly alter the direction of water flow, or change the strength and/or length scales of the pattern-forming feedbacks.
Pattern disruption is exemplified by the effects of roads (see figure 4a). Roads disconnect upslope and downslope vegetation bands by preventing the redistribution of water. Rapid alteration of the direction of water flow occurs along ridges and stream channels. The pattern near these locations is less likely to have a unique orientation or wavelength than in the mid-slope areas, and is more likely to have a large ΔΘ. The prominent ridgeline shown in figure 4c is locally surrounded by a region with ΔΘ≈π/2 and QΘ>0.75. This location provides an unambiguous example of a change in pattern orientation near no-flow boundary conditions, as predicted by McGrath et al. [17]. Other ridgelines are associated with large ΔΘ, but the low QΘ in these locations makes the interpretation of the observed ΔΘ ambiguous. Methods that identify a truly local metric of pattern properties, instead of the quasi-local metric used here could both help to resolve these ambiguities, and to extend the analysis into patterned areas less than 260 m wide.
Stream channel locations are also associated with rapid, and sometimes discontinuous changes in pattern orientation. Like the ridge shown in figure 4c, the stream channel shown in figure 4d is surrounded by a region where ΔΘ≈π/2 and QΘ>0.75, i.e. a unique pattern orientation that is perpendicular to the one expected from the slope orientation. There has been little exploration of the interaction of vegetation pattern with stream channels, presumably because the current paradigm of vegetation pattern models offers little reason to think that such interactions would be important. Stream channels occur in locations where surface runoff rapidly flows away and therefore cannot affect vegetation upslope from the channel. However, we find evidence that the distance from a stream channel alters the pattern properties on the Reakor association soil type. On these soils, the pattern wavelength, vegetation cover and disorder increase near the streams. These observations suggest that an additional mechanism could affect the patterning at this study site: for instance, the stream-channel boundary condition propagating back up into the hillslope.
We hypothesize that the ultimate cause of the changes in pattern morphology on the Reakor soils lies in the contrast in the soil depth between the Delnorte and Reakor associations: from 23 cm to over 2 m. We are not the first to propose an association between soil depth and the vegetation pattern structure. Depth to a silcrete hardpan is associated with a transition between sharp patterns (shallow hardpan) and diffuse patterns (deep or no hardpan) in Australia [45,70]. Strikingly coherent vegetation bands in Niger occur above a shallow ironstone hardpan [44]. The broad, diffuse banding near Fort Stockton studied by McDonald et al. [43] occurs on deep clays. McDonald et al. [43] cited unpublished research claiming that vegetated bands were associated with a local increase in the depth to the hardpan, similar to previous observations of increased vegetation density on deep soils in Australia [71].
Field evidence is needed to determine the exact mechanisms by which soil depth and proximity to streams could lead to the changes we observed in pattern morphology. Three scenarios illustrating potential mechanisms are shown in figure 7 and could form the basis for future field studies. First, shallow soils could promote lateral root extension by plants, exaggerating the effects of root competition [25,26] (figure 7a,b). Studies of Chihuahuan desert species confirm that there is intense root competition in the zone above the petrocalcic horizon, where root growth is concentrated [72].
Figure 7. Schematic of possible effects of a shallow impeding layer on pattern-formation mechanisms. (a) In the absence of an impeding layer, plants develop deep root systems with a smaller lateral root-zone radius, Rr, to exploit water from across the whole soil profile. (b) A shallow impeding layer restricts the vertical growth of roots, confining them to the surface soils, and promoting root competition laterally. (c) The shallow impeding layer also creates the potential for complex subsurface hydrology. Saturation of soils above the impeding layer may promote saturation excess runoff. Gradients—due to changes in the depth of the impeding layer, or in water potential across the impeding layer—can also induce saturated porous media flow. Provided that surface soils near plants can freely drain, the altered hydrology may exaggerate the positive feedback between plants and water. (Online version in colour.)
In the second scenario, we recognize that shallow surface soils have limited water storage and saturate readily. For example, the Delnorte association has as little as 2 cm of water storage capacity in the soils above the impeding layer. Unlike dry soils, which only generate runoff during very heavy rains, saturated soils shed all rainfall as runoff. If soils did not saturate near plants this runoff water could be trapped and infiltrate in these locations (figure 7c). Plant roots penetrate and thus may break up hardpans [72]. Root water uptake is also a driver of hardpan formation [73], and hardpans might thus form at greater depth beneath deep-rooted (woody) vegetation. Either mechanism could prevent the surface soil from saturating near the vegetated bands and allow it to store runoff.
The third scenario also relates to the potential for saturation to occur above the hardpan, allowing subsurface saturated flow to occur (figure 7c). The changes in pattern properties with distance to the stream suggest that water availability increases downslope: this requires subsurface storage, if not flow. A relationship between the pattern formation and subsurface flow could explain the local deviations of pattern orientation from slope orientation near the streamlines since the surface orientation might not correspond to the local water table flow direction.
All three scenarios would tend to strengthen the pattern-forming feedbacks on the shallow soils. They could be investigated using subsurface soil moisture sensors to observe shallow soil saturation; water isotope tracers to determine the water sources used by plants; observations of calcium ion concentrations in runoff which would provide an indicator of water contact with the petrocalcic horizon; and root excavations to compare morphologies in sites with different depths to the impeding layer. These studies could be valuable to provide more information about the hydrological role of petrocalcic horizons, which have received relatively little attention given their ubiquity in desert environments [73].
5. Conclusion
Large-scale analyses of variation in pattern morphology have provided broad confirmation of many theoretical predictions about vegetation patterning and its variation along climatic gradients [30], which are consistent with our observations at local scales. Our analysis shows that while pattern properties mostly vary on scales of 600–800 m, they can change much more rapidly around boundaries in topography or soil characteristics. By analysing the pattern on these fine scales, we identified deviations between hillslope and pattern orientation, soil-controlled changes in pattern wavelength, coherence and vegetation cover, and at least one likely example of the boundary condition effects predicted by McGrath et al. [17]. Two observations, echoed at multiple other sites, are not well explained by current theory: the observation of decreasing pattern wavelength with steepness, and the soil controls on vegetation pattern length scale and coherence.
The large time-scale separation between individual storm events and the time scales on which desert vegetation distributions change mean that an ongoing dialogue between empirical and theoretical studies is critical for understanding the dynamics of these ecosystems. Local information about pattern qualities, when combined with high-resolution information about the pattern substrate, is evidently a useful additional tool for analysis.
Despite the advances in understanding vegetation patterns in the past 10 years, theoretical models still require the use of effective parameters to describe feedbacks, soil–plant–water interactions and the resulting landscape fluxes. Linking theory and observation to make quantitative predictions, therefore, remains an outstanding challenge. Addressing this challenge requires improved observations of within-storm hydrologic processes and plant water use: observations that ecohydrologists are increasingly equipped to make due to developments in distributed sensing systems and tracer technologies. Computational tools for assessing three-dimensional soil moisture dynamics, land–atmosphere interactions and vegetation spread are also improving. By coupling these tools with detailed field measurements, there is potential to develop a detailed theoretical framework that can address the consequences of changing soil depth, root orientation, runoff generation mechanisms and subsurface flow processes on the overall dynamics and resilience of vegetation patterns.
Funding statements
S.E.T. and G.P. acknowledge support from NSF-EAR-1013339.
