A coupled vegetation/sediment transport model for dryland environments

Dryland regions are characterized by patchy vegetation, erodible surfaces, and erosive aeolian processes. Understanding how these constituent factors interact and shape landscape evolution is critical for managing potential environmental and anthropogenic impacts in drylands. However, modeling wind erosion on partially vegetated surfaces is a complex problem that has remained challenging for researchers. We present the new, coupled cellular automaton Vegetation and Sediment TrAnsport (ViSTA) model, which is designed to address fundamental questions about the development of arid and semiarid landscapes in a spatially explicit way. The technical aspects of the ViSTA model are described, including a new method for directly imposing oblique wind and transport directions onto a cell‐based domain. Verification tests for the model are reported, including stable state solutions, the impact of drought and fire stress, wake flow dynamics, temporal scaling issues, and the impact of feedbacks between sediment movement and vegetation growth on landscape morphology. The model is then used to simulate an equilibrium nebkha dune field, and the resultant bed forms are shown to have very similar size and spacing characteristics to nebkhas observed in the Skeleton Coast, Namibia. The ViSTA model is a versatile geomorphological tool that could be used to predict threshold‐related transitions in a range of dryland ecogeomorphic systems.


Introduction
Aeolian landscapes in dryland regions are dominated by the interactions between ecological processes and sediment transport processes [Baas and Nield, 2007]. Wind regime and sediment supply are known to largely control the formation of bare sand dunes [e.g., Wasson and Hyde, 1983;Werner, 1995;Momiji et al., 2000;Bishop et al., 2002], but the dynamic interactions between vegetation and sediment movement, which form part of a wider, multiscaled interplay of biophysical and anthropogenic factors (Figure 1), are less well understood [Schlesinger et al., 1990;Hesp, 2002;Breshears et al., 2003;D'Odorico et al., 2013;Nield and Baas, 2008a;Ravi et al., 2010;Yan and Baas, 2015;Li et al., 2017]. Predictive modeling of ecogeomorphic dynamics is crucial for improving our understanding and management of potential dryland landscape responses to both environmental and anthropogenic influences. Applications include mitigation of land degradation [Barbier et al., 2006;D'Odorico et al., 2013], climate change impact assessments [Thomas et al., 2005;Seager et al., 2007;Wang et al., 2009;Bogle et al., 2015], and investigating the effects of changes in land use [Levin and Ben-Dor, 2004;Middleton and Sternberg, 2013] in drylands.
(e.g., intensive agriculture and grazing, severe drought, and fire) can induce discontinuous transitions leading to hysteretic vegetation growth [Scheffer et al., 2001;Suding and Hobbs, 2009]. Bailey [2011] demonstrated that quasi-stable (vegetated and bare soil) population states can appear within a bistable region of a highly stressed system. Depending on whether the starting point is a full cover of vegetation or bare ground, Bailey [2011] showed that different vegetation pattern sequences can evolve in response to the same environmental stress.
Wind erosion models have played a key role in simplifying the complexities involved and have helped develop our understanding of wind flow and sediment transport dynamics on partly vegetated surfaces. Nevertheless, previous attempts to evaluate the impact of vegetation on erosivity and erodibility are challenging. Dune mobility potential over long timescales has been successfully modeled using simple wind-based indices [e.g., Thomas et al., 2005;Wang et al., 2009], but such an approach does not resolve sediment transport dynamics at the same scale as plant-flow interactions. Traditional shear stress partitioning methods [e.g., Schlichting, 1936;Raupach, 1992;Raupach et al., 1993;Marticorena and Bergametti, 1995] rely heavily on the dimensionless roughness density (λ) and therefore suffer from the inherent requirement for a given surface to be treated homogenously. In reality, the spatial distribution of vegetation strongly affects horizontal sediment flux, even at low roughness densities [Okin and Gillette, 2001;Gillette et al., 2006;Lancaster et al., 2010;Li et al., 2013]. Shear stress partitioning models may also only apply to the specific type of vegetation being investigated.
An alternative model by Okin [2008] emphasizes gaps in vegetation structure, rather than assuming a homogenous mean value of vegetation cover, and thus allows aeolian erosion thresholds to vary across a surface. This notion is supported by field observations that flux is not ubiquitous across an entire landscape during a transport event [e.g., Gillette et al., 2006]. The Okin model shows good agreement with field measurements at the plant/patch scale [Li et al., 2013;Sankey et al., 2013]. Li et al. [2013] derived a set of optimum model Journal of Geophysical Research: Earth Surface 10.1002/2016JF004096 parameters using a large multisite data set and found that the Okin model performed better than other parameterization schemes used in previous studies. However, while Okin's [2008] approach allows individual unvegetated gaps as well as entire landscapes to be simulated, it aggregates the total sediment flux produced over a given area. Moreover, the Okin model must be supplied with information on gap sizes and distributions, rather than allowing these to develop in a separate vegetation distribution model.
Computational fluid dynamics (CFD) modeling has been used to explicitly resolve turbulent wind flow and sediment flux around individual roughness elements to a high temporal and spatial resolution [e.g., Turpin et al., 2010;Dupont et al., 2013Dupont et al., , 2014. This is promising given the increasing evidence for the importance of turbulence in wind flow on semivegetated surfaces [e.g., Leenders et al., 2007;Mayaud et al., 2016c]. However, the computational expense of CFD limits its ability to simulate long-term erosion at anything larger than the field scale. Minimal dune models [e.g., Schwämmle and Herrmann, 2005;Hermann et al., 2008], which combine analytical descriptions of turbulent wind velocity with continuum saltation models, have been used to successfully simulate transitions between barchan and parabolic dunes [Durán and Hermann, 2006;Hermann et al., 2008]. However, such continuum models are also computationally expensive to run. A new modeling approach is therefore needed that recognizes the heterogeneous nature of semivegetated dryland surfaces, while allowing for landscapes to be simulated over large areas and long timescales.
The use of cellular automaton (CA) models, a class of numerical models consisting of a regular, discrete grid of cells operated on by predefined rules, has been particularly productive in the context of dune field patterning [e.g., Werner, 1995;Momiji et al., 2000;Bishop et al., 2002;Narteau et al., 2006]. The spatial nature of CAs allows distributed dynamics to be represented in a versatile way, which lends itself well to analyzing ecogeomorphic systems [Baas, 2007]. Importantly, CAs can be employed to capture full landscape-scale dynamics using only fundamental processes, without the need for intricate grain-scale processes of sediment transport. The Discrete ECogeomorphic Aeolian Landscape model (DECAL) [Baas and Nield, 2007;Baas, 2008a, 2008b], later extended by Yan and Baas [2017], provides an approach for simulating the dynamics of vegetation-dependent dunes such as parabolic dunes and nebkhas (coppice dunes), by incorporating feedbacks between sedimentation balance and vegetation growth. DECAL introduces a length scale to resultant landforms that remain otherwise dimensionless in bare sand CAs. However, the biotic component in DECAL in its original form only incorporated as a generic, unitless "vegetation effectiveness," which acts as a proxy for the aerodynamic effect of plants on sediment transport. Plant-plant interactions are therefore not explicitly modeled, and the effects of vegetation form and porosity on local-scale wind flow are not considered. Yan and Baas [2017] extended DECAL to explicitly relate the vegetation effectiveness to the size and growth of individual plants. This allows seasonal sediment transport and vegetation growth to be simulated, but seasonality is imposed directly in the model, rather than emerging within the vegetation system in response to environmental or anthropogenic stress.
In this paper we present a new Vegetation and Sediment TrAnsport model (ViSTA). The ViSTA model couples a CA vegetation distribution model with a CA sediment transport model, allowing explicit feedbacks among vegetation growth, wind flow dynamics, and sediment flux over a surface. The evolution of landscapes with different distributions, scales, and geometries of vegetated patches can therefore be simulated for the first time. Here we verify and validate ViSTA's performance using a variety of experiments. The aim of developing ViSTA is to predict transitions from vegetated, geomorphically inactive landscapes to more active unvegetated landscapes in globally significant dryland ecogeomorphic systems.

Model Setup
The ViSTA model consists of two coupled models that interact with each other over various timescales ( Figure 2a): (i) a vegetation model that simulates vegetation growth and declines in response to environmental stresses resulting from climate and land use changes, and (ii) a sediment transport model (composed of two modules) that moves sediment across the model domain in response to spatially varying wind speeds. The coupled scheme is formulated such that the distribution of vegetation (Module 1; Figure 2b) alters local wind flow characteristics (Module 2a; Figure 2c), thus impacting sediment flux patterns over the surface (Module 2b; Figure 2d), which in turn affects vegetation growth through ecological feedbacks (see Text S1 and Figure S3). The ViSTA model also incorporates several submodules that can be activated to simulate Journal of Geophysical Research: Earth Surface 10.1002/2016JF004096 herbivore/grazing, fire, and precipitation events. These events have a primary impact on the state of the vegetation in Module 1, although precipitation can also affect sediment transport directly by changing soil moisture levels. Since the focus of the work presented here is the aeolian component of dryland geomorphology, we do not explicitly model neither the movement of water across the landscape nor the associated fluvial sediment transport. However, the benefits of a modular approach mean that fluvial processes can be incorporated into future versions of the ViSTA model.
All modules in the ViSTA model rely on local neighborhood operations to produce dynamic responses from basic rules centered on each discrete part of the grid. All grid cells hold a variety of attributes (including sand depth, soil moisture and nutrient levels, and vegetation characteristics such as plant type, height, and porosity) that are altered by applying transition rules during each time step. The state of all cells at the end of that time step becomes the initial state of all cells at the beginning of the next time step. Cells interact with adjacent cells based on their local von Neumann neighborhood. Since sediment moves regularly across the ViSTA model domain, a wrap-around torus grid (i.e., periodic boundaries) is used to preserve sediment within a closed system. As such, there are no external sediment sources apart from the initial supply to the domain, and sediment is internally "recycled" between adjacent landscape surfaces (although there is the possibility within ViSTA to specify nonperiodic boundaries, in order to explore sediment supply effects [e.g., Eastwood et al., 2011]). While dryland landscapes will gain and lose sediment externally to the system, this is significant only over long (> 100 years) timescales [Bailey and Thomas, 2014]. Running the model on a large enough domain (≥ 10 2 m) ensures that edge effects are not propagated through the system in any meaningful way.
The model is implemented in the Python® programming language (Python 3.5.0 64 bits, Qt 4.8.7, PyQt4 (API v2) 4.11.4 on Darwin; Spyder® development environment), using code written by the authors. The model code is available upon request from the authors.

Vegetation Growth (Module 1)
The vegetation module consists of an adapted version of Bailey's [2011] CA model for simulating semiarid vegetation. In brief, a neighborhood is defined in the model that is used to calculate local plant-interaction effects, employing a similar approach to Thiery et al. [1995]. Grid cells accommodate plant biomass that is representative of three vegetation types: grasses, shrubs, and trees. Cells are either occupied (alive) or unoccupied (dead and therefore bare), and plant type and age are also recorded. A cell only accommodates a single plant type at any given time, so the spatial resolution of the module should ideally be of the same order of magnitude as real plants (i.e., 10 À1 -10 1 m). Seven factors govern a given cell's probability of survival: (i) neighborhood effects, (ii) response of vegetated cells to precipitation, (iii) cell biomass, (iv) cell age, (v) sediment balance (i.e., plant response to sediment erosion/deposition), (vi) grazing, and (vii) fire. Bailey [2011] originally formulated the dependence of generic plants on generic environmental stress and neighborhood effects as a function of cell age. However, since vegetation growth effectively depends on photosynthetic gains being higher than losses due to respiration, age-dependent growth does not account for situations where plants have not experienced fully optimal growing conditions throughout their lives and thus do not have a biomass that is in proportion with their age [Danin, 1996]. To correct for this ecological process, the majority of dependencies in Module 1 of the ViSTA model are formulated as a function of arbitrary growth units, such that individuals move along a predetermined nonlinear "growth pathway" to provide equivalent changes in biomass [e.g., Yan and Baas, 2017]. If a plant experiences optimal conditions (i.e., sufficient precipitation-derived moisture) at every module iteration, it will grow by the same number of growth units as the period of time (in months) represented by that iteration. If a plant experiences suboptimal moisture conditions, it will grow by a fraction of a full growth unit that reflects the available moisture. If precipitation conditions are particularly poor, plants can lose growth units and move back down their growth pathways. In other words, precipitation history determines the rate at which a plant travels along its growth pathway. This scheme means that plants can live for long periods but accrue relatively little biomass if conditions are harsh.
In the ViSTA model, the biomass of a plant (B) determines the strength of competition/facilitation that it exerts on its neighbors, as well as its sensitivity to competition/facilitation from surrounding plants. The neighborhood of a given cell is defined by five concentric shells in an extended Moore neighborhood, and each shell is assigned a coefficient (c i ), with positive values representing net facilitative effects at close distance (shells 1 and 2), negative values representing net competitive effects (shells 4 and 5), and zero being neutral (shell 3). The biomass-dependent contributions from all the occupied cells in a given cell's neighborhood are summed and combined with stress from incoming precipitation, (ρ, mm). This is combined with other stresses (due to accumulated drought (D S ), due to age (Z S ), due to sedimentation balance (Sed S ), due to grazing (G S ), and due to fire (F S )), which themselves depend on plant type. The total compound stress determines a given plant's likelihood of dying and thus being replaced by a new plant.
For Module 1, the temporal scale of plant-plant interactions can be adapted by changing the number of months represented by a single iteration of the module. Given the "pulse-activity response" regime observed in many dryland systems [Wand, 1999], it is important to allow plants to respond sufficiently quickly to individual rainfall and sediment movement events in the model. However, the temporal resolution of the model must also be coarse enough for realistic lags in vegetation growth [e.g., Milton and Dean, 2000; Journal of Geophysical Research: Earth Surface 10.1002/2016JF004096 Knight et al., 2004] to be preserved. Therefore, following a variety of preliminary tests, we set the standard vegetation update frequency (i.e., the frequency at which Module 1 is run) at 3 months in the ViSTA model, which fully captures seasonal variability. This matches the updating interval used by Yan and Baas [2017] in their extended DECAL model.
Text S1 provides a full description of Module 1.

Wind Dynamics (Module 2a)
The relative levels of erodibility and erosivity broadly determine the movement of sediment across the model domain. Sediment flux in the ViSTA model is formulated as a function of horizontal wind velocity (u) as opposed to the more commonly employed wind shear velocity (u * ), because several aeolian studies [e.g., Sterk et al., 1998;Schönfeldt and von Löwis, 2003;Leenders et al., 2005;Wiggs and Weaver, 2012] have shown that sediment transport is better correlated with wind velocity fluctuations than time-averaged wind shear velocity.
An unobstructed wind velocity (u ref ) is derived at each time step, which represents the mean flow over a flat, sandy surface with no elements. The stochastic nature of airflow [e.g., Rice et al., 1999;Böhner et al., 2003;Klose and Shao, 2012] is introduced in the ViSTA model by approximating a Gaussian frequency distribution (μ = u ref , σ = 0.1), from which a wind velocity is randomly chosen for each cell in the domain. This wind velocity is further adjusted depending on vegetation morphology and surface topography, as described below.
Zones of reduced wind velocity in the ViSTA model are formulated as rectangular "corridors" of recovery, following the theoretical approach of Okin [2008] (see Figure S4). All vegetation elements distributed on the surface are assumed to produce a downwind zone of affected wind velocity. Previous studies [e.g., Perera, 1981;Vigiak et al., 2003;Leenders et al., 2011;Gillies et al., 2014;Wu et al., 2015;Mayaud et al., 2016c] have shown that wake zone dynamics are primarily determined by the height, width, porosity, and species type of elements. Element height in the model is directly integrated into the formulation for recovery, and the impact of element width is implicitly accounted for in Okin's [2008] corridor approach. The ratio between the flow passing through the element ("bleed flow") and the airflow that diverges over the element is determined by plant porosity, which thus controls the extent to which wind speed is reduced through the element canopy and the rate of recovery in its lee [Mayaud et al., 2016c].
The type of vegetation (here grass, shrub, and tree) determines whether an element acts to increase or decrease wind velocity in its wake. Grasses and shrubs are collectively defined as elements with growing points that are mainly located at ground level, while trees are distinguished as having a distinctive trunk below a canopy [Leenders et al., 2007]. In the cases of grasses and shrubs, it has been shown in field and wind tunnel experiments [e.g., Bradley and Mulhearn, 1983;Vigiak et al., 2003;Leenders et al., 2007;Gillies et al., 2014;Wu et al., 2015, Mayaud et al., 2016c, 2016d] that the surface wind velocity in the wake (u surf ) is lower than in the absence of plants (u ref ) and exponentially recovers with increasing downwind distance. Therefore, in the ViSTA model a zone of reduced wind velocity in the wake of grasses and shrubs is parameterized using an exponential curve: where u 0 is the minimum wind velocity in the direct lee of the nearest upwind element, x h is the downwind distance from the nearest element in terms of element height, and b is a fitted coefficient. Physically realistic values for u 0 and b are derived using the field data of Mayaud et al. [2016b], such that the recovery curve is solely determined by the height of the element (h) and its optical porosity (θ i ) (see supporting information Text S2).
In the case of the wake zone in the lee of trees, u surf is higher than u ref in the immediate lee due to a "bottom gap" effect, which arises because of the compression of streamlines between the underside of the tree canopy and the ground [Gross, 1987;Kim and Lee, 2002;Dupont et al., 2014;Mayaud et al., 2016c]. This zone of increased wind velocity is parameterized using a logistic curve: Journal of Geophysical Research: Earth Surface where a and c are fitted coefficients. The tree trunk height dictates the extent to which airflow is compressed (i.e., coefficient a), while the height of the tree canopy determines the recovery distance back to the unobstructed wind velocity (i.e., coefficient c).
Topographic steering and flow field compression are known to play an important role in the morphology of aeolian bed forms such as dunes at a variety of scales [Lancaster et al., 1996;Wiggs et al., 1996;Walker and Nickling, 2002;Hesp et al., 2005;Baddock et al., 2011;Weaver and Wiggs, 2011;Bauer et al., 2012;Barrineau and Ellis, 2013;Chapman et al., 2012Chapman et al., , 2013. In the field, wind flow has been shown to accelerate in a broadly linear manner over the stoss slope of dunes [e.g., Mulligan, 1988;Weaver and Wiggs, 2011]. Consequently, a wind speedup/slowdown effect is introduced in the ViSTA model by using an additive factor that scales linearly with the length of a given wind path along a slope (see Text S2).
Many field and modeling studies [e.g., Wasson and Hyde, 1983;Werner, 1995;Bishop et al., 2002;Lynch et al., 2009;Rubin and Hesp, 2009;Jackson et al., 2011;Bauer et al., 2012;Chapman et al., 2012;Gao et al., 2015b] have emphasized the importance of wind approach angles for the development of dunes, due to the alteration of near-surface wind velocity and turbulence, and thus sediment transport magnitudes and pathways. Simulating different transport directions is therefore vital in any realistic model of erosion at both the plant and dune field scales [Nield and Baas, 2008a;Dupont et al., 2014]. Broadly, two methods exist for simulating different directions in cellular automata: (i) directly, which involves changing the direction slabs move across the grid, or (ii) passively, which involves the rotation of the resultant landforms in a fixed transport direction. Previous implementations of these methods have suffered from several limitations. Direct methods are restricted by the small number of discretized oblique directions possible or produce dunes that are elongated perpendicular to the wind [Nield and Baas, 2008a]. Passive methods require slabs to be artificially inserted or excluded after rotation to maintain volume continuity [Narteau et al., 2006], although edge effects can be minimized if the model domain is large enough or model time is relatively short. More details on these methods are provided in Text S2.
Here we propose a new method for directly imposing continuous wind and transport directions onto any cellbased domain, based on an algorithm developed by R. M. Bailey. Accurately determining the path of the wind and transported sand across the simulated surface requires the individual path to each cell in the domain to be calculated ( Figure 3a). In effect, this parameterizes the wind as a series of linear transport "corridors" rather than a unified "front" of airflow. First, the intercept (i) of the wind path with the edge of the domain from the target cell is determined trigonometrically: where r t is the row number of the target cell, c t is the column number of the target cell, and θ is the wind angle for the given time step. Then, the length of the wind/transport path that crosses each polled cell on its way to the target cell is calculated (Figure 3a). The path distance through a given polled cell can then be trigonometrically determined by identifying sets of local intercepts (x 1 * , x 2 * , y 1 * , y 2 * ) around the polled cell, such that where r p is the row number of the polled cell and c p is the column number of the polled cell. The location of these reference points relative to each other provides a spatial definition of whether the polled cell is located along the wind path to the target cell. If it is located along the wind path, the length of the individual (hypotenuse) path crossing through the polled cell can be calculated relative to the centerline of the target cell ( Figure 3b). Table 1 summarizes the calculations for determining the hypotenuse using the local intercepts for each possible wind path scenario. Determining the specific distance of the wind path crossing each cell allows wind speedup/slowdown to be accurately simulated given the vegetation properties of each cell. It also allows for transport to be probabilistically deposited along realistic "deposition pathways" downwind of each flux-emitting cell (see Text S3).

Sediment Movement (Module 2b)
Sediment transport dynamics are introduced into the ViSTA model using a transport model based on Werner's [1995] algorithm for simulating dune formation in arid environments. Discrete volumes of sediment are transported across the model domain based on probabilities of erosion (p e ) and deposition (p d ) and the  Werner [1995], whereby some cells are polled multiple times while others are not polled at all. Polling without replacement enhances the roughness of the resulting morphology, thus providing a more realistic simulation of natural bed forms [Nield and Baas, 2008a].
During all iterations of Module 2B, a function is applied to determine how much sediment is eroded from each cell. The total amount of sediment eroded from a given cell is calculated deterministically as a function of the cell's wind velocity, using the semiempirical flux relationship of Dong et al. [2003]:

Conditions Wind Path Calculation
Positive gradient where Q is the predicted mass flux for the given time interval (kg m À1 s À1 ), u is the mean horizontal wind velocity over the given time interval (m s À1 ), u t is the horizontal velocity threshold for sediment entrainment (m s À1 ), ρ air is the density of air (1.25 kg m À3 ), g is acceleration due to gravity (9.81 m s À2 ), and a D is an empirically fitted constant. Dong et al.'s [2003] model has been shown to perform well over a variety of averaging intervals using field-based wind velocity data [Leenders et al., 2011;Mayaud et al., 2016a].
The simulated flux is converted into an equivalent volumetric flux (m 2 s À1 ), which can be specified directly in the model. For the experimental runs presented here, we assume a bulk density for sand of 2000 kg m À3 [Bailey and Thomas, 2014]. In turn, the height of the sediment "slab" moving from any given cell can be derived from the volumetric flux. The p e value of a polled cell determines the proportion of the sediment slab that is eroded from it; p e varies from 0 to 1 depending on the moisture level of the sediment in the cell and its location relative to shadow zones. The final slab height is removed from the polled cell and distributed downwind along a deposition pathway specifically calculated for each source cell.
Transport and deposition in the model occur stochastically, with each cell being characterized by a different p d value. Destination cells already populated with sediment are assigned a p d value of 0.7, while bare destination cells are assigned a lower p d value of 0.4 (see Figures S22 and S23), to replicate the increased saltation distances on hard rock through more efficient rebound [Bagnold, 1941;Bishop et al., 2002;Baas, 2008a, 2008b]. If a cell is vegetated, the p d value increases proportionally to the plant porosity (see Text S3), to replicate the more effective trapping of saltating grains by denser plants [Wolfe and Nickling, 1993;Okin et al., 2006;Suter-Burri et al., 2013]. Surface moisture (e.g., after a precipitation event) also affects sediment transport dynamically by modulating the p e and p d values via a feedback function (Text S3 and Figure S25). The total volume eroded from a source cell is divided according to the p d of each downwind cell (until a set limit of downwind cells is reached), as well as the length of the wind path traversing each downwind cell. This approach allows for more realistic deposition gradation downwind of a source than previous algorithms, where the entirety of a sand slab is deposited into a single sink cell (see Text S3).
Finally, shadow zones exist in the lee of topography Baas, 2008a, 2008b], such that no erosion (p e = 0) and complete deposition (p d = 1) occur in regions within an angle of 15°to the horizontal surface. Avalanching occurs in the direction of steepest descent to enforce an angle of repose, induced by gravity, and this angle changes according to the vegetation occupation of the cell (30°for an empty cell and 40°f or a vegetated cell). If there is a tie between two or more neighboring cells for the angle of steepest descent, an avalanching direction is chosen at random from the tied cells.
A full description of the model, including a complete list of all the parameters and coefficients used in this study, is given in the supporting information.

Model Verification
Here we present results from selected experiments designed to test the behavior of various aspects of the ViSTA model. Additional verification tests are provided in the supporting information. Note that unless stated, the following experiments are based on model runs with a single plant type (shrub). Note that the wind flow and sediment transport (Modules 2a and 2b) were switched "off" for the experiments on stable solutions (section 3.1) and the impact of stresses (section 3.2).

Stable Solutions
When forcing conditions are constant (ρ = constant), the vegetation distribution module (Module 1), populated with a single vegetation type, equilibrates over time to a stable solution, where no periodic oscillations occur in the population density (p) (Figure 4). In all cases, sensitivity to precipitation (ρ) is relatively low above the highest levels shown, with spatial coverage remaining homogeneous and dense (p > 0.9). As ρ decreases to 0, lower density stable populations are observed, and spatial patterning spontaneously emerges. For shrubs (Figure 4, middle row), the patterning follows the predictable sequence of gap ➔ labyrinth ➔ isolated patch ➔ isolated spot ➔ bare, as derived in other spatial models [e.g., von Hardenberg et al., 2001;Rietkerk et al., 2004;Bailey, 2011]. In the grass case (Figure 4, top row), as ρ decreases to 0, some spatial patterning emerges at lower stable populations, although the classic patterning sequence is weaker than in the shrub case, owing to the weaker competition amongst grasses at large radial distances. Trees (Figure 4, bottom Journal of Geophysical Research: Earth Surface 10.1002/2016JF004096 row), which have similar facilitative/competitive characteristics to shrubs, display the standard patterning sequence from gap to bare surface, although the sequence emerges at higher ρ, since tree colonization is low in harsh precipitation regimes. The vegetation system in all cases displays distinct hysteresis [Bailey, 2011], whereby the semiarid system moves along different degradation and recovery curves depending on the initial population density ( Figure S10). The stable solutions for shrubs show good visual agreement with shrub occurrence observed in satellite imagery along a dryland rainfall transect in Niger [Bailey, 2011; see Figure S11]. Quantitative spatial metric analysis shows that the model produces emergent vegetation patterns with realistic patch dimensions and distributions (see Text S4 and Figure S12).

Impact of Stresses
In order to verify the combined impact of age stress (Z S ) and drought stress (D S )-which both vary between vegetation types-the module was run with different initial vegetation proportions and in different precipitation scenarios ( Figure 5). In these runs, precipitation was distributed evenly throughout the year. In a harsh precipitation regime (80 mm yr À1 ), domains initialized with a majority of either grasses or shrubs (Figures 5a and 5b) equilibrate within~50 years to an approximately equal distribution of grasses and shrubs. This occurs because grasses colonize bare cells more efficiently than shrubs and achieve peak facilitation/completion capabilities quicker, but shrubs on average tend to survive longer due to a lengthier life cycle and an ability to better withstand accumulated drought. When the model is initiated with a majority of trees (Figure 5c), the trees die off due to their low resilience to conditions of low precipitation, and the bare cells are rapidly colonized by grasses, before the shrub population begins to accumulate.
Shrubs benefit most from an intermediate precipitation level (400 mm yr À1 ), as can be observed from the encroachment of shrubs in grass-dominated domain (Figure 5d), and the maintenance of shrub dominance over time (Figure 5e). This is due to the smaller difference in recolonization abilities of shrubs and grasses at this precipitation level, which allows communities of shrubs to become established and rapidly out compete grasses for nutrients and water [Puigdefabregas, 2005;Stewart et al., 2014]. When trees initially dominate (Figure 5f), grasses out-compete shrubs for the remaining available space, and a stable tree-grass coexistence emerges (with a higher grass proportion at slightly lower precipitation levels). This last example effectively represents a "climatically determined savanna" [Bond et al., 2003], whereby grasses persist in conditions where tree cover is restricted by water limitation [Sankaran et al., 2005].

10.1002/2016JF004096
At a high precipitation level (1000 mm yr À1 ), trees eventually dominate in the equilibrium state, regardless of the initial vegetation type proportions. In these conditions, water availability is sufficient to allow trees to approach canopy closure such that grasses and shrubs are effectively excluded. Several studies [e.g., Sankaran et al., 2005;Lehmann et al., 2011;Staver et al., 2011] have shown that such disturbance-driven savannas require disturbances external to the system to maintain a diversity of vegetation types.
A range of human disturbances-namely agriculture, grazing, and fire-is a common feature of many drylands, and these can force the system to rapidly shift to an alternative stable state. Such disturbances often reduce vegetation cover, disturb soil structure and crusting, and create large patches of bare ground [Hagen, 1991;Belnap, 1995;Nash et al., 2004;Barbier et al., 2006;Miller et al., 2012]. In turn, this affects surface roughness, cover, and patch size distribution, which have important implications for the sediment transport rates and patterns [Okin and Gillette, 2001 The impact of disturbance on vegetation growth can be investigated in the ViSTA model by changing the fire stress (F S ) and grazing stress (G S ) imposed on the domain. The dynamic equilibria reached for ρ = 80 mm yr À1 (as shown in Figure 5) are used to initialize the model, and, with the same precipitation regimes, cyclical fire and grazing stresses imposed. Figure 6 shows the impact on vegetation type and cover of changing the frequency with which fires occur, in different precipitation regimes. In a harsh regime (ρ = 80 mm yr À1 ), frequent fire occurrence leads to a relative increase in grass populations and decrease in shrub populations, as grasses burn more readily but also recolonize bare cells preferentially in these conditions.
A bistable region exists at intermediate (ρ = 150-400 mm yr À1 ) precipitation levels, where the fire frequency determines whether the vegetation system becomes grass dominated or shrub dominated. A short (5 year) fire cycle in these conditions allows short-lived grasses to continually populate the surface, whereas longer (20-40 years) cycles cannot prevent the establishment of large shrub colonies that At high (ρ = 1000 mm yr À1 ) precipitation levels, water availability allows canopy closure to eventually occur and "deterministic forest" to develop [Staver et al., 2011], regardless of fire regime. However, the fire cycle does determine how rapidly the grass population declines relative to the tree population, such that more frequent fire events enable grasses to keep colonizing bare surface over a longer period. This supports the "disturbance-driven savanna" hypothesis [e.g., Sankaran et al., 2005], whereby fire and grazing buffer against the formation of closed-canopy forests, thus maintaining both trees and grasses in a mixed system [Bond et al., 2003]. Figure 7 shows the impact on vegetation type and cover of changing the frequency with which grazing occurs, at different stocking rates (with a constant precipitation regime, ρ = 150 mm yr À1 ). This precipitation regime represents conditions in the southern Kalahari Desert, a semiarid region that is heavily used by humans for agriculture and grazing [Thomas et al., 2005]. For the grazing experiments, the stocking rates varied between 0.015 and 0.06 LSU (large stock units) ha À1 , as these are typical for Kalahari rangelands [Blaum et al., 2009;Burgess, 2009], and cycling ranges from constant grazing to one grazing period ever 20 years. When grazing occurs constantly, grass populations rapidly decline within a few decades and spot patterning emerges, which may eventually lead to vegetation collapse [Rietkerk et al., 1997;HilleRisLambers et al., 2001]. As the time between grazing events increases, grass patches have more time to become established as they recolonize grazed surfaces, leading to a greater proportion of overall grass population (e.g., 30% cover for a grazing cycle of 20 years). Predictably, increasing the stocking rate results in more rapid declines in grass populations, as well as more dramatic responses in grass/shrub proportions following individual grazing events. Greater stocking rates also result in more dramatic population density oscillations, particularly in the early stages of grassland readjustment to shrubland.

Wake Flow Dynamics
In dryland landscapes composed of patchy vegetation, the wakes developing downwind of individual vegetation elements often interact with those of their neighbors [Wolfe and Nickling, 1993;Okin, 2008]. It is essential to verify that these interactions are correctly formulated in the model. The modeled effect of patch configuration and size on wind flow recovery is shown in Figure 8. In cases where vegetation elements are relatively isolated from one another (e.g., Transects A and B), wind flow recovers fully between elements, such that the approach wind velocity for each element is the same. Once wake interference flow occurs (Transects C and D), wind flow recovers almost fully, although it should be noted that at lower element porosities, the recovery would not be complete in these configurations. Transect E is an extreme example of wake interference flow, while Transect F is an example of skimming flow, because a sufficient number of contiguous plants form a patch within which wind velocity is reduced to zero. At the point of recovery, the cumulative porosity of the plants along Transect F results in a much lower minimum leeward wind velocity (u 0 ) than along the other transects. Consequently, the recovery of wind flow in the lee of the long patch is longer than in the lee of individual elements, reflecting wind tunnel [Burri et al., 2011;Youssef et al., 2012] and field-based [Mayaud et al., 2016d] observations that the lack of lateral airflow moving around patches results in extended recovery lengths.

Temporal and Spatial Scaling
Unlike scale-less sediment transport CA models, the inclusion of vegetation in the ViSTA model requires us to consider the spatial and temporal scales at which the landscape is simulated. While several geomorphic features such as bare sand dunes, river networks, and coastlines can be considered fractal like (and thus scale less) in nature [Tarboton et al., 1988], vegetation imposes a characteristic physical scale on aeolian dune morphology. This arises because of the fundamental relationship between photosynthetic potential and size range [Nield and Baas, 2008a]. Here we present tests focused on temporal and spatial scaling in Module 2 (i.e., wind and sediment transport characteristics). Verification testing of scaling in Module 1 (vegetation growth) is given in the supporting information (Text S5).
The temporal resolution of the ViSTA model can be adapted by changing the number of minutes of abovethreshold wind velocities represented by each iteration of Modules 2a and 2b. Figure 9 shows bed form morphology resulting from increasing the time of above-threshold wind per iteration (which effectively coarsens the temporal resolution of the model, since the total model runtime stays the same). Iteration durations of 60-300 min tend to produce relatively similar nebkha landscapes and sediment transport profiles. At iteration durations of 5000-10,000 min, sediment transport rates are much lower (Figure 9, bottom left panel) than at iteration durations of 60-1000 min, due to fine-scale sediment interactions not being represented as time blocks get too large. This results in essential nebkha-forming mechanisms being bypassed, as the full range of geomorphic processes are not fully represented by the longest iteration resolutions. The clear step change in sediment transport rates and bed form definition between intervals of 1000 and 5000 min suggests that erosion and deposition processes can be feasibly treated at relatively coarse temporal resolutions until~1000 min, allowing transport mechanisms to be resolved in a computationally efficient manner.
The spatial scale of the ViSTA model can be altered by modifying the representative size of a cell ( Figure 10). Dryland mesquite-type shrubs are typically on the order of 10 0 -10 1 m in diameter , and simulations using this vegetation type produce nebkha dunes of 2-10 m diameter with a grid resolution of 1 m [Nield and Baas, 2008a]. Reducing the cell size to below the resolution of single shrubs produces more closely packed nebkha forms. A breakdown in bed form development occurs at a resolution of 0.4 m, because vegetation elements effectively exceed their physical restraints by existing closer together, thus leading to the formation of connected sand patches around unrealistically small vegetation bunches. As the cell size increases to 2 m, nebkhas of a similar size to the 1 m resolution develop, but at a coarser resolution. At cell sizes of 5 m and 10 m, nebkhas exist as isolated, single sediment-covered cells and lose their recognizable forms. Across the model domains, total transported sediment during the first~50 iterations tended to be greatest at high (0.4-1.0 m) resolutions ( Figure S27, bottom). This likely results from many individual small cells being eroded more readily than equivalent single large cells, since shadow zones are more finely resolved at the highest resolutions. However, by the end of the model runs, the total transported sediment was broadly similar across all resolutions, suggesting that an equilibrium transport rate was reached across all spatial resolutions.  Figure 8a). Bottom plot in Figure 8b shows wind recovery downwind of the last plant of each transect, for comparison purposes. Stochasticity was removed from the wind velocity algorithm for these runs.
Journal of Geophysical Research: Earth Surface 10.1002/2016JF004096

Dune Simulation
Dunes in arid sand seas are primarily classified according to sediment availability and wind directional variability [Wasson and Hyde, 1983;Pye and Tsoar, 1990]. Characterizing dune processes in the reductionist terms proposed by Wasson and Hyde [1983] is a contentious approach [see Courrech du Pont, 2015], but it provides a relatively simple way of testing sand transport models. Therefore, the ability of the sediment movement module (Module 2) to replicate basic aeolian dynamics was tested by varying wind approach angles and the initial sediment volume in the domain. Figure 11 shows that when wind is unidirectional and sediment availability is low, moving sand slabs coalesce to produce elongated barchanoid structures that propagate at different velocities according to their geometry and volume [Hersen et al., 2004]. Smaller, faster-moving structures merge with larger ones and then decelerate, leading to the emergence of a crescentic barchan dune field. A dynamic equilibrium in dune size-distribution is eventually reached, owing to the periodic boundary conditions and the small size of the system, which results in dunes propagating into their own shadows [Narteau et al., 2006]. As the amount of sand increases in a unidirectional wind regime, the barchan dunes evolve into transverse ridges, which  [Rubin and Hunter, 1982;Ewing and Kocurek, 2010;Kocurek et al., 2010;Gao et al., 2015a], with the average size and wavelength of the transverse dunes increasing during the simulation.
In moderately variable wind directions (in our case, a bimodal regime composed of an obtuse divergence angle between the primary and secondary winds), linear dunes are observed. These long, relatively straight, parallel ridges simultaneously display extension in the resultant wind direction and lateral migration [Tsoar, 1983;Rubin and Hunter, 1987]. When a multidirectional wind regime is imposed onto a deep sand bed, a star dune-dune field forms. Star dunes exhibit a hierarchy of dune forms, limited at the lower end by the elementary dune length scale [Elbelrhiti et al., 2005;Zhang et al., 2012] and at the upper end by the convective boundary layer thickness [Andreotti et al., 2009]. In the ViSTA model, when a discrete conical pile of sand is exposed to a multidirectional regime [Zhang et al., 2012], a star dune always forms with a symmetry with respect to the imposed wind directions.
An emerging body of numerical and experimental work [e.g., Courrech du Pont et al., 2014;Gao et al., 2015b;Lucas et al., 2015] suggests that sediment availability not only controls dune type but also the dune orientation, due to a variation in the dune growth mechanism. Two modes for dune orientation, the "fingering" and "bed instability" modes, have been identified by Courrech du Pont et al. [2014]. In addition to the tests described above, we examined the ability of the ViSTA model to replicate the emergence of these two modes in bimodal wind regimes, following the approach of Gao et al. [2015b]. The experimental conditions are described in detail in Text S6. In conditions of high sand availability, the bed instability mechanism systematically led to the establishment of dunes aligned in one principal direction, which grew in height through pattern coarsening ( Figure S28a). When the sand source came from a single location, straight "finger dunes" developed over a range of divergence angles, elongating with a constant orientation ( Figures S28b and S28c). Text S6 provides detailed analysis of the experimental results. The diversity of patterns produced in these experiments supports the notion that sediment availability controls both dune type and orientation [Courrech du Pont et al., 2014;Gao et al., 2015b;Lucas et al., 2015].

Coupled Model Behavior
In order to verify the feedbacks between vegetation growth and sediment movement, the fully coupled ViSTA model was run with a v-shaped precipitation regime in different wind strength scenarios. In the model runs, precipitation decreased linearly from 800 mm yr À1 to 0 mm yr À1 over the first half of a 240 year simulation, before increasing linearly back to 800 mm yr À1 . The domain was exposed to either low (0 m s À1 ), moderate (7 m s À1 ), or high (12 m s À1 ) wind strengths.
Figure 12 displays the model results over the simulation period. As the precipitation decreased, vegetation cover decreased, allowing sediment movement to gradually increase. In the strongest wind regime, the high levels of erosion and deposition resulted in a more rapid and more dramatic decrease in cover, allowing even more transport to occur, and a single, large, migrating bed form to evolve. Compared to the low and moderate wind regimes, the conditions of the high-wind regime were more conducive to grass growth at near-zero precipitation. As the precipitation started to increase again, the various wind regimes resulted in very different vegetation cover and dominant vegetation type. In low-wind conditions, shrubs rapidly colonized the patchy surface, resulting in a near homogenous and densely populated shrubland. This arises because shrubs grow more optimally than grasses and trees in the model when the sediment balance is zero. In the moderate wind regime, the intermediate levels of erosion and deposition allow a mix of grasses, shrubs, and trees to survive, with individual patches of shrubs emerging in the landscape, as is often seen in southern African drylands, for example [Bhattachan et al., 2014]. In the high-wind scenario, the migrating bed form creates an environment whereby shrubs populate the interdune area, where sediment balance is near zero, while low-lying grasses tend to dominate on the dune crest. Although this is an extreme case with a single dune, such a spatially distinct shrub-grass distribution is commonly observed in many semistable dune systems such as the Kalahari [Bhattachan et al., 2014].

Summary of Model Verification
The verification experiments presented here and in the supporting information demonstrate the ability of the ViSTA model to accurately replicate a variety of ecogeomorphic features of dryland environments, including the formation of a standard sequence of vegetation patterns in response to stress, which is hysteretic in nature; dominance of specific vegetation types in different precipitation regimes, replicating field-based and remote sensing observations of savannah/grass transitions; dynamic response of the vegetation system to fire and grazing regimes; realistic wind flow recovery in the wake of various vegetation patch configurations; the formation of standard dune forms in response to changing wind variability and sediment supply; and the impact of feedbacks between sediment movement and plant growth on landscape morphology and dominant vegetation systems. In each case, model was run over a 100 m × 100 m grid with an initial vegetation cover of 90% (with initial grass, shrub, and tree proportions derived from equilibrium conditions observed in Figure 5) and an initial uniform sediment depth of 2.0 m, for a model equivalent of 240 years. Rainfall decreased linearly from 800 mm yr À1 to 0 mm yr À1 over the first 120 years of the simulation, before increasing linearly back from 0 mm yr À1 to 800 mm yr À1 over the following 120 years. In order to balance computational efficiency with simulation accuracy and physical length scaling, we define a standard set of vegetation and sediment update parameters for the ViSTA model, which can of course be adapted depending on the needs of the numerical experiment. Wind iteration resolutions are set at 300 min, the spatial model resolution is set at 1 m, and a minimum grid size of 100 × 100 cells is employed to minimize edge effects and prevent the development of large, self-affecting bed forms. Other standard parameters (e.g., avalanching and erosion/deposition probabilities) are provided in Text S5.
We now provide results from an experiment designed to validate the output of the ViSTA model. This involves driving the model with a known precipitation regime and wind data collected in the field, and comparing the resultant bed forms with empirical field measurements.

Simulating the Evolution of a Nebkha Dune Field
Nebkha dune fields are useful testing grounds for sediment transport models, as they are characterized by relatively simple vegetation and flow conditions [see Nield and Baas, 2008a]. Here data on vegetation, dune size, and wind flow from a real nebkha dune field in the Skeleton Coast National Park, Namibia, are used to validate equilibrium simulations produced by the fully coupled ViSTA model.

Field Experiment Setup
The field site was situated in the southern half of the Skeleton Coast National Park, Namibia,~5 km inland from the coast along the ephemeral Huab River valley (20°52 0 06″S, 13°29 0 24″E [see Mayaud et al., 2016b]). In this area, relatively small nebkhas (Zygophyllum stapfii shrubs) were distributed on a large, flat gravel plain with a fetch of >2 km, with fluvial deposits upwind supplying most of the sediment. The average grain size of the sandy substrate forming the dunes was 277 μm. Aerial photography of the field site was undertaken using an unmanned aerial vehicle (UAV), to characterize typical nebkha dune spacing (Figure 13a). Terrestrial laser scanning (TLS) was also performed in situ to analyze the characteristics of three specific nebkhas (N1, N2, and N3), which were chosen to reflect the range of shrub porosity and size observed at the study site. Nebkhas N1 and N2 were located~30 m from one another, and N3 was located~400 m away on the same gravel plain. The characteristics of the study nebkhas are summarized in Table 2. For all three nebkhas, a full scan was performed of the shrub and its associated sandy tail (Figure 13b), and another scan was performed after the removal of its tail with the use of a spade (Figure 13c). This allowed the optical and volumetric porosities of the shrubs, as well as the surface and thus total volume of sand in their leeside wakes (Figure 13d), to be accurately calculated.
Wind velocity at the experimental site was measured for a total period of 10 days in early September 2014 (sampling frequency: 10 s; averaging interval: 60 s), using a two-dimensional sonic anemometer (1.5 m height) located~1 km away from the study nebkhas, on the same gravel plain. Wind flow was broadly bidirectional over the study period (peaks at 160°and 320°), and wind velocities ranged from 2 to 12 m s À1 , showing a strong diurnal pattern with large peaks in the middle of the day, similar to the velocities measured at a nearby site by Mayaud et al. [2016a]. The wind flow data are provided for reference in Figure S29.

Model Setup
We ran ViSTA in an equilibrium (as opposed to transient) modeling exercise, to determine whether the model could equilibrate to resemble the present-day landscape when forced with contemporary climatic conditions. The domain was initialized with 90% vegetation cover, to verify the model's ability to cope with vegetation conditions that are unrealistically high for the Huab Valley in current environmental conditions. The vegetation module was updated every 3 months, and the coupled model was run for a total period of 25 model years, at the end of which dynamic equilibria in plant population and sediment transport could be observed. The shrubs in the model domain were randomly initialized with porosities between 35 and 70%, to reflect the range of porosities of the three shrubs that were characterized using the TLS. A fixed annual precipitation of 60 mm yr À1 was imposed throughout the run, reflecting the hyper-arid nature of the Skeleton Coast [Lancaster, 1985], which receives occasional advective fog precipitation as cool moist oceanic air moves onshore [Viles, 2005].   The threshold wind velocity at the field site was calculated to be 5.10 m s À1 , following the method of Ravi et al. [2012]. Wind velocities at the field site commonly exceeded this threshold during the experimental period, so for the model run, 20 above-threshold (7 m s À1 ) wind events, each lasting 200 h, were imposed on the domain each model year. This corresponds to Lancaster's [1985] observations that winds along the Skeleton Coast transport sand >40% of the time during the windy season.
The approach angle in natural wind flow is rarely exactly constant, with variations of a few degrees often observed even in supposedly unidirectional winds (as was the case at the Skeleton Coast site). Some multidirectionality in the wind regime is thought to be important for nebkha dune field development, owing to its impact on the distribution of wind stress around dunes as the lee side wakes around each plant constantly change [Gillies et al., 2014]. Since the ViSTA model can explicitly vary the wind direction by fractions of a degree from one time interval to the next, variations of up to ±5°from the principal wind approach angle were chosen randomly from the same distribution of wind approach angles measured at the field site (during periods of roughly unidirectional wind flow) and used to drive the Skeleton Coast simulation. Figure 14 shows the surface morphology and vegetation cover as the simulation progresses through time.

Simulation Results
After 5 years of model time, the decline of vegetation and emergence of labyrinth-type patterning due to precipitation stress allows increasing sediment transport to occur. Oscillations in the transported volume curve arise due to minor changes in wind direction, which to a certain extent rework nascent bed forms from one interval to the next. Sediment avalanching begins to occur at around four model years, as accumulations of sediment around vegetation strips exceed critical angles.
The combination of precipitation stress and sediment balance changes leads to the development of isolated vegetation patches by 10 model years, at which point the landscape begins to stabilize (as can be seen from the equilibration of the transport and avalanching curves, as well as spatial bed form and vegetation patterns). From this point, existing nebkha dunes continue to grow and starve other parts of the landscape of available sediment and provide ideal locations for further vegetation growth on the nebkhas themselves. The nebkha dunes are anchored bed forms but do shown some dynamism, with some migrating up to 5 m downwind in the final 10 years of the simulation (similar to nebkha migration rates of up to 1.5 m yr À1 reported by Necsoiu et al. [2009]).
Three-dimensional representations of the modeled landscape are shown in Figure 15. The fully vegetated domain (Figure 15a), which is an unrealistic condition in the Huab Valley, evolves to a landscape of fully developed nebkha dunes ( Figure 15) that is qualitatively very similar to that observed in the Skeleton Coast (Figures 13a and 13b). In quantitative terms, analysis of all the individual nebkhas in the model domain at the end of the simulation reveals that dune spacing ranged from 5 to~40 m (Figure 15c), which is on the order of spacing observed in Figure 13a and in the field. Final vegetation spot heights in the model ranged from 0.3 to 1.0 m, which corresponds well to shrub heights in the field (note that the maximum height a vegetated cell can grow to is predetermined in the model setup). Sandy tails from 3 m to 8 m in length were observed in the final modeled landscape, with volumes ranging from 0.20 to 2.30 m 3 , which is similar to the three nebkha volumes determined using TLS data (0.26-1.22 m 3 ).

Conclusions
1. A vegetation distribution model and a sediment transport model are coupled, allowing the ecogeomorphic evolution of dryland landscapes to be simulated in a spatially explicit way. 2. The model requires only a few data inputs (precipitation, wind velocity, and direction) to simulate landscape dynamics. 3. A novel method is introduced in this paper for directly imposing oblique wind and transport directions onto a cell-based domain. This allows modeled landscapes to be simply and accurately exposed to natural wind velocity variations. 4. The ViSTA model successfully simulates an equilibrium nebkha dune field in the Skeleton Coast, producing realistic dune size and spacing characteristics.
The ViSTA model is now ripe for further testing and validation. For instance, we envisage applying it to further investigate nascent theories on possible "modes" for dune orientation [Courrech du Pont et al., 2014], to recreate episodes of dune activation on semivegetated farmland, and to replicate mesoscale wind flow behavior based on long-term field data.
The coupled ViSTA model represents a versatile geomorphological tool which, if parameterized correctly, can be applied to simulate a variety of dryland environments. In turn, the simplicity of the model inputs should allow for extensive investigation of future climatic and anthropogenic disturbances to dryland regions. Identifying the proximity of populated dryland systems to ecogeomorphic thresholds is directly relevant to land management policies.