Toward a large-scale particle-based parallel simulator of Aeolian sand trans- port, including a model for mobile sand availability

We develop a numerical tool for particle-based simulations of Aeolian sand transport. Our model combines a Discrete-Element-Method for the sand particles with an efficient hydrodynamic description of the average turbulent horizontal wind velocity field over the granular bed, which has been developed in previous work and accounts for the two-way coupling of the granular and fluid phases. However, here we implement our model within the open source library LAMMPS for granular massively parallel simulations and incorporate a new grid coarsening scheme for the wind model. We show that our model quantitatively reproduces observed values of the steady-state (saturated) sand flux under various flow conditions. Furthermore, we model different conditions of mobile sand availability and find a strong dependence of the sand flux on this availability.


Introduction
Wind-blown (Aeolian) sand plays a vital role for a broad range of geophysical processes, including ripple and dune migration, rock abrasion and the emission of dust aerosols [1,2]. These processes are affected by the still poorly understood interactions of wind-blown sand particles with different types of soil, which can be fully or partially erodible. Such interactions and their effects on the sand flux are difficult to represent in theoretical models [3][4][5].
Therefore, there is a need to have an accurate, but also fast numerical model for particle-based simulations of Aeolian sand transport. Previous work has introduced a Discrete-Element-Method (DEM) for wind-blown sand, by two-way coupling the granular phase with an efficient hydrodynamic description of the average turbulent horizontal wind velocity over the surface [6,7]. However, here we present a model for massively parallel simulations, which incorporates an improved wind model and leads to enhanced computation performance by taking the vertical particle concentration profile in the Aeolian transport layer into account. Moreover, we present a method for DEM simulations under different soil conditions, ranging from a non-erodible surface to a fully mobile sand bed.

Numerical model
We develop our numerical model using the open source DEM library LAMMPS [8], which incorporates an environment for massively parallel simulations of granular systems (MPI implementation). However, this library is extended here to account for a hydrodynamic model and its coupling to the granular phase (see next section). * e-mail: skamath@uni-koeln.de * * e-mail: eric.parteli@uni-due.de  The sand grains are assumed to be spherical with physical attributes displayed in table 1. The contact force between two colliding particles encodes a normal and a tangential component, each comprising an elastic and a dissipative term (linear spring-dashpot model with parameters as in Ref. [9]), while the magnitude of the tangential force is limited through the Coulomb criterion (F t ≤ µF n , with µ denoting the coefficient of friction).
To produce the granular bed, we generate a loose packing (< 20% solid fraction) of about 65k sand grains of mean diameter d m = 200 µm and 20% polydispersity ( Fig. 1), and let the grains settle under the action of gravity. The dimensions of the simulation domain in units of d m read L x = 500 × L y = 8 × L z = 1000, with x and y denoting the along-wind and cross-wind directions, respectively, and z standing for the vertical direction. Boundary conditions are periodic in x and y and a frictional (reflective) horizontal wall is considered at the bottom (top).
Moreover, the granular bed is subjected to a fully turbulent atmospheric boundary layer wind flow. To model such flow, the vertical domain (z) is divided into horizontal layers of thickness (width in the vertical direction) equal to d m , each with velocity u(z) given by where u * is the wind shear velocity, κ=0.4 s the von Kármán constant, z 0 ≈ d m /30 is the bed roughness and h 0 is the bed height. To initiate transport, we produce a few grain-bed impacts that induce a granular splash thus entraining grains into the transport layer [7]. However, once transport begins and the grains accelerate, they extract momentum from the wind, so that the wind profile (Eq. (1)) must be updated during the simulation. The modified wind profile under consideration of the grain-borne shear stress τ p (z) is computed with the algorithm introduced in Ref. [7]. The corresponding workflow is displayed in Fig. 2.

Model validation
The steady-state (or saturated) sand flux q sat is computed as a function of the wind shear velocity and compared with experimental results [10]. The circles in Fig. 3 denote our numerical predictions for this flux, rescaled aŝ for different values of the Shields number θ, which is defined as The range of θ in Fig. 3 corresponds to 0.18 ≤ u * /(m s −1 ) ≤ 0.4, thus including values of u * characteristic of Aeolian dune areas [11]. We see from Fig. 3 that our numerical predictions agree quantitatively with windtunnel observations [10] (stars). Moreover, the dashed line in Fig. 3 denotes the best linear fit to our simulation results, i.e.,Q = 23θ − 0.12. From the x−intersect of this dashed line we find the minimal threshold Shields number for sustained transport, θ it ≈ 0.0064, which corresponds to the minimal threshold shear velocity u * it = 0.165 m/s. This value agrees well with the Bagnold's prediction that u * it is approximately 80% of the minimal shear velocity required to initiate saltation, u * ft ≈ 0.1 sgd m [1], which is around 0.206 m/s considering the mean particle diameter d m = 200 µm used in our simulations.

Model speedup scheme
To speed up our numerical simulations, we perform parallel computations with domain decomposition along the wind direction and a grid coarsening scheme for the vertical direction, which takes advantage of the boundary conditions inherent to Aeolian transport, as specified next.

Parallelization
To apply domain decomposition along the wind direction (x-domain), we have extended the MPI (Message Passing Interface) implementation of LAMMPS for the granular phase to allow for parallel computations including our wind and granular-fluid coupling schemes. Figure 4a shows that the vertical profile of the particle concentration does not change significantly with the number of processors concurring to the computation. The performance of the simulations as a function of the number of processors is discussed in the next section.

Grid coarsening
Our simulations reproduce the experimental observations [10,12] that, well above the ground, the volume particle concentration decreases roughly exponentially with the height. Moreover, as also found experimentally [12], we also observe a deviation from this exponential decrease near the bed (z 50 d m ). We can thus further improve the speedup of our simulations by increasing the thickness of the horizontal layers of the wind scheme with the height above the ground. Here, the vertical position z n of the n−th layer above the ground (n = 0, 1, . . .) is computed with the equation with A = hh/2 and B = (1 − hh/2), where the parameter hh dictates how fast the layer size increases with the height above the ground. The effect of hh on the simulation results is shown in Fig. 4b. Moreover, Fig. 5 summarizes the effect of hh and the number of processors on the computation performance.
As we can see in Fig. 5, for the system considered here (total simulation time = 10 s for all case scenarios) and taking hh = 0, we obtain a speedup of ≈ 10 by increasing the number of processors from 1 to 16. Moreover, we see in Fig. 5 that, by further considering the grid coarsening scheme along the vertical direction (hh > 0), we can further increase the performance of our simulations. Indeed, an increase in hh reduces the time (t u ) spent on the iterative computations of the wind field, specified in Fig. 2, by reducing the number of horizontal layers in Eq. (4). In simulations using 1 processor, the values of t u /t g for  However, further work is required to shed light on the simulation performance with larger systems including millions of particles. Furthermore, other grid coarsening schemes shall be considered in future work. For instance, we shall consider an exponential increase of the grid size with the height above the ground (z) as an alternative to Eq. (4), given the nearly exponential decrease of the particle concentration with z (Fig. 4).

Model for mobile sand availability
One of such boundary conditions is the amount of mobile sand available for transport, which affects dune shapes, inter-dune flux, and has still poorly understood effect on dust aerosol emission processes. We have introduced a module ("fix" feature) into LAMMPS for simulating rough beds through "freezing" the lowermost bed particle layers. The equations for the contact forces between the mobile and frozen particles (yellow and red particles in Fig. 6, respectively) are solved by considering that the latter have infinite mass and radius [13,14]. Moreover, interaction forces between frozen particles are not computed. In our model, the thickness of the mobile sand bed (h mob in units of d m ) is set, thus, by adjusting the number of frozen layers underneath. We find that the sand flux is strongly affected by h mob (Fig. 7). Indeed, it was shown in wind tunnel experiments of Aeolian sand transport that the sand flux over a rigid bed is much higher than that over erodible beds [15]. This behavior can be understood by noting that, over rigid beds, the grains reach greater heights as they experience a higher effective coefficient of restitution upon collision with the bed. These grains then achieve higher speeds resulting in higher flux values. However, using our numerical simulations, we can now investigate sand availability conditions that are intermediate between fully rigid and fully mobile beds, by changing h mob . Fig. 7 reveals an anomaly for transitional values of the mobile thickness layer, where the flux has a local minimum at low h mob . This phenomenon shall be now studied in depth through systematic simulations as a function of the flow conditions, using our DEM tool.

Conclusion and outlook
We developed a numerical tool for particle-based massively parallel simulations of Aeolian sand transport, which two-way couples the granular phase, computed with a Discrete-Element-Method, with the fluid phase, described through a hydrodynamic model for the vertical profile of the average turbulent horizontal wind velocity.
Our simulation predictions for the sand flux over an erodible sand bed as a function of the Shields number, as well as for the minimal threshold shear velocity for sustained transport, agree quantitatively with experimental observations. Moreover, simulations using our model for rigid beds led to much larger values of sand flux, a result consistent with wind tunnel observations [15]. However, for intermediate conditions of sand availability, we discovered an anomaly (local minimum) in the flux dependence on this availability, which shall be elucidated in future work.
Moreover, we investigated the speedup of our simulations upon application of parallel computations with a suitable domain decomposition and a new grid coarsening scheme for the wind model, which takes advantage of the vertical profile of the particle concentration in the Aeolian transport layer. However, Eq. (4) shall be replaced in future computations by an exponential increase for the grid size with the height above the ground, in view of the results displayed in Fig. 4. Based on the results presented here, our study shall be now extended to develop optimization strategies for large-scale systems with millions of particles. This is important, for instance, to elucidate the bed evolution in more detail and the granular processes underlying three-dimensional bedform instabilities.
The future application of our model has the potential to improve soil parameterization schemes in dust and climate models, which rely on the accurate representation of Aeolian surface processes [2]. Moreover, future work includes geomorphic simulations including a vegetation cover and its effects on the wind field and the sediment transport.