Dynamics and resilience of vegetation bands in the Horn of Africa

Bands of vegetation alternating periodically with bare soil have been observed in many dryland environments since their discovery in the Horn of Africa in the 1950s. Mathematical modeling efforts over the past two decades have sought to account for these bands via a self-organizing interaction between vegetation and water resources. Many model predictions for vegetation band response to climatic and human pressures, such as an increase in band wavelength, have yet to be observed, perhaps due to a lack of long-term observation. In this study we bring photographs taken in the 1950s over the Horn of Africa together with present-day satellite imagery to assess vegetation change over six decades. For many areas, changes in the vegetation bands are modest, with individual bands remaining identifiable between the images. Other areas have experienced substantial increases in land development since the 1950s, with many bands disappearing from the landscape. While wavelengths have remained approximately unchanged in all study areas, bands in the human-impacted areas have counterintuitively widened in the direction of slope, suggesting that band widths respond to human impacts on time scales relevant for remote monitoring. Our observations underscore the importance of human impacts to the resilience of dryland vegetation, and indicate novel future directions for both empirical and theoretical investigation.

Dryland environments cover over 40% of Earth's terrestrial area and are home to more than a third of the global human population [1]. Vegetation in dryland environments is often patchy in response to resource limitation [2]. Water redistribution can cause patchiness in the form of periodic patterning at length scales of tens to hundreds of meters [3,4]. Notable instances of such patterning are bands of vegetation separated by stretches of bare ground aligned perpendicularly to gradual slopes. Such bands have been reported at multiple sites in Africa, North America, and Australia [5], and comprise the majority of field observations of periodic dryland vegetation patterning [6]. Vegetation bands often occur on land used for pastoral farming in Africa [7] and can serve as crucial buffers against land degradation [8,9].
Due to the large spatial scales of vegetation banding, studies of the phenomenon have primarily relied upon remotely-sensed data and mathematical modeling. Aerial survey imagery taken in the 1940s and 50s over then British Somaliland was critical to the discovery [10] and initial characterization [11][12][13][14] of the phenomenon. These early studies and subsequent field work [15,16] identified conditions under which bands occur and laid out ecohydrological mechanisms that explain their formation. Modeling efforts since the late 1990s have focused on mathematizing these mechanisms through partial differential equation approaches [17][18][19][20]. Such models are able to qualitatively reproduce the appearance of the bands, as well as the slow uphill migration of bands estimated in early fieldwork [12].
Many model investigations have focused on the dynamics and resilience of vegetation in response to changes in aridity or rainfall conditions [4]. As aridity increases, banded vegetation is predicted to increase in wavelength [21][22][23][24][25], to undergo abrupt shifts in migration speed [26,27], and to become susceptible to break-up that results in a spotted state [28]. Rainfall fluctuations are predicted to affect migration speed and the widths of bands [24,26]. Grazing pressure, commonly viewed interchangeably with human impact, has been modeled as an increase to the plant mortality rate [22,29], and by disturbances * msilber@uchicago.edu to the state of vegetation. Such disturbances can cause an increase in band wavelength [30] or band break-up [31]. In many of these investigations, the simulated dynamics of environmental change and the resulting ecological shifts play out over the equivalent of decades to centuries.
Evidence for predictions related to the dynamics and resilience of vegetation bands has been limited by an irregular record of remotely-sensed observations from the pre-Landsatera. Wu et al. used aerial survey photography of tiger bush banding in Niger over three decades to observe vegetation fragmentation [32], a phenomenon that appears related to decreases in band width and break-up. Valentin & d'Herbès used similar imagery over Niger to observe that the widths of bands fluctuated in response to rainfall history over four decades [8]. Deblauwe et al. used declassified reconnaissance imagery to measure appreciable band migration over four decades in Texas, USA, and the Haud region of Somalia. They also observed band width fluctuation in response to rainfall in Texas [33]. Additional studies of dryland vegetation patterning on flat terrain found that the vegetation transitioned between qualitatively distinct states during a period of significant drought [34], and that human pressure hastened this type of transition [35].
In this study, we investigate the dynamics and resilience of vegetation bands in two distinct regions of Somalia ( Figure 1). We collected and precisely aligned digital scans of British Royal Air Force (R.A.F.) aerial survey photography taken in 1952 and multispectral satellite imagery taken in the last decade. Among our study areas are locations that remain relatively pristine, and also locations that experienced a dramatic increase in human pressure. Both regions have experienced multi-year fluctuations in rainfall and a warming trend over the last half-century. The aim of this study is to assess the nature of measurable change in vegetation banding on a multi-decadal timescale via a systematic visual comparison of imagery. We also identify band properties for which changes can be readily observed in remotely-sensed imagery. Our approach relies on comparative Fourier analysis to identify band wavelength and migration and transect-based measurements of band width.

