3D Direct Numerical Simulation on the Emergence and Development of Aeolian Sand Ripples

A sand surface subjected to a continuous wind field exhibits a regular ripple surface. These aeolian sand ripples emerge and develop under the coupling effect between the wind field, bed surface topology, and sand particle transportation. Lots of theoretical and numerical models have been established to study the aeolian sand ripples since the last century, but none of them has the capability to directly reproduce the 3D long-term development of them. In this work, a novel numerical model with wind-blow sand and dynamic bedform is established. The emergence and long-term development of sand ripples can be obtained directly. The statistical results extracted from this model tally with those deduced from wind tunnel experiments and field observations. A simplified bed surface particle size description procedure is used in this model, which shows that the particle size distribution makes a very important contribution to sand ripples’ final steady state. This 3D bedform provides a more holistic view on the merging of small bumps before regular ripples’ formation. Analyzing the wind field results reveals an ignored development on the particle dynamic threshold during the bedform deformation.


INTRODUCTION
Aeolian sand ripples are commonly observed in the deserts of Earth and some other planets, which exhibit regular patterns with wavelengths that vary from centimeters to tens of meters and amplitudes that vary from millimeters to centimeters. Ripples' morphological characteristics are wrinkle-like stripes perpendicular to the wind direction and almost symmetrical in the transverse direction. General sand ripples emerge rapidly in a continuous wind field and keep growing until they reach a saturated state. Its development process is always separated into two stages which are called the linear stage and the nonlinear stage. Starting from the phenomenological description developed by Anderson [1], researchers found that the linear stage is the first stage of ripple formation, in which the most unstable mode of the ripples' amplitude grows exponentially. After the initiation of the ripple instability, nonlinear behavior such as the coarsening process takes place [9,28,34] and the ripple growth rate slows down. The formation of the sand ripple bedform has been considered as an important research subject in planetary geology since the last century on account of the reflection of local atmospheric conditions and granulometric property. However, due to the difficulty of direct observation, the causes of many distinct ripple properties remain as unsolved questions. Thus, a proper model on the development of ripple morphology is required.
To research sand ripples, an eligible model should be able to correctly reproduce the rule of aeolian sand movement because the emergence and development of aeolian sand ripples are currently found to have a strong relationship to the sand transport process. Commonly, sand particles are transported by the wind in three major ways: reptation (or creep), saltation, and suspension. The saltating particles are those going forward by bouncing on the bed and splashing other ejected particles. These particles' trajectories are high enough to regain the impacting energy loss from the wind and support them to travel a long distance. According to Bagnold's seminal work [6], the formation of sand ripples is only influenced by the characteristic saltation length. However, his theory was disproved by experimental results [33,36], which indicate that the initial ripple wavelength is much smaller than the particle saltation length. Anderson then made a very important contribution to the research on aeolian ripples by introducing the importance of the reptation process [1]. He defined reptating particles as those sand particles motivated by the saltating process containing very low energy which can only support them to have one small hop. After one hop, a reptating particle dies immediately without inciting any new particles into the air. Anderson argued that the ripple formation mechanism was entirely attributed to the contribution of reptating particles.
According to Anderson's assumption, many theoretical models on sand ripples' formation have been established [9,17,28,34,37,40]. These theoretical models describe the continuous wavy bed surface as partial differential equations, which provide convenient methods to deduce the theoretical or numerical solution of bed form morphology. For many years of development, lots of refinements and development have been implemented on this kind of continuous model, making it suitable for more and more complex situations. However, no matter how advanced it is, some arbitrary simplification remains unchanged. By omitting hopping particles' trajectories, these models simplified the description of particle dynamics and only paid attention to the initial/final location of moving particles. Saltating particles in these models do not take part in the ripple formation directly but are treated as external energy sources just bringing energies into the mass exchange system. What is more, this kind of model involves a lack of wind field information, which makes the research on the relationship between wind velocity and ripple morphology unachievable.
These drawbacks of theoretical models make people study the direct method, which contains sufficient mechanism of particles' motion and takes into account the applied energy from the fluid field. In recent years, direct numerical simulations with particle dynamics have been frequently used for solving complex windblowing sand problems. Many are inspired by the work of Anderson and Haff [4]. They introduced hydrodynamic equations into the model and tracked every particle in the computation domain. Particles' movement and the interaction between air flow and particles are carefully treated by using properly simplified momentum equations. So far, just one direct simulation work has been carried out systematically on the research of aeolian sand ripples [12]. The discrete element method (DEM) used in this work considers all the forces applied on dynamic and static particles at every moment. Because of this feature, this work suffers a great limitation on computational resources. It simulated a very short time period which is only sufficient for the initial ripple forms' emergence. The subsequent turning point between the linear stage and the nonlinear stage cannot be directly obtained from their results, which makes them unable to take a further look at the transition between these two stages. Meanwhile, ordinary sand ripples not only develop in a timescale of tens of minutes but also require tens of thousands of particles to form a complex three-dimensional bedform. For DEM simulation, the particle number requirement of a 3D model is hard to achieve. Therefore, the simulation domain in their work is quasi-2D, which only has the length of one particle diameter in the transverse direction. Omitting the influence from the third dimension makes results less convincing, and the model will fail to explain many important 3D properties on ripple topology.
Be aware that in this work, we introduce a numerical model which just traces moving particles in the air and gains the postcollision velocity components of particles by solving the momentum equations in connection with Coulomb's law of friction. Particles dropped into or that escaped from the sand bed are converted to the elevation deforming of the local bedform. It greatly reduced the computational costs, making long-term simulation on ripple formation processes and 3D ripple morphology simulations achievable. Thus, using this model, we can perform better studies on the relationship between multifarious wind fields and ripples' development.

