A model for sedimentation of airborne particles

Download my dunes simulation package installer


Large open spaces, like deserts, offer a direct opportunity for the substrate particles (sand, silt, snow etc.) to move uniformly and form corresponding predictable dunes pattern. This regularity is not observed, however, under natural conditions and those terrains yield an unruly variety of landscapes, composed of different types of dunes (Bishop et al., 2002; Wasson and Hyde, 1983; Werner, 1995). This comes from the unsteady behaviour of the air flow when it gets closer to the surface, revealing secondary flow effects, and eventually turning into turbulent regime, which is too complex to simulate (Parsons et al., 2004a). And the dune shape is very responsive to the local wind pattern (Bishop et al., 2002; Warren and Allison, 1998; Wiggs, 2001; Zhang et al., 2000).

The author in Wahiba sands, Sultanate of Oman

Recent successful missions of Mars exploration vehicles (Crisp et al., 2003), discovery of specific dunes formations on Mars (Bourke et al., 2004; Fenton et al., 2003) and Venus (Greeley et al., 1992), as well as and specific wind regimes (dust devils and the dust cycle (Basu et al., 2004; Ferri et al., 2003; Toigo et al., 2003)) challenge researchers to go on further in developing models for desert dust transfer under variable global conditions, which may incorporate other forces like planetary gravitational and magnetic field force (Fenton et al., 2003; Spahn et al., 2003).

One of the Mars exploration vehicles, Hummer H3, undergoes thorough testing at Wahiba sands.