Human pressure and band degradation
We assessed changes in human pressure and vegetation banding over time through a systematic visual comparison of R.A.F. photography and recent satellite imagery (Methods, Supplementary Information). We found that band degradation ranged from partial to complete band loss in areas with the sharpest increases in human pressure, while bands in the other areas remain largely unchanged ( Figure 2). Roads and dirt tracks can be visually identified in both the aerial photographs and the satellite imagery, and their presence and qualitative appearance served as our primary proxy for inferring the extent of human pressure. Vegetation in both the aerial photos and satellite imagery contrasts sharply with the light background of bare soil, and bands are clearly identifiable. We defined degradation in this context as either the breakdown in regularity or the disappearance of banding.
Substantial road and track development occurred in much of the Sool Plateau, with most areas (SP1-SP4) transitioning from having either no roads or faint roads in 1952 to having roads or tracks that densely cover the landscape in the modern images ( Figure 2a). The settlement Dhahar was founded within SP4 after the 1952 photograph, and now supports a population of approximately 13,000 ( Figure 2c). We observed much less road and track development in SP5 and in the Haud (HD1-HD4) (Figure 2b). At many sites within the Haud, human-made structures visible in the 1952 images appear to persist into the current decade, suggesting no major change in land use over the intervening time (Supplementary Figure S5a,b).
Band degradation is prevalent in the human-impacted areas SP1-SP4 (Figure 2a). Bands have disappeared entirely from the landscape in large parts of SP3 and SP4. Only part of the band loss in these areas appears directly related to clearing for land development, since loss also occurs in areas without humanmade structures. Where faint remnants of bands are visible in SP4, degradation appears to have occurred without a visible change in band wavelength (Supplementary Figure S5c). In SP2-SP4, dense tracks often appear between bands (Figure 2d). Frequently we observed vegetation growing within roads and tracks, which suggests that these structures likely disrupt the flow of water on the landscape. In SP5 and HD1-HD4, individual bands often remain identifiable after six decades based on visible details of their morphology, and we observed no substantial band degradation.

Band wavelength and migration
Models predict that band wavelength should increase in response to sufficient increases in aridity or environmental pressure [21-23, 25]. We quantified band wavelength change in all study areas (  Table S2).
We also found that bands in all study areas migrate uphill to a detectable but modest degree, typically less than a quarter of the wavelength over six decades (Table 1, Supplementary  Table S2). The migration rates we measure are comparable to the values estimated by Hemming in 1965 [15] and measured recently by Deblauwe et al. [33] in nearby areas of the Haud. We found that band wavelengths and migration rates are significantly positively correlated (r = 0.57-0.86, p < 0.01) at all areas except for SP4, where substantial band degradation likely causes error in the measurement. This correlation also aligns with the findings of Deblauwe et al. [33] for an area in the Haud, and is consistent with the model prediction of an increasing relationship between migration rate and wavelength [22].
Band wavelength is predicted by models to vary with local slope, though the nature of this relationship can be parameter and history-dependent [37, 38]. We found a significant correlation between wavelength and slope at SP3 only (r = −0.34, Table S2). The sign of this correlation agrees with empirical findings in other studies [33,34,36]. Remarkably, we found significant negative correlations between migration and slope at multiple areas in the Sool Plateau and the Haud (Table 1, Supplementary Table S2). These findings contradict the intuitive model prediction of an increasing relationship between migration rate and slope [39].

Band widening in human-impacted areas
During the course of visual comparison, we noted that bands in some areas appear to widen over time in the direction of local slope. We quantified changes in band width in all study areas using automated transect measurements of individual bands pressure. The qualitative state of road and track cover, settlements, and vegetation banding was assessed visually in 1 km 2 boxes. Transitions between states observed in aerial photography and recent satellite imagery are shown here for all study areas. (a) SP1-SP4 in the Sool Plateau show a high degree of road and track development, and a moderate to high degree of band loss. A large settlement (Dhahar) developed within SP4, and is indicated with a black border. SP5 saw little increase in road cover, and no band loss or degradation was observed. (b) Areas in the Haud (HD1-HD4) show only a small increase of road and track development, and no substantial band loss or degradation is observed. (c) An example of band loss and degradation due to land development in SP4 (9. stant at the other study areas ( Figure 3).
We computed the ratio of band widths measured in recent imagery to the widths in 1952. The median ratio exceeded 1.2 at the human-impacted areas SP1-SP4 ( Figure 3a). The most substantial widening occurs at SP2, where the median ratio is 1.8. We measured band widths at SP2 using additional images taken in 1967, 2004, 2006, 2011, and 2013. We found that widths did not change between 1952 and 1967, and then nearly doubled  Figure S8). Since recent images were taken in a variety of seasonal and rainfall history conditions (Supplementary Figure S1), we conclude that the widening observed in SP1-SP4 is not an artifact of seasonality.