Particle Motion
Sand particles in this simulation have been assumed to be small spheres. Every particle in the air is tracked. Despite the wake influence of particle rolling, the equation of particle velocity components is simplified as follows: where m p is the particle mass, v is the particle velocity, f fp stands for the forces of hydrodynamical origin, g is the gravity acceleration, and f cp is the interparticle contact force. Naturally, f cp 0 if one particle is not in contact with another. For those particles within a contact pair, this value is not calculated directly here. The detailed dealing procedure on this term will be described in section 2.2.
The densities of the sand particle and the carrier fluid are ρ p and ρ f , respectively. As the density rate s ρ p /ρ f is large in this work and the particle diameter is much smaller than the Kolmogorov scale, the hydrodynamical force f fp here is dominated by the drag force of the fluid. For a particle with a certain diameter d, the drag force exerted on it is only influenced by the relative value between its velocity v and the fluid velocity u at its position: C d is the drag coefficient. We use the following approximation that C d ( C ∞ d + Re c p /Re p ) 2 [13], where Re p |u − v|d/] is Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 662389 the particle Reynolds number which is based on relative velocity. ] is the kinematical viscosity coefficient of the fluid. C ∞ d ≈ 0.5 is the drag coefficient of the grain in the turbulence limit (Re p → ∞), and Re c p ≈ 24 is the transitional particle Reynolds number.

Midair Collision and Surface Splash Function
Particles motivated from the bed are accelerated by the fluid alone in its projecting trajectory, and then some of them bombard the sand surface and entrain other particles. Momentum transfer frequently takes place via collisions among midair particles and the particles on the sandy bed during this process.
In this work, the collision event takes place while the centroid distance between a pair of particles is smaller than the sum of their radii. During the whole collision process, f cp should be calculated in no time to close (1). Since equations are solved in a dispersed manner while performing the numerical simulation, dv/dt in (1) becomes Δv/Δt. For the simplest linear spring model of contacting grains, the integrating during the contacting event requires the time step Δt to satisfy the limitation Δt Δt c ≪ T c ≪ π m p /k , in which T c is the total contact time and k is the spring stiffness [4]. Considering that k for quartz particles is very large, T c could then be a very small value. In DEM simulations, this is the most significant aspect influencing the computational cost. To avoid the influence from T c , in this work, the computational time step is only controlled by T f . T f 4sd 2 /(3]) is the relaxation time deduced from the rearranged form of (2), that is, f fp m p (u − v)/T f f (Re p ). For all the situations discussed in this work, T f is much larger than T c . Therefore, we define the time step as Δt Δt f < T f and Δt ≫ Δt c . By this definition, we can assume that the collision happens instantaneously in every time step.
Thus, in one iteration step, the computational process of particle movement is started by solving (2) using a Runge-Kutta method under the setting that f cp 0. Then, instead of solving f cp , the postcollision velocity is determined directly by the relative velocity between two colliding particles. Here, we omit the full calculation process and just provide the result for the post-collision velocity v 1 ′ on one of the particles after colliding [7]: where the variables with subscript 1 or 2 refer to the quantities of each particle in the collision pair and v 12 v 1 − v 2 is the relative velocity before collision; n is a unit vector from one particle center pointing to the center of the other one; ε and μ are microscopic restitution coefficients for the normal and tangential components, respectively; and q is a parameter which depends on the moment of inertia I 1,2 on the particle as shown below: Since the particles are assumed to be spheres, I 1,2 m 1,2 d 2 1,2 /10 and thus q 5/2.
The rotation of particles is not in the consideration of this work; then we have ω 1,2 0. (3) is reduced to the following: In terms of the effective restitution coefficients, we have the following: η m 1 /m 2 is the mess ratio of two colliding particles. The same procedure can be implemented to calculate v 2 ′ , and there is no need to repeat it here.
If one particle in the collision pair is at rest before the collision corresponding to the scene where an impactor hits the granular bed, (5) will be solved by setting v 12 v 1 or v 2 0. The variables with subscript 2 refer to quantities of particles on the bed hereafter. After considering all possible impacting statuses, a novel splash function is introduced by Lämmel [22], from which one can estimate the properties of the particles ejected from the bed surface, such as eject velocity, eject angle, and the number of ejected bed particles.
In this splash function, ejected particles are separated into two categories at first. The first one is called rebound particles, which are former injecting particles that then bounce from the bed after the collision event. It is quantified in terms of the mean restitution coefficient e and rebound angle θ 1 ′ for a given impact angle θ 1 , as shown below: where Particles entrained by the impactor are called ejected particles; the kinetic energy of them can be drawn from a log-normal distribution as shown below: where E 2 ′ is the kinetic energy of an ejected particle, E 1 m 1 v 2 1 /2 is the energy of the impactor, and we have the following: In here, E d2 m 2 gd 2 is the minimum transferred energy for a bed particle to be counted as ejecta. The number N of ejected bed grains is as follows: The ejection angles of all ejected low-energy particles are set to 90 + , and their initial positions are uniformly distributed near the impact point. As we are expecting a 3D model for particle entrainment upon an irregular bedform, there should be a particle trajectory component in the transverse direction (the y direction in this work). In others' work, a random distribution in the transverse direction of the rebound angle and the ejection angle was introduced into the model [21]. This artificial treatment will not be applied here because the uneven bed surface in this model introduces the y-direction particle velocity component automatically.
In (6) and (7), the restitution coefficients ε and μ characterize the relative normal motion-caused energy losses due to the deformation and the relative tangential motion-caused energy losses due to friction, respectively. ε can be treated as a constant parameter deduced from the particles' material property. For sand, ε 0.9. The assumption of μ is a little more complex. To derive the simple form expression of splash functions for particle-bed collision, we consider μ as a constant. It is natural to assume that the colliding particles roll past each other, that is, μ 0. For midair collision, as we know the exact velocity and location information of every collision pair, under the assumption of Coulomb friction, the tangential restitution coefficient μ can be deduced from the following [31]: where C f 0.4 [14] is the coefficient of friction, and v n and v t are relative velocities between two colliding particles in the normal and tangential directions, respectively.

