Abstract
A large number of numerical models have been developed to simulate the physical processes involved in saltation, and, recently to investigate the interaction between soil vegetation cover and aeolian transport. These models are generally constrained to saltation of monodisperse particles while natural saltation occurs over mixed soils. We present a three-dimensional numerical model of steady-state saltation that can simulate aeolian erosion, transport and deposition for unvegetated mixed soils. Our model simulates the motion of saltating particles using a cellular automata algorithm. A simple set of rules is used and takes into account an erosion formula, a transport model, a wind exposition function, and an avalanching process. The model is coupled to the stratigraphic forward model Sedsim that accounts for a larger number of geological processes. The numerical model predicts a wide range of typical dune shapes, which have qualitative correspondence to real systems. The model reproduces the internal structure and composition of the resulting aeolian deposits. It shows the complex formation of dune systems with cross-bedding strata development, bounding surfaces overlaid by fine sediment and inverse grading deposits. We aim to use it to simulate the complex interactions between different sediment transport processes and their resulting geological morphologies.
1. Introduction
Aeolian sediment transport is an essential process with strong implications for the science community. In geomorphology, landscape development by winds presents a rich mixture of complex systems at various spatial and temporal scales [1]. In geology, evidence of wind processes eroding, transporting, and depositing sediment exists throughout the geological record primarily as vast accumulation of wind-blow sand and loess [2]. In areology, aeolian transport shapes the Martian landscape and dust aerosols are of major importance to the Martian climate [3]. In ecology, soil dust emitted into the atmosphere plays a major role in many earth system processes [4, 5], including by providing limiting micronutrients such as iron and phosphorus to a variety of ecosystems [6]. In climatology, wind driven sediment transport in the atmosphere causes scattering and absorbing of both short wave and long-wave radiation, enhancing melting of snow packs and glaciers upon deposition (Painter et al. [7]), and is possibly affecting hurricane formation in the Atlantic Ocean [8].
The transport of soil particles by wind can be separated into several physical regimes: long-term suspension (<20 mm diameter), short-term suspension (20–70 mm), saltation (70–500 mm), and creep (>500 mm) [10]. Saltation is arguably the most important physical regime, because it occurs at the lowest wind speeds and causes the other modes of wind blown sediment transport [10]. A large number of numerical models have been developed to simulate the physical processes involved in saltation [11–15], (McEwan and Willetts [16]), and [17–22]. Recently, advances in numerical modeling have facilitated investigation on interaction between soil vegetation cover and aeolian transport [23–29]. In evaluating the ecogeomorphic landscape evolution under different wind regimes, these models make a significant breakthrough in the understanding of wind and soils interactions. However, these models were generally constrained to saltation of monodisperse particles while natural saltation occurs over soils containing a range of particle sizes. Consequently, the internal structure and composition of the deposited sedimentary layers were lacking despite a consensus on their geological and environmental consequences. In addition, to improve our understanding of Earth surface processes and resulting landforms, models need to be able to simulate the complex interactions between different sediment transports in various environments such as waves, winds, or rivers. To date, as far as we know, no model has been able to investigate these interactions.
In this paper, we present a comprehensive physically based numerical model of aeolian transport that can simulate saltation of soils consisting of particles of various sizes. The developed model is mainly based on the previous works of de Castro [23], Nishimori et al. [30], and Baas [26]. Then, the model has been integrated to the 3D stratigraphic forward model Sedsim that enables a variety of interdependent surface and basin-forming processes to be studied at both geological and engineering time scales. Sedsim consists of a series of separate modules, each designed to represent a different physical process affecting sediment deposition. The main fluid flow/sedimentation routine is linked to optional additional routines such as sea level change, subsidence, and wave transport. The background to, and operation of, Sedsim is described by Tetzlaff and Harbaugh [31]. The program was developed by Tetzlaff, with subsequent enhancements and additions by Martinez and Harbaugh [32], Wendebourg and Harbaugh [33], Dyt et al. [34], Tuttle et al. [35], Griffiths et al. [36], and Li et al. [37]. Since 1999, the carbonates, organics, storm, and slope failure (turbidite flow) modules have been developed and linked to Sedsim [38, 39]. In this paper, the aeolian transport theory is firstly summarized. Then, the aeolian module that has been added to Sedsim is described for unvegetated soils. Then, examples simulate the modifications of surfaces subject to aeolian transport. It illustrates the formation and evolution of typical dunes patterns depending on wind field and regime.
Finally, enhancements of the aeolian module are proposed to help aeolian deposits preservation in a geological context and focus on vegetated soils systems and water table level fluctuations.
2. Aeolian Transport Theory
2.1. Definition of Sand and Dust
The distinction between sand and dust is conventionally set at a particle size of 60 mm [1]. Nevertheless, in the model, soil particles will be considered as dust if they can be readily suspended in the atmosphere once airborne. The mean Lagrangian vertical velocity is , where is the friction velocity and is the von Karman constant [40]. The friction velocity can be calculated to sufficient accuracy from the wind speed measured at a specified wind reference height , using the logarithmic wind profile law stability corrections as given by Paulson [41] with the stability function expressed as with where is the zero-displacement height, is the aerodynamic roughness length of the surface, is the Monin-Obukhov length, and is sensible heat flux. Estimation of aerodynamic roughness length for different regions over the globe can be found in the literature [42]. When , upward dispersion of the particles dominates gravitational settling, so passive suspension is a good approximation for the particle motions; when , gravity dominates over atmospheric dispersion. In our model, we define dust particles as those with fall velocities , or diameters , where is the solution of For a spherical particle of diameter , is given by [43] where and are the particle and air densities, respectively, is the drag coefficient of a sphere, and is the Reynolds number at terminal velocity. For very small particles ( μm), and while for large particles ( μm), is constant at about 0.48 and . At intermediate values, this equation can be solved for by successive approximation, as in [43]. The choice of means that dust particles, once ejected from the surface, are likely to remain suspended for some time; this is consistent with our neglect of dust redeposition. The ratio of is found to be in the range of 0.2–1.25.
The definition implies that the distinction between dust and sand is no longer fixed but depends on the turbulence, characterised by . In a complementary way, sand in the model is defined as particles with and , where is the diameter of the largest particles, which can be mobilised in given wind conditions. This diameter is the solution of , where is the threshold friction velocity.
2.2. Threshold Friction Velocity
In the model, we consider the prediction of , the threshold friction velocity for a dry, bare (unsheltered) bed of uniform soil particles of size . Based on [2] could be approximated by the semiempirical expression where is the gravitational acceleration, is the particle-to-air density ratio, is an empirical coefficient, and are empirical functions, and (with the kinematic viscosity of air) is the threshold particle Reynolds number. Coefficient and functions and are determined by fitting to wind tunnel measurements of .
We express using a quadratic form of the function which fits with experiment values of Fletcher [44, 45], Iversen and White [46] and Cleaver and Yates [47]:
3. Model for Unvegetated Soil
Our model is based on the previous work of Nishimori and Ouchi [48] and Nishimori and Tanaka [9, 25]. It consists of horizontally extended 2D lattice (Figure 1). Each cell of the lattice corresponds to the area of the ground sufficiently larger than individual sediment grains. At each cell (), at each time step , a continuous field variable is allocated to denote the average height of the noncohesive sediment surface within the cell. Therefore, the evolution of at one time step does not express the movement of individual sediment grains. It rather describes the resulting surface height change after the collective motion of many grains during the unit time period much larger than the time scale of individual grains dynamics [9].
The dynamical aeolian process developed is divided into two types. (i)The inertial or advection process, which describes the accumulated transport effect by saltation over the unit time scale and the unit space scale of our non-cohesive grained model. (ii)The frictional (or diffusion) process, which takes into account the fluctuation effects around the average motion due to erratic wind directions, irregularities of the soil surface. this erratic “Brownian” motion is modified by the slope effect.
In the frictional process, we assume a local conservation law of the amount of sediment in the space mesh. Therefore, the discretized conservation law of is defined as where is the mass of sediment coming into site at time step while is the one leaving from at n, is the mass density of the bulk of sediment, and is the area of each site. Both the saltation mass inflow/outflow and the creep mass inflow/outflow contribute to the total mass inflow/outflow .
3.1. Transport Model
According to Almeida et al. [49] universal equation for saltation length can be obtained. After dimensional analysis and assumption made on the constants of proportionality, they arrived at the following expression, which is used here to predict the saltation trajectories for close to the threshold . The transport saltation length of grains leaving where is the characteristic length scale with where is the air viscosity ( kg/ms, gravity g = 9.81 m/s2, and air density = 1.225 ). represents the grain velocity necessary to escape from the sand [50].
Our aeolian transport model is based on the assumption that sediments are moving at the same speed as the wind velocity. Thus, we know for each time step the maximum distance covered by our eroded grains. If the transport saltation length is smaller than the cell length the sediment will stay in the cell where it was eroded. If not, we update the sediment trajectory according to the wind speed and direction in the transported cell. Using this wind-tracking method, we accurately follow the sediment transport over the simulated surface. This is really useful in case of a precise description of a wind regime over a region.
3.2. Formulations for Erosion
Owen [51] developed a theory for “transport-limited” saltation, applicable to a dry, bare surface of cohesionless, uniform sand particles: where is a coefficient of order 1 and is the air density. This equation has been shown to describe observations of saltation over loose sand surfaces. There is some uncertainty in the value of the coefficient , which is not necessarily constant and can depend on the ratio . In Owen's original formulation, , was found empirically to be Our aeolian model uses (10) with the formulation proposed in (11).
For each sediment size, we calculate the maximum saltation flux in kg/ms. Then using the simulation time step , we obtain the maximum saltation mass leaving for the considered sediment type : Thus according to the discretized conservation law (8), the maximum thickness of sediment d leaving cell during due to saltation is given by In our model, the bed is divided in several sedimentary layers. Each of these layers embeds several grain sizes and is also characterized by hardness, porosity, and permeability, which vary with time.
Considering that erosion occurs on cell , we proceed layer by layer. For each cell, we make the assumption that erosion stops if (i)the transport capacity of the wind is reached; depending on the wind speed, the surface sediment composition and the simulation time scale, one or several layers have been eroded partially or entirely and for each grain size d, ; (ii)in a considered layer, one or several sediment types cannot be removed by the wind. Thus, erosion stops even if erosion of other grain sizes is still possible in underlying layers and transport capacity by the wind for these sediments is still positive.
Thus, the effective saltation thickness for this considered cell is then given by: where is the thickness of sediment type eroded from the sedimentary layer due to saltation.
3.3. Wind Exposition and Avalanching
Besides the main sand transport process, two constraints are simulated in the model: shadow zones and avalanching.
A shadow zone is the area in the lee side of relief where wind flow has been slowed down sufficiently to suppress any further transport of sand. In the model, this is represented by a shelter zone downwind, in which all-passing sediments are trapped and erosion never occurs.
Avalanching is simulated to maintain the angle of repose of loosely packed sand, which is usually an angle of 33 to the horizontal. However, this angle must be put in perspective with the size of the simulation cells. It appears in Sedsim that for commonly used cell size (around 50100 m) the angle of repose must be defined much smaller to be able to reproduced realistic features. The value is usually set around 5 in our simulations. After removal or addition of a slab of sediment (erosion or deposition), the model assesses whether this angle of repose is exceeded and if so, moves neighboring slabs down the steepest slope (in the case of erosion) to reinforce this angle or moves the newly deposited slab down the steepest slope until it reaches a grid cell where it does not violate the angle of repose. Avalanching is the only means of sideways sand transport relative to the transport trajectory.
4. Examples
Previous section has defined the main frame of our aeolian module. We are now going to illustrate our model with two examples. As input, two modules (aeolian and expert aeolian) can be used depending on the level of understanding and the amount of data available for the simulated region.
The aeolian module is the simplest and required only an idea of the sedimentary sources at the boundaries and the wind regime: time and wind speed values on () coordinates. For a given time, a unique wind regime is defined on the entire region.
The expert aeolian module uses the logarithmic wind profile law stability corrections given by Paulson [41]. In addition, users could launch several wind regime maps representing an accurate description of wind circulation. Each of these maps embeds wind speed coordinate values at each grid point.
4.1. Monodisperse Particles Evolution due to Changing Wind Field
In this example, we use the aeolian module to simulate different wind fields. Only one type of sediment (diameter of 0.1 mm) is considered. The surface is an initial flat and square 5 km side area meshed with 50 m square cells. The wind regime varies during the 4 years simulation and black arrow in Figure 3 summarized the circulation pattern. The wind speed is set to 20 km/h during calm period and to 50 km/h during strong wind activity. Depending on the wind blowing direction, sources are activated along the south, west, and east side of the simulated area and consisted in 0.4 and 1 m thick sand slab for calm and strong activity, respectively. At each time step the thickness of sediment transported by to the wind from the source is replaced. Thus, the sand slab at the source remains constant over the entire simulation. Basement hardness is set sufficiently high and does not suffer from wind erosion.
Landscape construction and evolution in response to wind regime is simulated and presented in Figure 2. Our aeolian module reproduces most of the typical natural features such as transversal, barchan or seif dunes (Figure 3). The resulting morphological complexity shows how dune systems responds to wind strength and direction variations. It appears that there is not any specific wavelength in the dune field. Each dune has its own movement depending on the neighbourhood. Dune speed tends to accelerate when the downwind surface is flat but slow down close to another downwind dune and could finally catch up with this other one creating a larger and thicker dune. Dune heights over the simulation are ranging between 4 and 8 m. Maximum heights are reached when two dunes merged to form a wider one.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
According to our simulation, there is no real accumulation of any aeolian deposits in the region. A sand dune migrates quickly and except some cases of merging dunes there is not a lot of mixing over time between the different sand slabs. This is mainly due to the small scale of the simulated region in comparison to the grid cell size.
This test shows that when a sufficiently large amount of sand is supplied under strong winds, transverse dunes appear in the system. These dunes can connect to each other and grow larger as they move downwind (Figures 2(a) and 3(a)). The crests of the dunes extend perpendicular to the direction of wind.
When wind field shifts towards the NW or NE direction and decreases in strength (Figures 2(b) and 2(d)) barchan dunes are formed (Figures 3(c) and 3(b)). In more details, small parabolic dunes are formed just to the lee side of the sand source and quickly developed into barchan dunes. Barchan dunes arms length varies depending on the wind force, they extend further under stronger wind. Depending on the wind strength they will connect to each other and grow larger as they move to form transverse dunes (Figure 3(c)).
Under stronger winds, previous deposited sand is quickly transported out of the simulation and a new landscape appears with seif dunes formation. They develop along the wind direction and reach up to 8 m for the highest ones. In our simulation, the longer ones extend over 1.5 km (Figure 3(d)). They are sharpcrested and develop from barchans when the change of wind direction occurs. The new wind direction leads to the development of a new wing and the overdevelopment of one of the original wings.
4.2. Polydisperse Particles Evolution due to Changing Wind Field
We use the expert aeolian module with four types of sediment (diameters: 0.11, 0.09, 0.07, and 0.05 mm). The surface is an initial flat and 6.2 × 8.70 km side area meshed with 50 m square cells. Figure 4 shows the four types of wind direction that are used in the simulation, and Table 1 summarizes the organization of these wind fields during the modeled four years.
(a)
(b)
(c)
(d)
The maximum wind speed corresponding to storm winds is equal to 50 km/h and the average yearly wind is around 20 km/h. Northerly, northwesterly, and northeasterly winds are the three main active directions of the wind climate in the region. Sources are activated along the south side of the simulated area and consisted in 0.4 and 1 m sediment slab thickness for calm and strong activity, respectively. The sand slab is made up with 20, 30, 30, and 20% of finest, fine, medium and coarse sediments. Here also, basement hardness is set sufficiently high and does not suffer from wind erosion.
Figure 5 presents surface evolution and grain-size repartition at 1, 2, 3, and 4 years. It shows the construction of complex morphological features with evolution of the system from transverse to barchan dunes. During the first year, transverse dunes (less than 3 m height) form in the vicinity of the source. During their transgression these dunes develop in more complex transverse dune systems mainly due to both NW and NE wind flow. As shown in top of Figure 5, dune heights could reach up to 12 m during this year. During the second year, transverse dunes are still moving over the system and some barchan dunes form with dissymmetrical arms due to wind direction. Some of these barchan dunes connect to each other and form transverse dunes that could reach 18 m. During third year, a special dune feature with a thickness of about 25 m forms and results from the special pattern of the wind field in the region. The structure looks like a barchan dune with two long and thick arms pointing down the main northerly wind. However, the slipface is not parabolic due to a remaining transverse dune structure that is not completely connected to the other dune. The last year shows the evolution of the previously described dune, which now has an approximate height of 36 m. Some small transversal dunes (1 to 3 m thick) develop on the top of the previous structure forming a complex mixing system.
(a)
(b)
Figure 6 illustrates the internal structure of the resulting topography after 4 years. It focuses on the sedimentary layer and composition of this final dune, which has a height of approximately 36 m. The region is exposed to wind from a wide range of directions (Figure 4). Consequently, the dune that is formed records variable directions of sand movement. The resulting structure is close to star dune ones, which are characterized by slipfaces pointing in all directions and by arms spreading in opposing directions. Due to this typical morphology, star dunes often make little overall progress in terms of transgressing a surface. As shown in left images, coarsest grains are deposited near the upper part of the lee slopes. This result is consistent with the process of shear sorting during the transport, which causes coarser grains to rise to the surface, resulting in inverse grading. Inverse grading in some grainflow and wind deposits is so distinctive that it has been referred to as pin-stripe bedding [52]. With grainflow deposits, however, it is not always clear whether the fine material at the base represents infiltrated sediment during avalanching or a layer of wind transported laminae deposited between avalanches. The simulated organization is in good agreement with outcrop and seismic observations from various regions around the globe [53–55]. Right images from Figure 6 shows the development and overlapping of several dunes sequences leading to the construction of heterogeneous and complex deposits. The dune development involves lateral expansion in both a westerly and easterly direction. Figure 6 clearly shows cross-bedding dipping in both directions. While the main dune progression is dominantly toward the north. Bounding surfaces are preserved and correspond to interdune migration surfaces. The figure shows that fine sediment deposits, which have smoother dip cross-beds, overlay most of these bounding surfaces.
(a)
(b)
(c)
(d)
(e)
(f)
4.3. Comparisons with Other Studies
The main contribution of our work to aeolian sediment transport modelling is the ability to simulate polydisperse particles saltation in three dimensions and the resulting complexity of dune sedimentary facies at regional scale (see Figure 6). Our model is currently successful in generating relatively realistic dune landscapes and it appears that a set of well-chosen rules is of enough accuracy to understand the fundamental processes resulting from polydisperse particles saltation at regional scale level. It seems that there is no need to incorporate all the complex details present on a smaller scale. For example, secondary wind flow patterns induced by relief are apparently not essential to the development of realistic dune topography. As suggested by Baas [26], this type of approach provides a connection between descriptions of a more geologic or physiographic nature [19] and small-scale deterministic processes [20, 21]. The level of description of the physics used in the model is far from being perfect but seems sufficient to investigate the impact of aeolian dynamics on a landscape-scale level. Other works such as the one from Kroy et al. [18] proposed a two-dimensional description of the shear stress based on turbulent boundary layer calculations that can be applied to simulate large-scale desert topographies. Their results help to better understand the mechanisms responsible for longitudinal shape and aspect ratio of dunes and heaps. However, using these results for 3D problems is still under question when dealing with nontransverse dune fields.
In addition, the aeolian module, presented in this paper, is coupled to the stratigraphic forward model Sedsim. It may allow both landscape reconstruction and heterogeneities characterization of a complete geological history for regions that undergo multisedimentary processes due to various environments (marine, coastal, fluvial, and aeolian).
5. Discussion
The results previously described clearly show the potential of this approach for simulating different realistic dune patterns under the influence of wind fields. The great contrast in appearance between the landscapes of Figures 2 and 5, for example, shows the large impact of wind regime (direction and strength) and sediment supply on the developing morphology. It shows how a change in parameters results in a fundamentally different landscape. The driving force in the system is the transport of sand by wind (it corresponds to the perturbation energy given to the initial stable system). This energy is transformed in potential energy created by the elevation differences due to the formation of dune morphology. Individual sediment slab trajectories are unpredictable due to the complex interactions between the erosion and deposition rules for the different grain sizes and the morphology in the landscape.
The created potential energy is dissipated through the process of avalanching which turns the potential energy into kinetic energy and finally into frictional heat after coming to rest at a lower elevation. Thus, dunes could be considered as dissipative structures formed during a disequilibrium situation and which are used, in nature, to absorb this perturbation energy in order to find a new equilibrium profile. This process leads to the formation of complex structures, which are evolving rapidly in a geological point of view. In this context, further developments need to be undertook to simulate the way aeolian deposits could be preserved through time. Two critical features discussed below seem to play a critical role in the preservation of aeolian landforms: vegetation covers and water table level fluctuation.
There are strong relationships between ecological and geomorphic processes in the development and preservation of aeolian landforms [56]. In order to bring vegetation into the model, an additional variable describing the vegetation effectiveness at each cell site need to be added. This term refers to Baas [26] definition and can be interpreted as a coverage density or a frontal area index and has a value between 0 and 100%. The vegetation effectiveness is supposed to be a function of sedimentation balance and is defined with growth function curves for typical vegetation types. The vegetation effectiveness affects the erosion and deposition process, but not the intermediate transport trajectory. This simple approach is usually used to simulate the influence of vegetation on the threshold shear velocity required to initiate and sustain sand transport rates [57–59]. These functions are deliberate simplifications that can only be described in terms of general characteristics, such as steepness (i.e., rate of response to burial or erosion) and the location of the peak (the optimum growth with respect to the erosion/deposition balance).
Water table level fluctuation is another fundamental mechanism that enables the accumulation and long-term preservation of wet aeolian systems [60]. It may be either absolute, resulting from a climatic shift to wetter conditions, or relative, whereby the sediment accumulation subsides through a constant water-table level. A shallow water table, either at the depositional surface or within its capillary fringe, ensures that moisture at the interdune surface is at least sufficient to raise erosional threshold values to the point at which the surface is largely protected from deflation, thus increasing preservation potential [61]. Wet aeolian systems characterized by migrating dunes and a rising water table will therefore preserve climbing sets of dune-interdune strata. The angle-of-climb of the aeolian bedforms and adjoining interdune flats is controlled by both the time-averaged rate of water-table rise and the rate of downwind migration of the bedforms. Modeling these fluctuations will need further enhancements of the code and rely in a good understanding of the physical process. Our proposition consists in defining water level curves as the ones used for sea-level variations. In a coding point of view, this solution is relatively simple to carry out as our model Sedsim already embeds a sea-level fluctuations module. However, on the user side, these curves need to be defined precisely for each particular study depending on the geological and climatic context and on the data sets available for the considered study of interest (such as wells analysis or seismic profiles).
6. Conclusion
This paper presents a three-dimensional multigrain aeolian sediment transport model coupled to other transport processes. Sedsim Aeolian module is capable of simulating realistic-looking dunes development at various space and time scales and could interact with previous modules such as fluvial fluid elements or wave circulation. Simulated dunes fields could be of different types transverse, longitudinal, seif or barchan. These complex landforms develop as a result of interactions between wind climate (function of both speed and direction), grain size, and mixed sediment transport processes, modeled by simple, local rules in a cellular automaton algorithm. In the future the model will be used to improve our global understanding of earth surface systems processes interactions and to investigate both environmental (climate change, desertification, and coral threats) and geological (aeolian/fluvial processes interaction, aeolian sands reservoir rock formation) issues.