Discussion
In many areas, remarkably little about the vegetation bands has changed. We found no evidence of systematic changes in wavelength despite apparent increases in environmental pressure in the region over the last 65 years (Supplementary Section S1.1). Moreover, individual bands remained largely identifiable between images over time, which allowed us to observe modest uphill migration. Our findings suggest that large-wavelength (∼150 m) bands should migrate appreciably (∼5 m) over the span of a decade. Given sufficiently accurate image alignment, migration of this magnitude should be detectable in current high-resolution satellite imagery, and may in principle be brought to bear on relevant model predictions [22,27]. Additionally our findings suggest a counterintuitive relationship between migration and slope, one that may be stronger than the relationship between wavelength and slope. Further theoretical and empirical investigation of this relationship may shed light on the connection between topography, hydrology, and the emergent scales of vegetation banding.
In other areas, the most profound changes in vegetation banding are associated with increases in human pressure. This mirrors the findings of previous empirical studies [32,34,35], and underscores the importance of human impacts to the resilience of dryland vegetation. Apart from land clearing for development, we observed widespread degradation in areas with few human-made structures, suggesting more subtle but important forms of impact. Roads and dirt tracks now densely cover the landscapes of many of the study areas, and have likely affected the flow and availability of water to the vegetation bands [41]. Biomass harvesting, here in the form of Acacia cutting for charcoal production [42], has no doubt also played a role in the degradation we observed. Both changes in road cover and vegetation harvesting can potentially be monitored remotely, and such monitoring efforts should be well-informed by theoretical investigations that account for these and other forms of human impact.
Surprisingly, we observed a widespread and persistent increase in band width at the human-impacted areas SP1-SP4. In previous studies of Niger and Texas, USA, researchers observed sizable fluctuations in band width on sub-decadal time scales in response to rainfall variation, with favorable rainfall conditions resulting in wider bands [8,33]. In the Sool Plateau, we did not observe appreciable band width fluctuations in response to multi-year rainfall variation, nor did we find evidence that rainfall conditions have become more favorable in the region over time (Supplementary Section S1.1). The fact that widening is localized to Sool Plateau areas SP1-SP4, that it occurs to different degrees in each of these areas, and that it is not observed at nearby area SP5 strongly suggests that non-climatic factors have driven the apparent changes. We simulated the model by Klausmeier [18,43] to examine how changes in band width might be achieved through factors other than rainfall (Methods, Supplementary Information).
Acacia cutting for charcoal production has been prevalent in the Sool Plateau since at least the 1980s [42], and has likely caused a decrease in woody biomass within the vegetation bands in many areas. A shift in composition from woody to grass biomass plausibly increases the effective transpiration rate, biomass yield per unit water, dispersal rate, and mortality rate of the vegetation. We found that individually increasing transpiration, yield, and dispersal rates resulted in band widening, while increasing mortality rate had the opposite effect ( Figure 4). Of these changes, only increases in dispersal rate result in both band widening and a decrease in peak biomass (Figure 4a), a result consistent with our observation of a modest inverse relationship between band width and vegetation index (Supplementary Section S5.3, Table  S3). We found that band widening alongside diminished peak biomass can be achieved through a simultaneous increase in transpiration rate, biomass yield per unit water, dispersal rate, and mortality rate parameters (Supplementary Figure S10), and conclude that a shift in species composition is one viable explanation for the changes we have observed in the Sool Plateau.
Band wavelength has been a prime focus of theoretical and empirical investigations of vegetation pattern resilience. The dependence of wavelength on model parameters can often be studied analytically, and it is conceived that wavelengths are easily measured in remotely-sensed imagery. However, as we and others [36] have found, the spacing between bands is often quite irregular, making imprecise both the notion of wavelength and its measurement. Moreover, because a range of band wavelengths may be stable over a range of environmental conditions, wavelength changes in model scenarios of environmental change are history-dependent and discontinuous [22][23][24][25]. These shortcomings make wavelength changes a poor signal of ecosystem regime shift, and explain the dearth of evidence for such changes in this and other studies.
We have provided new evidence that vegetation band widths change on observable time scales, and we argue that they represent an underutilized window into the response of dryland a b +10% biomass yield per unit water initial +10% mortality +10% transpiration +50% dispersal vegetation to climatic and human pressure. Our model investigation suggests that band widths change continuously in response to parameter variation, and that these changes should be observable on sub-decadal time scales (Figure 4b). Band width is straightforward to measure in remotely-sensed imagery, and comparisons over time need not depend sensitively on the quality of image alignment. Future theoretical investigation of this pattern property will be important to establishing its utility to dryland vegetation monitoring.

Regional information
We studied areas within the Sool Plateau and Haud pastoral regions of Somalia. Both regions are generally characterized by an arid climate (aridity index = 0.04-0.1) [44]. Due to a lack of continuous rainfall station monitoring in and around our regions of study, we assessed the historical regional climate using 20th Century Reanalysis [45] and the CPC/Famine Early Warning System Dekadal Estimates datasets. Mean annual rainfall in both regions ranges between 100-300 mm. We found no evidence that rainfall conditions have improved in either region in recent decades, and we identified a warming trend in average yearly temperature of 1-2 • C over the last half-century. Regional soils are claylike and prone to crust formation, resulting in low permeability and surface water runoff following high-intensity rainfall [15,42]. Hemming found that soils are wetter beneath bands in the Haud, indicating greater soil permeability in vegetated areas [15]. Vegetation bands in both regions are dominated by Andropogon kelleri grasses [10,15]. Bands also contain a mix of trees and shrubs, most notably Acacia bussei. In recent decades, Acacia bussei has diminished in abundance in the Sool Plateau due to cutting for charcoal production [42]. Disruption of traditional grazing patterns has resulted in overgrazing in many areas of the Sool Plateau, including Dhahar (SP4) [42].

Data
We studied approximately 260 km 2 of imagery within the Sool Plateau and 200 km 2 of imagery within the Haud. Study areas were chosen based on a combination of factors; in particular, we wished to include areas with different development and degradation outcomes, areas with recorded soil and floristic information based on field studies, areas in geographically distinct regions, and areas featuring well-defined banding. Study area boundaries are defined by our choice of British Royal Air Force (R.A.F.) aerial survey photography, which comprise our earliest image datasets. Aerial survey photographs were taken in 1951-52 over broad areas of British Somaliland, and specific photographs were scanned on request by the Bodleian Library at the University of Oxford. We also studied declassified reconnaissance satellite imagery taken in 1967, and Digital-Globe imagery for dates spanning 2004-2016. Resolution of imagery used in this study ranges between 1.4-2.4 m/pixel. Satellite images taken between 2004-2016 containing red and near-infrared channels were used to compute a Soil-adjusted Vegetation Index [46]. We manually georeferenced R.A.F. scans and the 1967 reconnaissance image using visually identified control points. We estimated alignment error to be approximately 1-2 pixels.
We estimated gradient for our study areas using the Shuttle Radar Topography Mission Global 1 arc second elevation dataset [47]. Because of the noise characteristics of the dataset and the low relief of our study areas, we used a second-order finite difference operator with noise-suppressing properties to estimate gradient and slope [48].

Visual comparison
We assessed changes over time through a systematic visual comparison of imagery. We developed a graphical user interface in MATLAB for comparing images. For each area, we split both R.A.F. scans and recent imagery into 1 km × 1 km boxes, and evaluated qualitative features within these boxes. We evaluated the extent of roads and dirt tracks, which served as our primary proxy for human pressure. We also evaluated the extent of banding and defined degradation in this context as the breakdown in regularity or disappearance of banding over time.

Fourier analysis
We measured changes in band wavelength and band migration by modifying the Fourier window method by Penny et al. [36]. The method measures band wavelength and orientation in a sliding window using a 2D FFT, and computes a uniqueness metric based on the unimodality of the power spectrum. We modified the method to also estimate band migration by calculating the phase difference of the common dominant Fourier modes between corresponding windows at two different time points. In all analyses, we discard data points which correspond to sites without banding using a manually-drawn mask. We additionally discard data points with uniqueness values below a threshold.

Band width measurement
We quantified band widening in the direction of slope using by gathering image intensity profiles along transects through the bands in the direction of slope. We used the same transects for multiple images at the same study area. We then fit a simple plateau function to each intensity profile to extract band width. Band width was measured along multiple parallel transects, and data points with high variance in measured widths were discarded.

Model simulations
We simulated a 1D version of the model by Klausmeier [18,43] to estimate the sensitivity of band width and peak band biomass to parameter changes. We obtained the initial parameter set from the values and ranges given in [18]. Parameters stated in [18] to differ between grasses and trees are set at intermediate values so that the spatial scale of banding resembles those in our regions of study. The time scale of migration was similarly tuned using the downhill water flow rate parameter. We simulated the model using the ETDRK4 pseudospectral scheme [49].

Statistics
We assessed the significance of correlations between quantities computed from Fourier and band width analyses using a two-tailed paired t-test corrected for spatially autocorrelated data [40], implemented in SpatialPack for R [50].

Data availability
The imagery data that support the findings of this study will be made available on a publicly-accessible repository.

Supplementary information Dynamics and resilience of vegetation bands in the Horn of Africa
Karna Gowda, Sarah Iams, and Mary Silber 1 Contents S1 Regional information 1 S1

S1 Regional information
We studied imagery in areas located within the Sool Plateau and Haud pastoral regions of Somalia. Sool Plateau study areas are located approximately 50 km west of Gardo, and Haud areas are located approximately 40 km south of Las Anod. Areas were chosen for this study based on a combination of factors; in particular, we wished to include areas with different development and degradation outcomes, areas with recorded soil and floristic information based on field studies, areas in geographically distinct regions, and areas featuring well-defined banding. Detailed information about areas and imagery is given in Section S2.1 and Table S1.

S1.1 Climate
Both Sool Plateau and Haud pastoral regions are characterized by an arid climate (aridity index = 0.04-0.1) [1]. Rainfall in Somalia is bimodally distributed between the Gu season, spanning Apr. Due to a lack of continuous rainfall station monitoring in and around our regions of study, we assessed the historical regional climate using climate reanalysis and remotely-sensed rainfall estimation datasets. The 20th Century Reanalysis (V2c) dataset assimilates surface pressure observations, sea-surface temperature, and sea ice extent into a global climate model to obtain a reconstruction of Earth's climate spanning 1871-2011 [2]. The V2c dataset is available at 6-hour temporal and 2 • spatial resolution, and uncertainty estimates can be derived from 56 replicate model simulations. The CPC/Famine Early Warning System Dekadal Estimates (RFEv2) dataset uses satellite microwave sensing and ground station observations to estimate total rainfall over the African continent for dates spanning 2000 to present. RFEv2 data is available at daily intervals and 0.25 • spatial resolution.
To assess the rainfall conditions surrounding our imagery datasets, we obtained annual total rainfall estimates from the V2c and RFEv2 datasets ( Figure S1a-d). In the absence of ground confirmation, we exercise caution in interpreting the V2c estimates for the 1940s-60s, and conclude that there is no evidence that rainfall conditions have improved in either region in recent decades. We speculate that rainfall conditions surrounding the 1952 and 1967 datasets were quite favorable. We also speculate that conditions have either declined or reverted to a regional mean in recent decades.
We assessed rainfall conditions for the recent imagery in greater detail by calculating seasonal rainfall totals from the RFEv2 dataset ( Figure S1c-d). In the Sool Plateau, the images used in this study were taken in a variety of seasons and rainfall history conditions. The 2004 image was taken shortly after the return of rains that followed a very severe multi-year drought. The 2006, 2011, and 2013 images were taken amid more typical rainfall conditions. The 2016 image was taken during a period of drought, which is ongoing at the time of writing. In the Haud, the images used in this study were taken in years (2012 and 2016) with robust rainfall during the wet seasons. We examined regional temperature history using surface temperature estimates from the V2c dataset. We computed the average yearly temperature, defined as a yearly average over the daily midpoint between minimum and maximum temperatures ( Figure S1e). We identified a distinct linear warming trend between 1960 to the present of 1-2 • C.

S1.2 Vegetation
Field investigations in the 1950s-60s found vegetation bands in the Sool Plateau and the Haud to be dominated by Andropogon kelleri grasses [3,4]. The bands also usually contain a mix of trees and shrubs. Long-lived Acacia bussei trees often populate the bands. In many of our study areas, bands occur on relative low-ground (e.g., in channels), while a more uniform cover occurs on the surrounding higher ground. Greenwood attributes this difference to the clay content of the soil: higher-ground areas have a lower clay content, and thus can support a greater density of vegetation [5].
In recent decades, Acacia bussei has diminished in abundance in the Sool Plateau due to cutting for charcoal production, as it is considered the most lucrative species for this purpose [6]. Disruption of traditional grazing patterns has resulted in overgrazing in many areas of the Sool Plateau, including Dhahar (SP4).

S1.3 Soil
Soils studied around areas with vegetation bands in the Haud [4] and Sool Plateau [6] are claylike and finely textured. In both regions soils appear prone to crust formation and soil pore plugging, resulting in low permeability and surface water runoff following high-intensity rainfall. Hemming found that soils are wetter beneath bands in the Haud, indicating greater soil permeability in vegetated areas [4]. Soils in the Sool Plateau thinly cover a limestone bedrock, which is exposed in some parts of our study areas.

S2.1 Imagery
We studied approximately 260 km 2 of imagery in areas of the Sool Plateau and 200 km 2 of imagery in areas of the Haud. Study area boundaries are defined by our choice of British Royal Air Force (R.A.F.) aerial survey photography, which comprise our earliest image datasets. Aerial survey photographs were taken in 1951-52 over broad areas of British Somaliland, and are archived at the Bodleian Library at the University of Oxford. The aerial photographs used in this study was scanned on request by the Bodleian Library using British Ordnance Survey maps to identify images.
The coordinates of study areas and additional information about imagery used in this study are given in Table S1. R.A.F. images were scanned at a nominal resolution of 1.4-2.5 m/pixel. We obtained more recent imagery through the USGS and DigitalGlobe Foundation. We purchased declassified reconnaissance satellite imagery 1 taken in 1967 from the USGS Earth Explorer site. We downloaded freely available Orbview-3 imagery taken in 2005 from the USGS Earth Explorer site. We were granted Quickbird-2, WorldView-1, and WorldView-2 imagery for dates spanning 2004-2016 by the DigitalGlobe Foundation.
We manually georeferenced aerial survey photograph scans in ArcMap 10.3 against the ArcGIS World Imagery layer using the WGS84 Web Mercator coordinate system (EPSG:3857). Because vegetation bands migrate over time, we could not match scans with geospatial coordinates using the appearance of the bands themselves. Instead we relied upon apparent geological features, such as limestone outcrops, and geometrically distinct clusters of individual trees or shrubs that persisted over time. Aerial survey photographs were matched using no fewer than 10 control points per image, and were aligned by fitting a projective transformation. Control points were stored in a tab-delimited file. A projective transformation is overdetermined for greater than 4 control points, and the root mean squared error (RMSE) of the transformed control points served as our estimate of alignment error.
To estimate the effect on RMSE of adding additional control points, we used a resampling procedure that calculates the alignment RMSE for different subsets of the control points. For an image that was aligned using n control points, we computed the RMSE for permutations of 5 ≤ k ≤ n control points ( Figure S2). The average RMSE values over the permutations were then computed for each value of k. In this procedure, if the total number of such permutations ( n k ) exceeded 10 3 , a random sampling of 10 3 distinct permutations were used. Otherwise, all permutations were used. We then fit the resulting curves by the function ak/(1 + bk), wherek = k − 5, to extrapolate the saturating value of the average RMSE curve. In all cases the saturating RMSE value was comparable to the resolution of the imagery, suggesting an alignment error on the order of 1-2 pixels. A reconnaissance satellite image was also manually georeferenced in ArcMap 10.3 using a third-order polynomial transformation with 18 control points. The image covers a much broader area than the aerial photographs, and due to distortions arising from the imaging methodology a projective transformation did not produce a suitable fit 2 . RMSE of this alignment is 0.94 m, which is on the order one pixel. DigitalGlobe imagery was pre-georeferenced, and was manually shifted to align more precisely with the ArcGIS World Imagery layer.
Most recent satellite imagery used in this study contain data sensed at different frequency channels. The red, green, and blue channels were used for visualization, and the red and near infrared channels were used for computing the Soil-adjusted Vegetation Index (SAVI) [8], an index of photosynthetic activity: N IR is the near-infrared intensity value, and R is the red intensity value. The parameter L is used adjust for exposed soil surface in low-vegetation cover scenarios, and is often used in place of Normalized Difference Vegetation Index (NDVI) in dryland vegetation inference. SAVI is equivalent to NDVI for L = 0. For all analyses, we used a conventional value of L = 0.5.

S2.2 Elevation
We used NASA Shuttle Radar Topography Mission Global 1 arc second (SRTMGL1) elevation data for our upslope migration assessment and comparison of pattern properties with slope. Datasets were obtained from the USGS website 3 . Datasets are packaged in 1 • latitude × 1 • longitude tiles, and were loaded, georeferenced, projected onto the WGS84 Web Mercator coordinate system (EPSG:3857) in MATLAB 2016b. This allowed for elevation data to be matched with imagery. Though the SRTMGL1 dataset has a nominal resolution of 1 arcsecond (∼30 m/pixel near the equator), the true resolution is closer to 45-60 m/pixel due to the manner in which data was collected [9]. The data also contains speckle noise which is autocorrelated at a length of 1-2 pixels, and also random error, both of which together result in average vertical error of approximately 4 m in areas like the Sahara Desert. To eliminate autocorrelated errors, we subsampled the data to 3 arcsecond (∼90 m/pixel) resolution.
Noise in the SRTM data presents a challenge to gradient estimation in areas of low relief, such as our regions of study, where in banded areas vertical change can be as little as 1 m per 500 m of horizontal change.
To compute gradient fields, we used a second-order accuracy finite difference stencil with noise suppressing properties [10]. As an example, a 5 × 3 noise suppressing gradient operator as defined in [10] is where h is the discretization step size. Convolving this operator with the data array produces an approximation of the partial derivative in one direction. Operators with noise suppressing properties discussed in [10] can be computed for arbitrarily large stencil size. Using a finite difference operator allows for a straightforward propagation of i.i.d. normal errors in the elevation data through the calculation of gradient and slope. We estimated the magnitude of errors in the elevation data by computing the standard deviation of residuals from a median subtraction: Slope is obtained from the magnitude of the gradient vector, and to leading order the gradient error value is also equal to the error propagated to the slope calculation.
We tested the sensitivity of slope calculations to varying stencil size s (which yields a (2s + 1) × (2s − 1) operator). We note again that the truncation error of the finite difference operator is second-order for any s. Intuitively too small a stencil size will have high measurement error, and too large a stencil size will result in oversmoothing. In Figure S3 we show the 25th, 50th, and 75th percentiles of the slope values at each study area computed over an interval of s, and indicate one standard deviation of propagated error around these values. We conclude that slope values are not sensitive to stencil size when s ≥ 15, and we use 15 (which gives an operator of size 31 × 29) for all slope and gradient calculations. We confirmed by visual inspection that this stencil size produces smooth gradient fields that match hydrological features visible in the imagery (e.g., hills and channels).