Aerodynamic Entrainment
Static particles on the surface are not only motivated by impactors but also, for aeolian sand drifting, are directly entrained by the turbulence fluid. Aerodynamic entrainment plays a significant role in initiating a continuous sand flux. As the saltating process is gradually enhanced, the airborne shear stress on the sand surface (wall stress) τ fw decreases and eventually becomes smaller than the aerodynamic entrainment threshold. In this work, the objects we studied are all under the circumstances of saturated sand flux, which means that the aerodynamic entrainment can be ignored during the simulation. Although bed particles in saturated sand flux can not be directly lifted by air flow, their bombardment entrainment can still be slightly influenced by it. As shown before, the number and energy of ejected particles are controlled by the minimum transferred energy. This value varies while a particle on the surface receives the shear force from the wind. Figure 1 shows the free-body diagram of a particle exposed to the wind shear stress. If this particle can be barely lifted by the wind, the airborne shear force F s τ fw πd 2 /4, and the counter force F N from the adjacent support particle and the gravitational force m p g should be in the equilibrium state. In this state, τ fw is equal to the aerodynamic entrainment threshold stress τ ft . Being aware of these, we get the following equilibrium equation: where ϑ stands for the angle from these two particles' centerconnecting line to the horizontal surface. Otherwise, if the equilibrium state cannot be fulfilled, an effective mass m eff is then introduced into the model as follows: Substituting (17) into (18), we have the expression of the effective minimum transferred energy E eff as shown below: In practice, during simulation, E d2 in (14) and (15) is replaced by E eff . The aerodynamic threshold τ ft is deduced from Shao's work [32]. FIGURE 1 | Free-body diagram of a surface particle and its adjacent neighbor exposed to the wind shear stress.
Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 662389 A similar method has been used in Lämmel's previous work [23], leading to convictive results. By employing it, reptation will be enhanced. The bed surface becomes "looser," allowing more low-energy particles to move a tiny distance and leading to a higher frequency of midair collisions.

Hydrodynamic Governing Equations
In this study, the hydrodynamics is described using an RANS equation: where (·) denotes a temporal averaged quantity. For simplicity, in what follows, we note u u for the averaged fluid velocity. T represents the stress, which contains both viscous stress and Reynolds stress; F p is a forcing term describing the counter force exerted by the fluid-accelerated particles; ϕ p is the volume fraction of particles; and ∇p is the pressure gradient.
In the Eulerian field, ϕ p is obtained by averaging the information of Lagrange particles in grid cells. For a computational grid cell with particles in it, the following equation is applied: where N p is the total particle number in this cell, V c is the cell volume, and V p,i is the volume of the ith particle. The symbol 〈 · 〉 stands for ensemble averaging.
For the steady and homogeneous fluid field which we studied in this article, the inertia and horizontal stress gradients of the fluid are neglected. Regarding x as the streamwise direction and z as the vertical direction, (20) becomes the following: where F p,x (z) and ϕ p (z) are the particle counter forces in fluid direction and particle volume fraction, respectively, in a thin layer at a specific altitude z. ϕ p 1 means that the grid is full of particles. This will not happen because when the value of ϕ p is larger than the volume fraction of the bed particles, the hydrodynamic equations will not be calculated in the grid. Instead, the corresponding wind velocity will be set to 0 directly. A Prandtl's mixing length model with the kinematic turbulent viscosity ] t l 2 m zu/zz is used to close the RANS momentum equation [35]. τ f can then be expressed as follows: The mixing length scale l m is provided by the following: where κ 0.42 is the von Karman's constant and u * is the fractional velocity to describe the original wall shear stress of the grain-free fluid field (i.e., the fractional velocity far away from the sand bed).
Integrating (23) alone with the z direction, we obtain a drag partition equation as follows: τ p (z) is the grain-borne shear stress obtained by the following: By combining (24)-(27), the fluid field can be calculated in a straight manner in every time step.

