Models of dune field morphology
Back to main page

Dune Field Simulator

Click here for the software of the model described on this page.

Brad Werner (1995) presents a cellular automaton that is able to produce many of the typical morphological features of aeolian dune fields and yet is surprisingly simple. Below are satellite (Google Earth) images of real dune fields (left) compared with model dune fields (right).

Fig 1. Dune fields, real (Google Earth, left) and simulated (right)

Werner's model was extended by many researchers. I've created software to implement Werner's model and several extensions, described below.

Werner's Model

In Werner's model, a dune field is represented by a grid of square cells, each containing a stack of slabs representing the height of sand at that cell. A cell is in a 'wind shadow' if a stack on cells upwind of it would cast a shadow upon more than half of the top of the stack on that cell, when the angle of lighting is 15 degrees (approximating the extent of the recirculation cell downwind of a slipface). Thus in figure 2 below there are four cells that are in wind shadows (cells with top surfaces more than 50% in shadow), with wind from the left.

Fig 2. Shadows cast at 15 degrees by sand slabs. Four cells are more than half shadowed.

Sand transport is modeled by repeatedly selecting a cell at random and, if the cell has sand and is not in shadow, taking a slab from that cell and moving it downwind a fixed number of cells (the saltation hop distance). Saltating grains are more likely to bounce from hard surfaces than from sand or if entering a wind shadow, so the sand is deposited at that new cell with a probability 'p' that depends on whether the new cell is bare (p = 0.4 of deposition), already sand-covered (p = 0.6), or in wind shadow (p = 1, ie., always deposited). If the deposit does not occur (a bounce), that sand is moved farther downwind by the hop distance and again the 'dice are tossed' to see if the slab will deposit. This repeats until eventually it is deposited, and then another grid cell is randomly selected.

Slabs that reach a boundary of the field are wrapped to the opposite end (closed periodic grid). If deposition or erosion creates a slope greater than the angle of repose, the appropriate surrounding slab(s) are moved downslope one cell (avalanched), expanding outward as necessary. To make it easy to check for the need to avalanche, sand slabs are considered to have an aspect ratio of 1:3 (height:width), so that a difference of 2 units in height is a reasonable angle of repose (tan-1 2/3 = 33.7 degrees). An avalanche can occur in any direction, and this mechanism is the only opportunity for sand to move laterally. The lateral edges of the grid are joined by wrapping to the opposite side (periodic boundaries), which is made easier for the software by restricting the grid size to power of 2.

Below is the main loop for my interpretation of Werner (1995).

    public override void Tick() {
      for (int subticks = LengthDownwind * WidthAcross; subticks > 0; subticks--) {
        int x = rnd.Next(0, LengthDownwind);  // get coordinates [w, x] for a random cell
        int w = rnd.Next(0, WidthAcross);
        if (Elev[w, x] == 0) continue;    // if the cell is bare, get another cell
        // if (Shadow[w, x] > 0) continue;   // if the cell is in shadow, get another cell; not in Werner(1995)
        erodeGrain(w, x);                 // remove slab from this cell (also do any needed avalanching)
        while (true) {                    // repeat until slab is deposited or lost
          x += HopLength;                 // jump downwind by saltation hop length
          if (x >= LengthDownwind) {      // if past end of grid...
            if (openEnded)                // exit loop if open ended (discard slab)
              break;
            x &= mLength;                 // otherwise, wrap to start of grid (grid length is a power of 2)
          }
          if ((Shadow[w, x] > 0) ||       // if new location is in shadow, or 
              (rnd.NextDouble() < (Elev[w, x] > 0 ? pSand : pNoSand))) {  // if slab sticks (not a bounce)
            depositGrain(w, x);           // then deposit slab (also do any needed avalanching)
            break;
          }
        }
      }
    }

(Click here for more source code and detail.)

Fig 3. Moore (left) and Von Neumann (right) neighbourhoods

When checking a cell for avalanche conditions, the model could examine the slope angle of four immediately-adjacent cells (a Von Neumann neighbourhood) or all eight surrounding cells (a Moore neighbourhood); the alternatives are shown in figure 3. With a 3:1 aspect ratio, the slope angle is greater than 33.7 degrees when the difference in height is greater than 2 slabs (for adjacent cells, tan-1 3/3 = 45 degrees, and for corner cells, tan-1 3 / (3 * √2) = 35.3 degrees).

Also, when checking neighbouring cells for slope, the checking can be in a fixed sequence each time ('deterministic') or in a random sequence ('stochastic'). For a Moore neighbourhood, whether checked in a fixed or random sequence, it makes sense to always check the four steeper neighbours first.