S3.1 Protocol
We assessed changes over time at study areas via a systematic visual comparison of imagery. Roads can be visually identified in both the aerial photographs and the satellite imagery, and their presence and qualitative appearance served as our primary proxy for inferring the extent of human pressure. Vegetation in both the aerial photos and satellite imagery contrasts sharply with the light background of bare soil, and bands are clearly identifiable. Degradation was inferred through either the breakdown in regularity or disappearance of banding.
We developed a graphical user interface (GUI) in MATLAB 2016b for visually comparing images (Figure S4). The GUI allows the user to select two imagery datasets for comparison, a georeferenced R.A.F. photograph and a more recent image. Images used for visual comparison are indicated in Table S1. The recent image is projected onto the intrinsic coordinate system of the R.A.F. photograph, so that the data can be cleanly divided into non-overlapping 1 km × 1 km windows. The GUI simultaneously displays corresponding 1 km 2 windows of the R.A.F. photo and more recent imagery. Additionally the GUI displays a false color overlay of the two images, which was used to assess whether migration occurred. The GUI plots the local slope direction vector (computed as described in Section S2.2) on top of the overlay, which allows the user to visually assess whether the migration is in the upslope.
For each image window, the user is prompted to enter whether regular banding is present and whether a dense settlement is present (more than 5 structures in close proximity) using checkboxes. The user can select the extent of apparent road cover via a dropdown menu. Additionally, the user can check boxes to indicate whether band widening in the slope direction is apparent, and whether it appears that the same roads or settlements are present in both images. The user can enter comments for each image and the overlay. If the recent image contained red and near infrared channels, the Soil-adjusted Vegetation Index (SAVI) can be displayed in place of the RGB image (see discussion of SAVI in Section S2.1). The GUI selections are automatically saved to a MATLAB .mat file. When the user has finished assessing the window, the user can then navigate to different windows in the dataset using buttons.
The checkboxes are ternary; with banding, for example, a fully-checked state is taken to mean distinct banding, a half-checked state is taken to mean indistinct banding, and an unchecked state is taken to mean no banding. The road cover dropdown can take on one of four states: "no roads," "faint road(s)," "clear road(s)," and "clear, dense road(s)." The last state, "clear, dense road(s)," is taken to mean that a large number of well-incised, clearly visible roads cover a large portion of the window.

S3.2 Highlighted examples
In many sites of the Haud study areas, we observed that human-made structures appeared to persist from 1952 to the present. In Figure S5a,b, we show two examples of such structures.
In addition, we observed that in some areas of SP4, bands appeared to degrade without apparent change in wavelength. We show an example in Figure S5c. We verified that the pictured bands have significantly lower vegetation index (SAVI) values than nearby bands in the study area, and are plausibly degraded.

S4.1 Protocol
We quantitatively assessed changes in band wavelength and band migration at study areas using a modification of the Fourier window method by Penny et al. [11]. Penny et al. developed the method to compute spatial maps of local wavelength and orientation from imagery over banded areas in Fort Stockton, Texas, USA. In a manner analogous to a short-time Fourier transform, the method measures wavelength and orientation in a sliding window using a 2D FFT. Vegetation banding typically contains sufficient irregularity to complicate the inference of dominant wavelength and orientation from a 2D power spectrum. The Fourier window method addresses this issue by binning power, radially for estimating wavelength and angularly for orientation, and by computing a weighted average among the contiguous bins with largest power. The method computes a uniqueness metric for both wavelength and orientation based on the distance between the maximal peak and the nearest peak with 75% of the maximal power, if present. The metric equals one if the maximal peak is the only powerful peak present, and approaches 0 as distance to the nearest powerful peak increases. In order to exclude short-wavelength noise and long wavelengths which are under-sampled for the given window size, the bins are only computed for a specified minimum and maximum wavelength interval. As presented in [11], the pattern irregularity issue is also addressed by averaging measurements over overlapping windows. We do not perform the latter step for the analyses in this study.
Penny et al. provide MATLAB code for their method, which we modify for our analyses. We modify the main routine to take as input two images that have been resized to the same dimensions. Computations are then performed on square windows. To reduce aperiodicity effects, we apply a 2D Hamming filter (a bellshaped function that decays to zero away from the center) to each window. We note that this filter also has the effect of giving more weight to the central area of the image, focusing analysis on this area. Wavelengths and orientations are then computed using Penny et al.'s routine. Our most substantial modification is to also estimate band migration between the two windows by calculating the phase difference of the common dominant Fourier modes. This procedure assumes that migration resembles a spatial translation of the dominant Fourier modes in the direction of their wave vectors. It also requires that band migration is no more than half a wavelength, which we had previously verified through visual inspection.
To measure migration, we computed the cross FFT and cross power spectrum of the window pair: imfft1 = fft2 ( im1 ); %FFT for window 1 fftpower1 = imfft1 .* conj ( imfft1 ); %calculate power spectrum imfft2 = fft2 ( im2 ); %FFT for window 2 5 fftpower2 = imfft2 .* conj ( imfft2 ); %calculate power spectrum imfft12 = imfft1 .* conj ( imfft2 ); %cross FFT fftpower12 = abs ( imfft12 ); %calculate cross power spectrum We compute the migration distances by multiplying the phase of each mode in the cross FFT by the wavelength of that mode. We bin the cross power spectrum both radially and angularly, and identify the bin with maximal power. Mirroring the Penny et al. binning method, we identify the adjacent bins that have at least 75% of the maximal power. We then calculate a weighted average of migration distance across these largest bins, where the power of the bins serve as the weights. This is taken to be the estimated migration distance for the window pair. We also compute a migration uniqueness metric in a manner similar Penny et al.'s wavelength uniqueness metric. We applied this methodology to all study areas using the image pairs indicated in Table S1. In order to perform the windowing on a rectangularly-oriented dataset, we transformed the recent imagery onto the intrinsic coordinate system of the aerial photograph. We downscaled image pairs to a resolution of approximately 2.5 m/pixel to reduce computation time. We then applied two layers of preprocessing to emphasize the vegetation bands and de-emphasize other features in the imagery: we subtracted a coarsely gaussian-blurred version of the image to eliminate large scale variations in pixel intensity (such as darkening near the borders of the aerial photographs), and we applied a manually-tuned threshold to create a binary image of the vegetation bands.
For each binary image pair, we computed wavelength, orientation, and migration maps for the three square window sizes: 384 pixels (∼1 km), 512 pixels (∼1.3 km), and 768 pixels (∼2 km). After applying a Hamming filter, about 4-8 vegetation bands can be sampled in the central area of a 512 pixel window (3-5 bands for a 384 pixel window, or 6-10 for a 768 pixel window). For the binning procedure, we set the minimum wavelength to 10 pixels (∼25 m), and the maximum wavelength to one-fourth of the window size (∼240 m for a 384 pixel window, ∼320 m for a 512 pixel window, and ∼480 m for a 768 pixel window). In order to balance computation time and even-sampling of the data, we set the step length of the sliding window to be one-fourth of the window size, resulting in adjacent windows that overlap in 75% of their area.
After wavelength, orientation, and migration maps were computed for an imagery pair, we transformed the measurements from units of pixels to units of meters. We then manually drew a mask on the imagery and applied it to the measurements in order to exclude measurements from areas without vegetation bands. Additionally we excluded measurements with wavelength or migration uniqueness metrics smaller than 0.75. We found that in some areas, a window size of 384 pixels was too small to detect the largest wavelengths. The results of 512 and 768 pixel windows did not differ strongly, so we used the 512 pixel window computation for the results reported in the main text and Table S2.

S4.2 Wavelength change
Given spatial maps of wavelength for a pair of images, we computed change maps where elements are given by where W i,j 1 is the wavelength in the first image at position (i, j), and W i,j 2 is the wavelength in the second image at (i, j) ( Figure S6). Typical ranges of the computed changes at each study area are given in Table S2. Note that for computing change maps, we have manually masked out areas with no banding, and we have also excluded measurements with a wavelength uniqueness metric smaller than 0.75. We made the latter choice to reduce the incidence of falsely detected changes between the maps, reasoning that measurements in areas with multiple dominant band wavelengths are error prone. Typical change ranges between 0-10% for all study areas except SP3 and SP4, where change ranges between 0-20%.
Our previous visual inspection suggested that there were no obvious systematic changes in wavelength at any study area, except perhaps for those associated with isolated instances of band loss in human-impacted areas. We visually reinspected areas where measured change was greater than 25% in magnitude. In some cases, it appears that these detected changes occur due to the loss of an individual band, often near evidence of human activity ( Figure S7). In Figure S7b-e, significant road cover appears in the interband areas of the recent imagery, and are likely related to the loss of bands in these areas. In most cases, however, we saw no clear indication for the detected wavelength changes, and attributed these false detections to wavelength measurement error that arises due to the irregularity of the banding.

S4.3 Wavelength-migration and wavelength-slope correlations
We computed the correlation between wavelength and migration distance at all study areas, which are shown in Table S2. Since measurements were computed in overlapping windows and are geographically proximate, a large degree of autocorrelation is expected. Standard significance tests for uncorrelated data depend on the sample size of the data, and a simple count of measurements misrepresents the true degrees of freedom in autocorrelated data. To assess the significance of the wavelength-migration correlation, we used a paired t-test for which sample size is corrected to account for spatial autocorrelation in the data [12]. The test is implemented in the library SpatialPack for R [13]. We exported wavelength, migration, and spatial position from MATLAB to CSV format, imported it into R, and used the modified.ttest function for significance calculations. We find that all correlations between wavelength and migration distance are significant at a 5% level except for SP4, which is the most heavily impacted study area.
We similarly computed slope-wavelength and slope-migration correlations, and the results are reported in Table S2). We used slope values that are closest to the center point of the window corresponding to the wavelength/migration measurement. The correlation between wavelength and slope has been empirically investigated in [11] and [7]. We found no significant correlation between slope and wavelength for our study areas except for SP3 (r = −0.34, p = 0.04). Interestingly, we found several significant correlations (at the 5% level) between migration and slope (SP2, SP5, HD2, and HD4). These correlations range from -0.34 to -0.21. To our knowledge, a correlation between these quantities has not yet been reported in the empirical literature.

S5.1 Protocol
During the course of visual inspection, we observed that bands in some areas appeared to widen over time in the direction of slope. We quantified this effect at all study areas by measuring and comparing band widths in both aerial photographs and more recent imagery. To do this, we segmented the aerial photograph to identify bands, gathered image intensity profiles along transects through the bands in direction of slope, and fit a simple step function to extract band width.
We eliminated the large-scale background variations in pixel intensity in the aerial photographs by subtracting a coarsely gaussian-blurred version of the image. We then applied a manually-tuned threshold to create a binary image of the vegetation bands. We passed this binary image to the regionprops function in MATLAB 2016b, and extracted the areas, centroids, major and minor axis lengths, and orientations of all connected components. We applied manually tuned criteria on the areas and ratio between major and minor axis length to isolate the vegetation bands in the binary image.
We then drew linear transects through the centroids of the stripes in the direction of the minor axis. We visually confirmed that the minor axis direction serves as an effective proxy for the slope direction. For all study areas, transect lengths are approximately 100 pixels. The transects were drawn so that 25% of the transect lies downslope of the centroid, and 75% lies upslope. We did this so that the same transect could be used for both the aerial photographs and the more recent imagery, accounting for band migration upslope. In order to obtain replicate measurements for estimation of error and variance, we drew eight additional transects transverse to the original (four on either side). These transects are spaced approximately 4 m apart, which precludes double-sampling of pixels by adjacent transects. We then used these transects to extract pixel intensity profiles of the aerial photography and more recent imagery. We converted color images to grayscale before extracting intensities using the rgb2gray MATLAB function. Whenever nearinfrared channels were available in the data, we computed the SAVI values for the image and extracted the SAVI intensities along transects as well.
We fit the intensity profile with simple plateau-like curves using MATLAB's nonlinear least squares curve fitting function, lsqcurvefit. The curve has the form which approaches a piecewise constant function with levels b 1 and b 2 and breakpoints at b 3 and b 4 in the limit as α → ∞. We use α = 500. To fit this function, we rescaled all transects to lie along the interval x ∈ [0, 1]. For our data, the squared error cost function typically had many local minima, and so the result was sensitive to the initial guess for the parameters b 3 and b 4 . We fit each intensity profile using 10 uniformly random initial guesses for b 3 and b 4 (such that b 3 < b 4 ), and chose among these the result with the minimum squared error. We then used the value w = b 4 − b 3 as the measured width of the band along the particular transect. If the transect had an accompanying SAVI profile, then the median SAVI value s along the interval b 3 ≤ x ≤ b 4 was recorded.
Since each band was measured using multiple parallel transects, we obtained multiple measurements of w for each band. We used a threshold on the standard deviation (σ < 0.1) of these measurements to exclude data points where bands may have substantially degraded or disappeared, or where the measurement is likely poor for some other reason. After applying the threshold, we calculated the mean w for each remaining band at each time point. When SAVI data is available, we calculate the mean s for each band as well. We verified that the measurements of width ratio and width ratio-SAVI correlations presented below are not sensitive to the standard deviation threshold chosen.

S5.2 Sool Plateau measurements
We observed appreciable increases in band width in SP1-SP4. We measured widths at multiple time points in these areas to assess when the widening may have occurred and whether it is a seasonal phenomenon. In Figure S8, we plot the distribution of band widths at SP1-SP4 over time. At SP1, SP2, and SP4, for which we have reconnaissance imagery taken in 1967, we observe that band widths changed little between 1952 and 1967. Widths are then larger in the recent imagery (onward from 2004), and do not return to their 1952/1967 widths. We conclude that band widths increased sometime between 1967 and 2004, and that this widening is not a seasonal effect.

S5.3 Width ratio-SAVI correlations
In SP2-SP4, we observed significant (at the 10% level) negative correlations between width ratio and the mean SAVI value measured along the band (Table S3). Where multiple time points were available for recent imagery, we used the images most proximate to a period of significant rainfall (Figure ??), e.g., 06/10/2004 for SP1-SP3.

S6 Model simulation
To explore how bands widen in response to parameter variation in a PDE vegetation model, we simulated a modified 1D version of the model by Klausmeier [14] with a water diffusion term added by Ursino [15]: Descriptions, units, and values of the parameters used are given in Table S4. The parameter set used for K99 is based on the values given in [14]. Parameters which are stated in [14] to differ between grasses and trees (M , J, and R) are set at intermediate values so that the spatial scale of banding resembles the scales in our regions of study. Water flow rate V was also approximately tuned to the time scale of migration in our regions of study. We set the mean annual rainfall parameter to 160 mm, which is within the range of typical rainfall levels in our regions of study. We simulated K99 using the ETDRK4 explicit pseudospectral scheme [16], with 2048 grid points, a 1000 m domain, and a time step of 0.1 years.
We perform a sensitivity analysis to estimate the linear response of band width and peak band biomass to changes in the parameter set. We began each simulation using the parameter set shown in Table S4 with an initial state of the uniform equilibrium value plus small-magnitude spatial noise. We evolved the initial state to 10,000 years to obtain an equilibrium migrating patterned state. We then began a set of perturbation simulations, where in each we perturb one value in the parameter set listed in Table S4 by a percentage between 5 and 100% that is manually tuned to produce a 5-10% response in width ratio. We then evolve the initial equilibrium patterned state by 50 years. The resulting patterned states are pulselike, and we measured widths by thresholding using a small value (10 −2 ) and counting the lengths of the resulting connected components. We stored the maximum values for each band as well. We show the width ratios and height ratios for all parameter perturbation simulations in Figure S9, where the ratios are computed by dividing band widths or band peak values in perturbed simulations by the widths or peak values from the initial patterned state. Increasing A, J, R, D n , and V and decreasing L, M , and D w results in band width increases. All parameter changes which increase band width also increase peak biomass, except for D n and To simulate a scenario where band species composition shifts from woody to grass biomass, we simultaneously increase J, R, and M by 10% and D n by 50% ( Figure S10). Although increasing mortality by itself reduces the band width, the simultaneous increase of these four parameters results in band width increase. Peak biomass decreases as well.          Table S2: Band properties measured using a modification of the Fourier window method by Penny et al. [11]. Ranges shown are the 25th and 75th percentiles. Wavelengths W 1 were measured in the R.A.F. aerial photography datasets, and W 2 were measured in recent satellite imagery datasets.
Slopes were computed from the SRTM 1 arc-second elevation dataset. Significance of correlations was assessed using a t-test corrected for spatial autocorrelation [12], and p values, t values and degrees of freedom are given in parentheses.   Table S4: Parameters for the model by Klausmeier [14]. 22