Topography Geometry
The bedform surface in this model is triangulated and expressed by a series of key points. We propose to establish a dynamic topography model considering bombarding particles. Thus, the key point location is variable and decided by the vicinal real-time material exchange.
To simulate the terrain, the key points in this model can move vertically and are digitalized as their elevations. We separate the bed surface into several subregions in order to estimate these elevation values. As shown in Figure 2A, the subregion is usually defined as a rectangular area (dashed line box) around every key point and plays a similar role to a sand trap. By counting all particles that drop into it or escape from it within a computational time step, the elevation change of each subregion in this time step is deduced from the following equation: where h is the elevation of a subregion, N s is the net number of the particles trapped in this subregion, and ϕ b is the average volume fraction of particles inside the bed. Considering that all the sand particles are assumed to be spheres, ϕ b is set to 0.6 in this work. A is the area of one subregion. Then, we assign the elevation of each subregion to its central key point. What should be clarified here is that the terrain generated within every calculating step is not simply a ladder-shaped surface as it is in Anderson's simulation model [2]. Connecting all the key points of the surface generates a vivid microtopography with triangular slopes. As the sketch shown in Figure 2B indicates, the relative altitude of every particle to the surface h p is equal to the distance from the particle center to the slope surface right beneath it. The splash event occurs when h p is smaller than 0 m.
What is more, the splash function on a slope is still a matter of research and no ready-to-use model has been proposed. During the splash procedure in this model, the coordinate system of particle movement is steered toward the direction along the slope surface. In other words, for Eqs 8-15, all values of moving particles are converted to the values relative to the local bed slope. This kind of procedure seems arbitrary, but it is supported by Yin et al.'s recent work [39]. They found that for a constant impact angle relative to the bed, only the relative rebound angle shows a small decrease with the increasing bed slope angle. Other relative properties of rebound and ejected particles are nearly independent of the bed slope.
Another aspect one needs to pay attention to here is the angle of repose for every small bump and pit. We apply a natural method to regulate the behavior of particles interacting with the bed, that is, giving out the slope at every impactor's drop point to see whether it exceeds the repose angle. If so, and if the nearest key point is a salient point, the deposited particle will roll down the slope and settle at the lowest neighbor subregion. Similarly, for erosion situations, particles that impact at a pit with a sheer slope will entrain ejectors from higher neighbor subregions. The repose angle in this work is set to 30 + for original sand particles.

Particle Size Distribution
We can calculate polydisperse situations without changing the main frame of the particle movement description. The diameter of every particle in the air is recorded. Meanwhile, a simplified procedure is added to simulate the diameter distribution of particles on the bed surface. As we know, the whole computational domain is separated into server subdomains. Every subdomain can be considered as a sand trap which contains a specific number of particles describing the local surface elevation h. Now, we separate each subdomain into several "diameter bins." One diameter bin represents a specific particle diameter range. Particles with various diameters are sorted and dropped into corresponding bins automatically. By defining the sufficient number of diameter bins, one can describe the size distribution of bed particles. However, not all particles constituting the sand bed take part in the splash event. Only those near the surface influence the splash process. Therefore, here we introduce "affectable depth" into every subdomain standing for the depth from the surface, in which particles can be influenced by the impactor. Naturally, there is an "affectable layer" representing a layer with affectable depth covers on the sand surface. Each subdomain has its own affectable depth. Thus, for a specific subdomain, the development of surface particle size distribution has been converted to the development of the affectable depth and the diameter bins' volume ratio within this depth.
There are three kinds of situations while dealing with the changes in the size of diameter bins and the affectable depth. Here, we use the sketch in Figure 3 to show the dealing procedures for these situations. There are 12 units in Figures  3, 4, and units in a row describe the development process within one subdomain in a specific situation. Rectangular bars in different colors represent diameter bins of different particle diameter ranges. Area of bars corresponds to the bulk volume of particles within these bins. For simplicity, we set the number of the diameter bins to 3, and they stand for surface particle sizes around d 2,1 , d 2,2 , and d 2,3 , respectively. The affectable depth is represented by ψ, and the magnitude of ψ is distinguished by subscripts. ψ 0 is the initial value of ψ, which is also the lowest limitation of the affectable depth. Dashed boxes point out the affectable layer. Thus, in this sketch, the area ratio of rectangular bars within the dashed box corresponds to the surface particle size distribution of a single subdomain.
For the net erosion situation with ψ 0 as the original affectable depth before the collision event, since this affectable depth cannot be smaller anymore, particles deeper inside the sand bed with initial size distribution will replenish the newly spared empty space of the affectable layer when erosion happens. This process is   vividly portrayed in Figure 3A by the downward shifting of the dashed box and the following redistribution of bars inside the box. Therefore, we name these two steps as the "shifting step" and the "redistribution step." If the situation performs as net deposition, the shifting step becomes a "resizing step." In the resizing step, the arbitrary original affectable depth ψ 1 will increase to a larger value ψ 2 to make the space within the affectable layer sufficient to contain the extra deposited particles. This can be sketched as the resizing of the dashed box as it is shown in Figure 3B. After some time of development, the affectable depth at some locations will be larger than ψ 0 because of the previously mentioned net deposition situation. Net erosion happens at these locations, which leads to the third situation. The whole process of this situation is sketched in Figure 3C. It is a combination of those methods mentioned before. The resizing step will be carried out first. The affectable layer becomes thinner to meet the remaining surface particle quantity. The affectable depth stops decreasing when it reaches ψ 0 , and then the shifting step will take over.

Simulation Procedure
The following list gives the simulation procedure in a complete time step: (i) If this step is the first step, calculate the initial wind field using (24)-(26) with τ p 0. Then, release particles into the computational domain randomly. (ii) Evaluate (2) for every particle to obtain f fp and deduce the new location and velocity of particles from (1) with f cp 0. (iii) Search all particles and find out every collision pair. Renew the velocity of collision particles using (5). (iv) Calculate the effective minimum transferred energy E eff using (19). Find out particles below the surface and their impact location. Derive the properties of rebound and ejected particles from (8)- (15). Define particles which cannot jump higher than 1d as dead particles. (v) Obtain the elevation change of every subregion using (28). N s is derived from the number of the dead particles and the number of the ejected particles. (vi) Calculate the surface particle distribution according to the three steps described in section 2.6. (vii) Obtain F p from (22) and renew the wind field using (24)- (27). (viii) Go to step (ii) to start the next iteration.

Particle Transportation Flux on Flat Surface
We test the particle transportation performance of this model by comparing the simulation results of particle transport flux to the experiment results and the validated laws.
In this section, we use monodisperse particles. The computational domain is set as a 4000d × 10d × 1200d rectangular volume which has a fixed flat sand bed. During calculation, periodic boundary conditions are used in the x and y directions for both particles and the fluid field. As shown in Figure 2A, the mash of this model is curved along the bumpy surface and refined exponentially down to the height close to the bed, and then Δz becomes a constant near the bed surface because of the linear increase in wind velocity in the viscous sub-layer. The lowest grid is on the surface, and the smallest value of Δz just above the surface is controlled by the viscous length ]/u * , whose smallest value is in the order of 0.1d.
As we know, the steady and homogeneous particle transport conditions are characterized by three dimensionless numbers: the density ratio s ρ p /ρ f , the Galileo number Ga s g d 3 /], and the Shields number S ρ f u *2 /(ρ p gd), where g (1 − ρ f /ρ p )g is the buoyancy-reduced gravitational constant [26]. In this work, cases with various situations have been calculated to make the test more convincing. Among these multifarious cases, the value of s is set between 125 and 2098, Ga ranges from 9 to 304, and S varies from 0.01 to 1.
Sand particle transport can be induced by a bunch of randomly separated triggering particles. These particles bombard the sand surface, causing a chain reaction of bouncing and ejecting, which eventually leads to a saturated particle transport; the amount of sand particles trapped on the bed is equal to the amount of newly ejected ones. The saturated state is always satisfied during a short time period. In this work, the timescale is nondimensionalized by d/ g and represented by t * . The duration time of all cases is longer than t * 12000, which is sufficient for them to reach the saturated state. The dynamic threshold confirming particles which continue transport can be deduced then, and it is expressed as S d (threshold Shields number) or u * d (threshold friction velocity). One will be able to deduce the particle transport flux of a system after it becomes saturated. In this work, the mass flux Q of particles per unit width is calculated by using the following: where P stands for a collection of all moving particles above an area Δ. For aeolian sand transport, the relationship between Q and the Shields number S has been studied for many years. Q was first proposed to be the 3/2 power of S [6,20]. However, in this decade, wind tunnel experiments and numerical simulations indicate that this scaling law between Q and S is more likely to be linear for situations with smaller S [8,16,18]. To analyze the simulation results, here we introduce a nondimensional flux as follows: Figure 4 demonstrates the development of Q * versus S. We find that our simulation results roughly satisfy the scaling laws deduced from those previous works. Q * grows linearly on S for the cases with S < 0.1. When S increases to a value larger than 0.1, the scaling law imminently changes to Q * ∝ S 3/2 .
Recently, a unified expression of non-suspended sediment flux both in water and in air was proposed [26], which reveals the mesoscopic mechanism that controls the flux quantity and avoids the subsection description in a wide range of S. By this expression, the relationship between Q * and S can be written as follows: where c M and μ b are constant parameters. c M correlates with moving particles' fluctuation energy dissipation. μ b is a bed friction coefficient that characterizes the geometry of particle-bed rebounds. In this work, c M 2.3 and μ b 0.61. These results are close to those from Pähtz's work, in which c M 1.7 and μ b 0.63. From Figure 4, we notice that (31) fits very well with our results, which means that the simplifications of our model have not influenced the veracity of the aeolian sand transport results. What is more, when comparing the simulation results of this work with those of the 2D DEM simulation [12], one will find that our result is much closer to the data from wind tunnel observation. This demonstrates the superiority of the 3D simulation. Since for moving particles in 2D models, energy obtained from the wind will not be reallocated to the particle movement in the transverse direction, 2D simulation overestimates the horizontal sand flux in the stream direction. Another plausible reason for the different performance on flux prediction is the clearer definition of the sand bed in our model. In DEM simulation, every particle in the "static" grain bed is making tiny moves at all times. It is difficult to tell the moving particles from the static ones, especially near the solid-gas interface. Thus, the bed surface becomes a blurry concept. Sand flux in this kind of model contains a nonnegligible contribution from the particle creep inside the bed, which causes the overestimation.

Ripples Formed by Monodisperse Particles
To test the ability on reproducing the ripple development stages, in this study, we use our model simulating sand ripple emergence from a flat surface in different wind velocities. The computational domain is set to 4000d × 200d × 1200d, which has an 80d high deformable sand bed. Nine cases under typical wind-blowing sand conditions are calculated. The nondimensional properties of them are set to s 2098 and Ga 38, and S ranges from 0.008 to 0.07. The longest physical duration time of the simulations is set to t * 720000. All the other settings are the same as those used in the flat surface simulations mentioned before.
After computation, by analyzing the result of bedform topology, one can extract the bedform profile at every time step. Figure 5 shows the evolution of the bedform profile. Four subpictures correspond to the bedform at the wind velocity u * /u * d 1.7, u * /u * d 2.2, u * /u * d 2.8, and u * /u * d 3.4, respectively. From them, one can find that tiny periodic structures emerge rapidly from the flat sand surface. Some of them merge with others afterward. These tiny structures are the initial wavelength of aeolian sand ripples, and the merging process is called coarsening. Figure 6A shows the growth curve of the ripple amplitude. To calculate the average amplitude, one should first cut the computational domain into slices along streamwise grid lines, deducing several cross sections. Then, we perform autocorrelation C(Δλ, t) 〈h(x, t)h(x + Δλ, t)〉 − 〈h(x, t)〉 2 on the ripple profiles of every cross section, averaging the results of all slices to get C(Δλ, t). Finally, the average amplitude A can be deduced from A(t) 2 2C(0, t) . The linear stage is a period in which the amplitude A increases exponentially versus time. Small bumps rapidly emerge from the flat surface at this stage and immediately grow to regular ripples. After that, coarsening processes take over ripples' development and slow down ripples' growth on the amplitude. The x-coordinate value of the first peak on the C(Δλ, t) profile represents ripple wavelength λ at time t. As the results shown in Figure 6B indicate, there is a plateau in wavelength at the end of the linear stage for each wind velocity. According to Andreotti's  Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 662389 10 suggestion [5], this plateau indicates sand ripples' initial wavelength λ int . To avoid arbitrary treatment, in this work, we define λ int as the average wavelength in t * 100000 ∼ 150000. This time period contains the previously mentioned plateau of all simulation cases. We compare the deduced λ int with the experimental result in Figure 7. For clearer comparison, we assume that all these wavelengths obey the scaling law introduced by Durán [12]. Then, they can be nondimensionalized as λ * int λ int g/(u 2 d ρ f /ρ p ). We find that λ * int deduced from our model is very close to the experimental  result, which grows linearly with frictional velocity and almost vanishes at u * u * d . After the linear stage, the ripples' growth rate slows down. Checking the increase in wavelength in Figure 8, the simulation result of the wavelength exhibits a trending to a longer plateau, but generally, the growth rate of it holds as a constant. Wavelength data from field observations [38], from wind tunnel experiments [2,5], and from continued models [9,40] are consistent with the power law λ(t) ∝ t χ L , in which the growth rate 0.27 < χ L < 0.5. χ L in this work roughly fulfills this restriction and varies at the neighborhood around 0.5. Similar to experiments, χ L increases slightly with u * . However, this small increase is inconspicuous in this work. Amplitude grows coordinately with wavelength. Here, we check the ripple index (RI), which represents the ratio between wavelength and amplitude. The initial RI is a very large value with a quantity bigger than 100. During the linear stage, the RI decreases rapidly like it is shown in Figure 9 and eventually becomes a relatively constant value after regular ripples emerge. For aeolian ripples observed in the field, the RI usually ranges between 15 and 20 with a standard value of 18 for sand ripples and 15 for granule ripples [40]. This limitation roughly coincides with our RI results between t * 50000 and t * 300000, which are shown in the inset of Figure 9. For RIs deduced after t * 300000, they become smaller than the expected value. It is because of the unstopped growth of the ripple amplitude, which will be discussed in the next section.

Ripples Formed by Bidisperse Particles
Ripples' growth rate should eventually become unperceivable after a long time of development. The existence of aeolian sand ripples' saturated state has been proven by wind tunnel experiments [5]. In previous monodisperse cases, the simulation time is long enough to reach the saturated state. However, according to Figure 6, the "true" saturation has not been achieved at the end of the simulation. The incapability of reproducing the saturated bedform may be because of the lack of restriction on the increasing ripple scale. There are two plausible factors. The first one is that the homogeneous wind field used in this work cannot reproduce the complex distribution of τ f near the surface after the ripple grows to a nonnegligible size. This treatment loses the retroactions from the bedform. The second one is that the monodisperse particle is unable to reproduce the particle size distribution on the bed surface and in the air. Experiments, observations, and theoretical models show that differences in particle size distribution influence aeolian sand ripples' morphology [19,25,37].
The second imperfection can be conquered by using polydisperse particles. In this work, for simplicity, we use bidisperse particles to test the influence on ripple formation of particle size distribution. Similar to the computational settings in the last section, three cases with the same wind velocity are calculated under typical wind-blowing sand conditions. Case 1 is a monodisperse case with a particle diameter d which acts as a comparation. Case 2 is designed as a situation with an average diameter d. The diameter of those two kinds of particles used in case 2 are d + d/2 and d − d/2, respectively. Dynamic thresholds of particles in case 2 are all smaller than the given wind velocity. Case 3 is calculated with particles not smaller than d. The particle diameters in this case are d and 2.5d. It represents the situation in which the coarser particles' dynamic threshold cannot be reached by the wind field.
Ripples' amplitude development in all three cases is demonstrated in Figure 10, in which variables are nondimensionalized by the average particle diameter d a of each case. From it, one will surprisingly find that the cases with bidisperse particles can reach a relatively steady amplitude after a distinct rapid growth stage. This phenomenon proves the assumption mentioned before that the armoring from the coarser particles contributes to the stabilization of the ripple morphology. Whether the coarser particles could be lifted by the splash process or not, the effect of armoring can always be observed on ripples' development.
What is more, we can extract the average diameter distribution of surface particles as it is shown in Figure 11. D in it represents the local average diameter of particles on the bed surface. Particles in different sizes are regularly distributed in the computational domain for all bidisperse cases. Coarser particles prefer to gather on the crests of ripples. Meanwhile, fine particles tend to deposit on the trough. This coincides with the observation results from wind tunnel experiments and field observations.
The general view on the mesoscopic particle dynamic during ripples' development indicates that bigger particles are pushed to the crests because of the bombardment of saltating particles on the stoss slope. Our model has the ability to reproduce this process, which is very important for the simulation on some specific ripple forms, such as the so-called mega ripples.

Ripple Morphology in Transverse Direction
We can see from Figure 5 that larger ripples move slower, and smaller ones catch up and merge with them. Small bumps continually merge after their emergences. For a clearer observation, we run a wider monodisperse case, which is 4000d × 800d × 1200d. This case is only run at one specific wind velocity. Figure 12 shows the bedform topology from t * 20000 to t * 80000. The color represents the magnitude of DR, which stands for the deposition/erosion rate at that location, that is, DR Δh(x, y)/(dΔt). As we can see in these pictures, irregular structures gradually emerge from the flat surface and then strand to each other in the transverse direction, becoming relatively regular patterns. The initial wavelength appears very early at t * 40000 while the bedform is chaotic in the 3D view. Amplitude growth slows down until t * 80000 after recognizable ripples come out marking the end of the linear stage. However, ripples cannot become perfectly parallel lines after the linear stage in this wider computational domain. Defects such as disconnections and Y-shape branches arise just like features observed on the field or in the wind tunnel.
Deposition and erosion are obviously influenced by the bedform. In Figure 13, by comparing the deposition/erosion rate in the dashed circles with the corresponding surface topology, one can see that before the bedform becomes regular, deposition frequently happens at the gap between   Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 662389 14 staggered bumps, making them join together. At some locations, the deposition rate in the gaps is even bigger than those behind the adjacent bump. Then, as it is shown in Figure 12, on the ripple surface, deposition always takes place at the lee side, while erosions are all located at stoss slopes. From trough to crest, the erosion effect is reduced. It is because saltating particles bombard at the stoss slope, inducing reptating particles to climb up toward the crest. Near the top of the ripples, reptating particles jump over the ridge and deposit at the lee slope, making the deposition effect more remarkable close to the crest. This mechanism exists throughout the entire lifetime of sand ripples and motivates ripples to migrate along the wind direction.

Wind Profile Upon Ripple Surface
From the wind field profile corresponding to saturated aeolian sand transport, one can see whether the counter effect received from particles performs correctly. Meanwhile, differences of wind profile corresponding to various underlayer surfaces will be revealed. We will test the simulation performance on wind field calculation in this section. Although the fluid simulation in our model is simplified as a one-dimensional RANs problem, some average properties extracted from the simulation cases are still comparable with the experimental results.
The wind profiles on the flat surface and on the ripple surface are shown in Figure 14. Cases with the same s and Ga in different wind velocities are picked out from pervious sections. The wind profile on the ripple surface is deduced from the wind field within the time period when the initial ripple emerges. The sand flux weakens the wind velocity near the sand bed. An obvious focal point deduced from the wind profile of the ripple surface indicates a constant height of the saltation layer for different wind velocities, which is also obtained from wind tunnel experiments [24]. Furthermore, when extracting the aerodynamic roughness z 0 from the wind profile far from the bed, as it is shown in the inset of Figure 14, z 0 can be fitted linearly. It means that z 0 fulfills the equation of constant focal point theory as shown below: where z F and u F are the height and wind velocity at the focal point, respectively. We can see that z F deduced from our model is very close to that from Ho's experiment [15], while u F is closer to that in Durán's DEM simulation result [10]. Since (32) is derived after taking the existence of a unique focal point as a prerequisite, this result also points to the constant height theory on the particle transport layer. In fact, according to the simulation results, (32) holds all the way during the ripples' development process, which is shown in the inset of Figure 15.
Comparing z 0 deduced on the ripple surface to those on the flat surface, it exhibits a slight decrease, especially for results at smaller wind velocities. As we know, in the presence of particle transport, the aerodynamic roughness z 0 deduced from the wind profile far from the bed has no relation to the topography of geometrical roughness and the size of the viscous sublayer. The main contributing factor to z 0 is particles' trajectories. Thus, the counterintuitive decrease in z 0 for small wind velocities is not directly related to the surface topology but comes from the change of particles' movement manner.
To study the statistical properties of bouncing particles' trajectories, z F and u F are the values one should pay attention to. z F describes the bouncing height of particles, while u F reflects length as a function of the bedform topography and particle properties, which is very important for the research of mega ripples and dunes. Compared to the complete DEM model, our model has a clear-cut definition of particle status. There are no blur definitions in this model between the moving particles and the static bed. Meanwhile, the simplicity of it allows us to add more complex mechanisms into the model to study more influential aspects of the sediment transport, such as complex turbulence structures or the cohesion. This model reproduces the morphology development of aeolian sand ripples qualitatively and quantitatively. The relationship between ripple wavelength and wind velocity fits well with the wind tunnel experiment data. What is more, due to the simplified treatment on particle size distribution calculation, one can deduce the final stable wavelength of aeolian sand ripples, which has not been deduced using any direct simulation models before.
The inverse grading in aeolian ripples has been reproduced as well. Although Anderson and Bunas studied this phenomenon using their cellular automaton model [3], some new results can be deduced from our novel model in the future. This is because compared to the automaton model, our model has realistic saltation trajectories, while the cellular automaton model simplifies the saltating particles as a series of randomly distributed impact beans. Durán et al. found that the trajectories of saltating particles can be modulated by the wavy surface [12]. These resonant saltation trajectories contain the information of the ripples' geometric scale. On the other hand, our model can calculate polydisperse problems with arbitrary particle size distribution. This advantage gives us an easier way to perform careful research on the relationship between the inverse grading and the particle size distribution of the original bed.
This model provides a 3D view on observing ripple formation. Small bumps merge with each other and migrate along the wind direction as soon as they emerge from the flat bed. It is caused by the inhomogeneous distribution of the erosion/deposition region. Deposition always happens at the lee side of ripples, while erosion takes place at the stoss side.
Lastly, by testing the wind profile on the flat surface and on the ripple surface, we find that this model can simulate the average wind field on the ripple surface as well. It reproduces the focal point of different wind profiles and shows the development of particle transportation during bed form deformation. The dynamic threshold of bouncing particles increases to a certain value and holds at it at the early stage of ripple formation. Then, it increases linearly with time after the initial ripple is formed.
This model performs well on simulating aeolian sand ripples' development on Earth. Furthermore, since all equations contained in this model are derived from basic physical laws, it can be directly used in the research studies on planetary geomorphology as well. However, due to the simplification on the wind field and the surface particle size distribution, this model cannot deal with big structures such as meter-scale mega ripples and kilometer-scale sand dunes because we believe that the bigscale turbulence structure may affect the particle trajectory profoundly. Expecting this model to be the most feasible method to reveal the immanent properties of mega ripples in planetary geography research, we will continually improve it in future works.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.