The choice of neighbourhood and sequence affects morphology. Figure 4 below shows 256 x 256 portions of 512 x 256 dune fields at T=500 with identical starting conditions (average sand height of 3 randomly distributed) and with the four neighbourhood variants. A Von Neumann neighbourhood tends to create narrower, more angular dunes.

Werner's model appears to use a deterministic Moore neighbourhood (eight cells) when checking for angle of repose violations, judging from the size and shape of barchan dunes in his figure 2, copies of which appear at the right side of the figure below (differing scales, presumably).

Fig X. The effect of neighbourhood type (output of Werner's model at right; from Werner (1995))

Sand that bounces into a wind shadow is deposited with certainty (p = 1), so it follows that sand in a wind shadow would be unerodable. However, Werner doesn't mention that being part of his model and he doesn't publish code fragments.

Dune fields in Werner's model typically evolve from uniform or random initial conditions -- roughly the sane hieght of sand everywhere -- into bedforms that steadily decrease in spatial frequency by merging and increasing in size until there is just one large transverse ridge moving downfield.

If the initial sand supply is low (eg., less than 3 units of sand per cell, average), barchan-like dunes form initially and then merge into transverse ridges. If the sand supply is higher, ridges form immediately, evolving from high defect levels through fewer and fewer ridges. Most morphologies 'mature' to transverse ridges within about 1000 cycles.

The two video clips below show examples of the main two evolutionary patterns. The first scenario has a low initial sand supply (2 units on average, randomly distributed), and the second has a higher supply (10 units on average, randomly distributed). These run my interpretation of the Werner (1995) model, plus no sand erosion in wind shadows and a determistic Moore neighbourhood (eight cells).

The dune field is shown plan view (ie., from above) in the main window, with wind blowing from the left. Above the dune field is a cross-sectional profile (unequal horizontal and vertical scales) of the top of the dune field at the triangular marker. The scales of both the dune field and the cross-section were set to auto-adjust to current conditions, so they change as the field evolves. The shaded areas in the cross-section window indicate wind shadows.

That a few simple rules can give rise to such interesting morphology is impressive.

Below is the same configuration but with greater initial sand height (10 units instead of 3). Note the absence of a barchan dune phase. With even more initial sand, the ridges are more linear, sooner.

These samples are 512 cells x 512. Dune fields with smaller dimensions evolve similar patterns in a fraction of the time.

The videos help emphasize that Werner's model evolves toward a single transverse ridge. The interesting stages, as usually captured in figures of journal articles, are transient. If they could have enough time and space, dune fields might evolve as described by the model, but wind ripples do not seem to continually merge, growing in size and reduce in frequency -- instead, they seem to quickly reach a steady state. The model extensions discussed below address that.

The symmetrical cross-section of 'mature' bedforms in Werner's model is largely due to avalanching in both the upwind and downwind direction. If the model is run with avalanching prevented in the upwind direction, the bedforms are asymmetrical. Figure 5 shows cross-sectional profiles (streamwise) of models with normal avalancing (in the direction of the steepest-slope) and with avalanching prevented in the upwind direction.

Fig 5a. Werner's model with normal avalanching
(vertical and horizontal scales unequal)
Fig 5b. Werner's model with avalanching upwind prevented
(vertical and horizontal scales unequal)

Jon Pelletier (2007) publishes an interpretion of Werner (1995) that includes protection from erosion in shadow zones and Von Neumann neighbours. Pelletier provides code in his book (be sure to check his web page for code corrections) and also online at the Community Surface Dynamics Modeling System repository.

Momiji's Extensions

Fig 5a. Werner (1995) generates roughly-symmetrical dunes
(vertical and horizontal scales unequal)
Fig 5b. Momiji (2000) generates more realistic dune profiles
(vertical and horizontal scales unequal)

Hiroshi Momiji et al. (2000) built upon Werner's model, calling for no erosion in wind shadows and introducing 'wind speed-up' on the stoss surface proportional to surface height, which helps produce dunes with more realistic asymmetric cross-sectional profiles that have much less tendency to grow forever in size. Figure 5 to the right shows two simulations, one of Werner model (top) and the other of Momiji's model (bottom). The dune field is 512 cells long by 64 cells wide, wind blowing to the right. The top half of each sub-figure shows the cross-sectional profile, streamwise, through the dune field (plan view), at the triangular marker. The shaded zones in the cross-section view are shadow areas. There is vertical exaggeration on the cross-sectional view.

For cells with sand heights signficantly below the dune field average, Momiji's wind speed-up algorithm produces windward (reverse) sand transport, perhaps simulating transport by recirculation cells. This limits how deeply the dunes can erode.

Steven Bishop et al. (2002) extends Momiji (2000), with more detailed documentation of their model. Bishop implies their model uses a Moore neighbourhood (eight cells).

Wind speed-up, including wind reversal in troughs, tends to limit the size of dunes and thus greatly delays evolution toward a single transverse ridge. After over a million cycles with deep initial sand on a 1024 x 1 dune field, there are still a half dozen bedforms, whereas dunes in Werner's model coalesce in about 150,000 cycles.

Baas's Extensions (DECAL)

Andreas Baas (2002) extends Werner's model to accommodate the effect of vegetation. Jumping L cells per saltation hop might bypass a cell with vegetation, so in his model sand jumps only one cell per saltation hop. Joanna Nield & Baas (2008) (the DECAL model) says this results in streamwise compression of the model space but "does not contribute anything additional to the dune field patterns". Cells are selected randomly as in Werner's model, but are only allowed to be selected once per cycle, which they call "polling without replacement" (whereas Werner's model may select a particular cell repeatedly or never in one cycle). This is said to enhance the roughness of the resulting morphology.

Their model checks for avalanching using a Von Neumann neighbourhood (four adjacent cells), resulting in characteristic V-angled barchans.

Below, in figure 6, the streamwise compressing effect is explored. Figure 6a shows Werner's model (plus no erosion in shadows) with L=5, Von Neumann neighbours, initial sand height of 100, at T=200. Figure 6b shows the same arrangement but L=1. Figure 6c shows the figure Xa compressed to 48% (using Photoshop) such that the bedform frequency roughly matches that of the bedforms in figure 6b.

Saltation hop length of 5 vs 1: Sand height 100, T=200

Fig 6a. Saltation hop L = 5
Fig 6b. Saltation hop L = 1
Fig 6c. Figure 6a (L=5) compressed to 48% using Photoshop

Reducing the saltation hop length to 1 from 5 appears to compress the dune field by about half. The differences are more pronounced with limited sand (initial average height = 2), as shown in figure 7 below.

Saltation hop length of 5 vs 1: Sand height 2, T=200

Fig 7a. Saltation hop L = 5
Fig 7b. Saltation hop L = 1
Fig 7c. Figure 7a (L=5) compressed to 48% using Photoshop

Adjusting the cross-field dimension might be appropriate too, to account for the difference in transport angle when a slab is avalanched one cell laterally after a hop length of 1 cell compared with a hop of 5 cells.

Pelletier's Extensions

Jon Pelletier (2009) enhances Werner's model along the lines of Momiji et al (2000), but instead of adjusting the saltation length based on height, Pelletier's model adjusts the probability of erosion and deposition at each grid point, at each time tick, according to the local sand flux, which is dependent on bed shear stresses of the then-current topography. This results in bedforms that evolve into a dynamic steady-state train of bedforms (unlike earlier models where the bedforms coalesce into one transverse ridge), and suggest a linkage between physical parameters and bedform size.

References

Baas ACW, 2002. Chaos, fractals and self-organization in coastal geomorphology: simulating dune landscapes in vegetated environments. Geomorphology 48(1-3): 309-328.

Bishop SR, Momiji H, Carretero-Gonzalez R, Warren A, 2002. Modelling desert dune fields based on discrete dynamics. Discrete Dynamics in Nature and Society 7: 7-17.

Momiji H, Carretero-Gonzalez R, Bishop SR, Warren AW, 2000. Simulation of the effect of wind speedup in the formation of transverse dune fields. Earth Surface Processes and Landforms 25, 905-918.

Nield JM, Baas ACW, 2008. Investigating parabolic and nebkha dune formation using a cellular automaton modelling approach. Earth Surface Processes and Landforms 33, 724-740.

Pelletier JD, 2008. Quantitative Modeling of Earth Surface Processes. Cambridge University Press.

Pelletier JD, 2009. Controls on the height and spacing of eolian ripples and transverse dunes: A numerical modeling investigation. Geomorphology 105, 322-333.

Werner BT, 1995. Eolian dunes: computer simulations and attractor interpretation. Geology 23, 1107-1110.

Software

Source code for the above model is described and available here.

Author

Jim Elder
Ottawa, Canada
2009 Oct
contact

Creative Commons License
This work is licensed under a Creative Commons Attribution 2.5 Canada License.

Evolved morphology from initially moderate (h=10) sand in lower portion of the field, low (h=1) elsewhere