There exist a number of theoretical models on sediment deposition and erosion, as well as on dunes formation (see e.g. (Anderson, 1988; Bishop and Momiji, 2002; Bishop et al., 2002; Bitzer and Harbaugh, 1987; Crave and Davy, 2001; Nishimori et al., 1998; Parsons et al., 2004a; Parsons et al., 2004b; Sorensen, 2004; Werner, 1995; Wiggs, 2001). The physics of local dust movement is well understood (Anderson, 1988; Raupach and Lu, 2004); so, many of those models employ robust mathematical apparatus of computational fluid dynamics ( see e. g (Parsons et al., 2004a) and references therein).

However, all those local mechanics models have common drawback — they rapidly become too complex for the fields of larger scale, as the computational time required grows exponentially with number of particles in the system.

The goal of this paper is to describe an algorithm for modelling the global, planetary scale wind vector field, which can be further overlaid and combined with existing local-physics models when necessary. As the specialised supercomputers emerge, similar to Earth Simulator (Sato, 2004), this task becomes more achievable than ever. The importance of large scale aeolian transport modelling has been well emphasized (Livingstone, 1988, 1989, 1993; Raupach and Lu, 2004; Wang et al., 2003; Wiggs, 2001). Of all the challenges of up scaled dust transfer modelling, the proposed algorithm is dealing with surface wind velocities and, to lesser extent, with nonlinearity of grains movement. The main objective of the model is to predict sediment movement in the wind vector field of any real or desirable complexity.

The description of model

The landscape formation starts on an empty rectangular grid with altitude at each cell is equal to zero whatever units. The surface obeys periodical boundary conditions, whereby a particle crossing the lattice boundary is transferred to the same position on the opposite edge, the corresponding 3D representation of this surface is effectively being a toroid. The model resembles to some extent approaches of Werner (1995) and Bishop et al. (2002), but here the particles are being withdrawn from an infinite stock and placed under a “cursor”, which moves over the lattice following the distribution of probabilities in each individual cell of the grid. The last is determined by the local wind vector. From given definitions, the main driving force of the model is a Markov process on 2D surface (Borodin and Salminen, 2002).

The use of rather abstract notion a “cursor” emerges from the well-known quasi-ergodic hypothesis (Berry et al., 2000), widely used in chemical physics. In brief, it can be formulated as following: a single particle of a system, observed for an infinite period of time, will reveal all the possible behaviours of infinite number of particles, observed in some time chunk. Hence, tracing the path of only one particle in a stochastic force field will give one a plot of densities of probability to find a particle at given location in any arbitrary moment.

It would have been possible to account for a probability of particle anchoring as described by Werner (1995) and Bishop et al. (2002), and what would comply with the experimentally observed saltation mechanism (Iversen and Rasmussen, 1999). However, as a certain wind vector field is laid over a surface patch, the terrain height on that patch manifests the very same anchoring probability, but in an indirect way. Different rules of such an anchoring would yield then different terrain heights, without affecting their shape, representing additional yet obsolete scaling factor. Currently, the model uses 100% binding rule, which means that the terrain height in given cell is equal to the number of  “cursor” hits.

Usage of an infinite supply stock of identical particles allows one to investigate only influence of the wind field on the dunes shape, effectively eliminating the impact of sand availability, which is known to be a competitive factor in determining the dune type (Wasson and Hyde, 1983; Wiggs, 2001). The last feature could be incorporated in the future model advancements.

Simulation begins with the cursor placed in the center of grid (Fig. 1). Each cell has 8 neighbours, so there are 8 possible directions of the succeeding move. Once the cursor gets in a certain cell, the altitude there is increased by 1 unit. Given equal probabilities for either direction, the process represents common Brownian motion (Borodin and Salminen, 2002), corresponding terrain for which is shown on Fig. 2.



Figure 1. Sample qualitative redistribution of particle transfer probabilities under the action of local wind velocity (white arrow). 


Figure 2.Dune field generated by sheer Brownian motion. Terrain is composed of 11520000 particles in zero wind vector field. For best perception part of the area is evenly “flooded”.

No dunes dynamics is represented in the current state of the model, so the terrain formations are static and their height just manifests of probability of the particle to pass through the certain surface spot. The same applies to the wind field — it is fixed and does not change as the dunes grow.

Redistribution of probabilities

Let us presume that there is a “wind” in each individual cell. This wind blows in some direction (0–360º) with strength on a 0–220 whatever units scale (see appendix A). If there is “dead calm” in the cell, then, as described above this cell represents trivial case with all the probabilities of cursor move equal to 1/8, which is henceforth called the basic or initial probability (P0). Now, let us represent this P0 as a physical spring with length of 1/8 whatever units obeying the Hooke’s law, with force constant k, and being able to change its length (equal to probability) under the action of external force, determined by the wind. This spring is capable of shrinking and expanding under the action of external force, which is calculated as a scalar product of wind vector and unary vector, corresponding to one of the transfer directions:


where is the force vector,  — wind vector,  — unary vector, defining direction of transfer. Then


where Pi is the probability for the particle to go in the i-th direction.

The symmetry of the cells and properties of Hooke’s law make sure that the sum over all the 8 transfer probabilities is unity, forcing the “cursor” particle to obligatory leave the cell.

The input for the program represents a graphic PNG file with each pixel composed of three colour components: red, green and blue. As those components may take values in range 0–255, they internally manifested as vectors of corresponding length and oriented as shown on Fig. 3. The wind vector is, therefore, a sum of those colour components:




Figure 3.Sample derivation of wind velocity (white arrow) from red, green and blue components of a graphic pixel.

Given a wind vector field the model in its current state does not discriminate all possible types of particles motion (like flow, saltation, creeping and wafting (Iversen and Rasmussen, 1999; Raupach and Lu, 2004)). Current rule of probabilities redistribution (Equation 2) is a kind of a generic one. It does not include any physical information about the system. But, future extensions of the model would allow to tune this rule more accurately, even for each individual grid cell to express impact of particle size fractions and external non-equilibrium conditions. Itself, the non-zero probability of the grains to move against the local wind direction represents nonlinearity of the transfer process (Raupach and Lu, 2004). This nonlinearity may stem from the existence of more complicated phenomenae like the flow of non-spherical and/or charged particles interacting with fluctuating planetary magnetic field. For example, it was shown that stochastic particle movement (diffusion coefficient, D) depends on particles size (r) as Dµ r-4, so the diffusion can quickly override the directional physical forces as the particles get smaller (Spahn et al., 2003). To the date all this manifold of non-linear movements is manifested just in one parameter, k, of the model.

The description of program

The program consists of two intertwined parts: the core, written in Ada language (compiler GNAT 3.15 from Ada Core Technologies), and user interface written in Visual Basic 6 (Microsoft). The quality of random number sequence had been enhanced by employing the Ada implementation of Mersenne Twister algorithm (Matsumoto and Nishimura, 1998).

The interface consists of graphic window, where the generated landscape is displayed; messages window, where the program outputs all text messages; and control panel with all available controls.

The controls are as following:

“Load winds from PNG” button is to load the predefined wind vector field from Portable Network Graphic file, where each pixel is described by red, green and blue components in 0–255 range. After loading the wind field program begins to render the landscape automatically with initial default values for k = 0.95 and sea_level = 100

“Sea level” control serves to better visualise the results, by changing the depth of “flooded areas”.

“Force constant” control is to change the parameter k of the model.

It is necessary to press Enter after changing either sea_level or k values to render the image anew.

“Grayscale/color representation” and “Load palette” controls also serve to change the visualisation settings. When “color representation” is selected, the user may choose from several palettes supplied.

“Watch the point running” checkbox is the option to actually see how the cursor particle is moving under the action of the winds. In this case, however, the completion of landscape may take unbearably long time.

“Events scheduled” shows how many user commands are waiting their turn to be executed.

“Quit” terminates the program or stops simulation when “Watch the point running” is selected.

Every time a new PNG file is loaded the program creates a subdirectory named “PNGfilename_MAPS_PAPER_SUPPLEMENT”, where pictures of all generated landscapes are to be written in Windows bitmap format (BMP). Name of every file consists of predicate “gray” or “color”, depending on the visualisation colour map, followed by sea level value, underscore, and force constant k, multiplied by 100.

Sample simulations

The sample nontrivial simulations of dune formations are depicted on figures 4 through 10. The corresponding source wind fields in PNG-files can be found in the supplemental.


Figure 4. Ridges forming along the wind path when uniform wind field is overlaid (White arrow is the direction of the wind, equal for each cell).

Figure 5. Simulation of linear dunes.

Figure 6. Simulation of barchan dunes.

Figure 7. The effect of force constant k on crescentic dune formation. a. k = 0.10; b. k = 0.50.

Figure 8. The effect of force constant k on network of star dunes. a. k = 0.25; b. k = 0.50; c. k = 0.90.

Figure 9. Simulation of complex dune formations.

Figure 10. Sample simulation resulting in a pier-like dune.

The first logical step after simulation of general Brownian motion (Fig. 2) is to put a uniform wind field onto the grid (Fig. 4). In this case the ridges start to form, aligned with the wind direction. This may look as a heresy to an experimental geomorphologyst, for the dunes are well known to align perpendicularly to the wind direction (see e.g. (Livingstone, 1988; Wiggs, 2001). This situation stems from the model’s abstraction level, because the wind field does not change as the landscape is being built up, so, even in the lee-side of the grown dune the wind vector remains the same.

When the wind vector begins to vary from cell to cell one gets exceptionally distinct picture of dunes formations evolving. The model performed well in simulations of linear (Fig. 5), barchan (Fig. 6), crescentic (Fig. 7), star (Fig. 8) types of dynes, as well as mixed landscapes (Fig. 9).

The star dunes were the easiest to reproduce, whereas barchans did not look quite perfect. These results are consistent with earlier findings (Wasson and Hyde, 1983; Wiggs, 2001), stating that barchans forms when there is a shortage in sand supply and winds are unidirectional, but star dunes grow up if there is plethora of sand, yet winds are mostly variable. The bi-directional wind hypothesis (Livingstone, 1988, 1989, 1993; Wang et al., 2003; Wiggs, 2001) worked well in the proposed model, consistently generating linear dunes (Fig. 5), in spite of this type of dunes also requires restricted supply of particles (Wasson and Hyde, 1983; Wiggs, 2001). The “bimodal wind” has been introduced by winds of uniform strength facing each other head-on at various angles. The angle between winds did not significantly influence the shape of dunes.

One peculiar case deserves attention, namely, the ability of the model to simulate landscapes that look artificial, like the pier-like dune shown in Fig. 10. This illustrates the ability of random processes to generate structures with ambiguous geometry, which can be useful for explanation of long-discussed observations of artificial formations on Mars (Van Flandern, 1998, 2002).

The effect of force constant, k

The model has but one parameter, force constant k, which encompass all the possible non-linearities of particles transfer. Therefore, it is necessary to demonstrate its influence on the landscapes.

Two cases are shown in Fig. 7 and 8 for crescentic and star dunes respectively. If one keeps the “sea level” equal for simulations with different k, it is clear that the lover constants make the dune boundaries more fuzzy and eroded, which is not surprising, because as k goes to null, the overall process becomes more and more Brownian. Also, the height of the dunes varies less with decreased k.

The current state of model clearly demonstrates, therefore, that even a little variations in the ratios between the transfer probabilities determines certain landscape form, whereas the exact amplitude of those probabilities has a little effect on it.

Numeric definition of the R, G and B color components for given wind vector field

The basement to use colour components for the wind vector definitions is the ability to build and use the complicated input easily perceivable by the humans. This is important when dealing with complex phenomenae, like turbulent flows. The obvious drawback here is a redundancy of representation when a wind vector is laid out as a sum of three colour component vectors. This redundancy must be worked around in some way.

From current arrangement of ,  and  vectors (Fig. 3) it follows that the wind vector  is defined as:


It is possible first to exclude the red component totally and solve the trivial system of two equations. The boundary conditions on colour components imply that they must stay in the range of 0–255. If, upon solution of those equations, the G and/or B components go out of that range, then it is necessary to raise gradually the R component until G and B get in the range.

Another drawback is that, in spite of the RGB components may lie anywhere between 0 and 255, the continuous radial coverage for resulting wind vector is possible only when actual wind force is not exceeding 220 units. So the actual field data should be scaled accordingly to this upper cutoff.


Anderson, R.S., 1988. Simulation of Eolian Saltation. Science 241 (4867), 820-823.

Basu, S., Richardson, M.I., Wilson, R.J., 2004. Simulation of the Martian dust cycle with the GFDL Mars GCM. J. Geophys. Res. 109 (E11), 1006.

Berry, R.S., Rice, S.A., Ross, J., 2000. Physical Chemistry, 2nd ed. Oxford University Press, New-York, 1064 pp.

Best, J., Kostaschuk, R., 2002. An experimental study of turbulent flow over a low-angle dune. J. Geophys. Res. 107 (C9), 3135.

Bishop, S.R., Momiji, H., 2002. Just Deserts. Mathematics Today 38 (4), 111-113.

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

Bitzer, K., Harbaugh, J.W., 1987. DEPOSIM: A Macintosh computer model for two-dimensional simulation of transport, deposition, erosion, and compaction of clastic sediments. Computers & Geosciences 13 (6), 611-637.

Borodin, A.N., Salminen, P., 2002. Handbook of Brownian Motion - Facts and Formulae, 2nd ed. Birkhauser Verlag, Basel, 672 pp.

Bourke, M.C., Bullard, J.E., Barnouin-Jha, O.S., 2004. Aeolian sediment transport pathways and aerodynamics troughs on Mars. J. Geophys. Res. 109 (E7), 7005.

Crave, A., Davy, P., 2001. A stochastic "precipiton" model for simulating erosion/sedimentation dynamics. Computers & Geosciences 27 (7), 815-827.

Crisp, J.A., Adler, M., Matijevic, J.R., Squyres, S.W., Arvidson, R.E., Kass, D.M., 2003. Mars Exploration Rover mission. J. Geophys. Res. 108 (E12), 8061.

Fenton, L.K., Banfield, J.L., Ward, A.W., 2003. Aeolian processes in Proctor Crater on Mars: Sedimentary history as analyzed from multiple data sets. J. Geophys. Res. 108 (E12), 5129.

Ferri, F., Smith, P.H., Lemmon, M., Renno, N.O., 2003. Dust devils as observed by Mars Pathfinder. J. Geophys. Res. 108 (E12), 5133.

Greeley, R., Arvidson, R.E., Elachi, C., Geriger, M.A., Plaut, J.J., Saunders, R.S., Schubert, G., Stofan, E.R., Thouvenot, E.J.P., Wall, S.D., Weitz, C.M., 1992. Aeolian features on Venus: preliminary Magellan Results. J. Geophys. Res. 97 (E8), 13319-13345.

Iversen, J.D., Rasmussen, K.D., 1999. The effect of wind speed and bed slope on sand transport. Sedimentology 46 (4), 723-731.

Koziy, L., Maderich, V., Margvelashvili, N., Zheleznyak, M., 1998. Three-dimensional mode of radionuclide dispersion in estuaries and shelf seas. Environmental Modelling & Software 13 (5-6), 413-420.

Livingstone, I., 1988. New models for the formation of linear sand dunes. Geography 73 (319), 105-115.

Livingstone, I., 1989. Monitoring surface change on a Namib linear dune. Earth Surface Processes and Landforms 14 (4), 317-332.

Livingstone, I., 1993. A decade of surface change on a Namib linear dune. Earth Surface Processes and Landforms 18 (7), 661-664.

Matsumoto, M., Nishimura, T., 1998. Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans. on Modeling and Computer Simulation 8 (1), 3-30.

Mihailovic, D.T., Koci, I., Lalic, B., Arsenc, B., Radlovic, D., Balaz, J., 2001. The main features of BAHUS — biometeorological system for messages on the occurrence of diseases in fruits and vines. Environmental Modelling & Software 16 (8), 691-696.

Nishimori, H., Yamasaki, M., Andersen, K.H., 1998. A simple model for the various pattern dynamics of dunes. International Journal of Modern Physics B 12 (3), 257-272.

Parsons, D.R., Walker, I.J., Wiggs, G.F.S., 2004a. Numerical modelling of flow structures over idealized transverse aeolian dunes of varying geometry. Geomorphology 59 (1-4), 149-164.

Parsons, D.R., Wiggs, G.F.S., Walker, I.J., Ferguson, R.I., Garvey, B.G., 2004b. Numerical modelling of airflow over an idealised transverse dune. Environmental Modelling & Software 19 (2), 153-162.

Raupach, M.R., Lu, H., 2004. Representation of land-surface processes in aeolian transport models. Environmental Modelling & Software 19 (2), 93-112.

Sato, T., 2004. The Current Status of the Earth Simulator. J. of the Earth Simulator 1 (1), 6-7.

Sirca, A., Rajar, R., Harris, R.C., Horvat, M., 1999. Mercury transport and fate in the Gulf of Trieste (Northern Adriatic)—a two-dimensional modelling approach. Environmental Modelling & Software 14 (6), 645-655.

Sorensen, M., 2004. On the rate of aeolian sand transport. Geomorphology 59 (1-4), 53-62.

Spahn, F., Krivov, A.V., Sremcevic, M., Schwarz, U., Kurths, J., 2003. Stochastic forces in circumplanetary dust dynamics. J. Geophys. Res. 108 (E4), 5021.

Toigo, A.D., Richardson, M.I., Ewald, S.P., Gierasch, P.J., 2003. Numerical simulation of Martian dust devils. J. Geophys. Res. 108 (E6), 5047.

Van Flandern, T., 1998. The "Face on Mars" at Cydonia: Natural or Artificial? Bulletin of the American Astronomical Society 30, 1453.

Van Flandern, T., 2002. Artificial structures on Mars. Bulletin of the American Astronomical Society 34, 784.

Wang, X., Dong, Z., Qu, J., Zhang, J., Zhao, A., 2003. Dynamic processes of a simple linear dune - a study in the Taklimakan Sand Sea, China. Geomorphology 52 (3-4), 233-241.

Warren, A., Allison, D., 1998. The palaeoenvironmental significance of dune size hierarchies. Palaeography, Palaeoclimatology, Palaeoecology 137 (3-4), 289-303.

Wasson, R.J., Hyde, R., 1983. Factors determining desert dune type. Nature 304 (5924), 337-339.

Werner, B.T., 1995. Eolial dunes: Computer simulations and attractor interpretation. Geology 23 (12), 1107-1110.

Wiggs, G.F.S., 2001. Desert dune processes and dynamics. Progress in Physical Geography 25 (1), 53-79.

Zhang, W., Qu, J., Dong, Z., Li, X., Wang, W., 2000. The airflow field and dynamic processes of pyramid dynes. J. of Arid Environments 45 (4), 357-368.



Please send your questions to: