A Numerical Research on the Relationship between Aeolian Sand Ripples and the Sand Flux

: Theoretically, the sand ﬂux will not change after the wind-driven sand particle transport reaches the saturated state. However, it has been found in many wind-tunnel experiments that the sand ﬂux will gradually decrease with time in long-term particle transport duration and will eventually reach a new stable state. In this work, we used numerical simulations to study the source of this kind of decrease and found it is caused by the sand ripple on the bed surface. The ripple index showed a strong correlation to the sand ﬂux, and it decreased during the initial stage of the ripple formation. With a simpliﬁed theoretical model, we found the linear relationship between the Shields number and the particle transport load holds. However, the slope of this relationship and the dynamic threshold of particle entrainment decreased with the ripple index. As the sand ﬂux scales linearly with the particle transport load, we ﬁnally derived an expression that describes how the sand ﬂux on the ripple bedform varies with the wind strength. From this expression, we found the sand ﬂux increases with ripple index, and it was easier to be inﬂuenced by the ripple bed form in small wind strength.


Introduction
While subject to a continuous wind, small particles in the granular bed can always be transported by the means of suspension, saltation, or reptation. As illustrated in Figure 1, comparing to the suspension particle who drifts a very long distance without touching the ground (usually in the scaling of kilometers), saltation and reptation particles interact with a local bed surface more frequently, i.e., they bounce forward along the wind direction. Saltation particles and reptation particles dominate the magnitude of the total particle flux in the wind-blow-particle situation [1]. The only difference between them is that the former inherits more energy from the wind, making it rebound after every impact and moving forward continuously. Because of these reasons, we realize that the local topology of the bed surface can seriously influence the aeolian particle transport. On the other hand, the bed surface consists of grains that can be modified by the particle transportation as well. Among all the aeolian-caused bed forms, the sand ripple is the most observed one. Its morphological characteristics are stripes perpendicular to the wind direction and almost symmetrical in the transverse direction. The ratio between the ripple's wavelength and the amplitude is called the "ripple index." For ripples on earth, this value ranges from 15 to 25, with a standard value of 18 for sand ripples and 15 for granule ripples [2].
The influences on sand flux from the aeolian sand ripple surface have been seldom studied systematically. The incipient ideas come from the observation on wind tunnel experiments in which the initial sand bed is always flat, since the ripple bedform emerges and grows quickly after particles' transport starts. It is obvious to speculate that this kind of surface deformation will cause changes on the trajectory of hopping particles, which in turn would adjust the transport flux. By assuming this, Rasmussen and Mikkelsen [3] observed the sand flux changing upon a pre-rippled sand bedform in a wind tunnel and found the sand flux was constant during the first 30 min; then, it decreased after 40 min. After 75 min, the flux had dropped to about 75% of the initial rate. However, this result was deduced from a uniform time sampling with vertically stacked traps, which has poor accuracy on the temporal resolution and mass flux quantity. What is more, as this work used polydisperse dune sand as the base material, influences from the surface particle size distribution were added into the results. In order to reveal more on the long-term stability of mass flux upon erodible sand surface, Rasmussen et al. then tested the flux changing of the uniform sand by using a particle image velocimetry (PIV) technique [4]. Starting with a flat surface, sand transport first increased to reach the saturated state. Then, sand flux continuously decreased until the ripple surface developed to a fully developed state. The total decrease quantity of mass flux was between 75% and 90% for all the situations they covered. They believe this decrease comes from the slight increase in ripples' from drag, but this theory has not been verified yet. Wang et al. in their wind tunnel experiments also found an obvious sand flux variation during the development of sand particle transport [5]. However, they did not relate it to the ripple formation at the bed surface. Tong and Huang directly simulated the sand particle transport over a ripple surface by coupling every particle's movement to the large-eddy simulation wind field [6]. They brought a very different result that the ripple bedform in their work enhances the value of sand flux. This distinct result may be because of the unreal ripple amplitude, which is as high as 2 cm. The bumps in their work are about ten times larger than that of the normal aeolian sand ripples on earth. No matter for wind tunnel experiments or for numerical simulations, these bumps perform more like sand dunes than sand ripples.
A sketch of the main particle transport forms in the aeolian particle transport system.
Then, we turned to look at the studies on sand ripples' development. It has been intensively studied for many years. The seminal work was established by Bagnold who predicted that the emergence and development of sand ripples are caused by nothing but the interaction between bouncing particles and the erodible bedform [1]. Since then, many experimental works [7][8][9][10], theoretical analyses [2,[11][12][13][14][15][16], and numerical simulations [17][18][19] were done to study the development of sand ripples. However, within all these works, no one has studied the relationship between the ripple morphology and the sand particle transport.
On these bases, we find that the study on how the ripple bedform influences the sand flux is very important. The sand flux observed on a specific ripple bed may perform differently to that observed on a flat surface. Moreover, the sand ripple has different development states, which will cause different impacts on the sand flux. Because of the lack of the attention on the surface bedform while studying the sand flux, the observation result may become very different if we compare one experiment to another. Thus, we need to obtain an expression to describe the relation between the ripple topology and particle transport statistics. It can help us to prevent the ripple-caused inaccuracy.
In this article, we focus on the numerical simulation of particles' transport upon sand ripples in aeolian systems, which will help us learn more about the importance of the tiny sand ripple structures. Moreover, this study theoretically revealed the source of the observed decreasing sand flux during wind-tunnel experiments.

Methodology
In brief, the numerical model used in this work is a simplified discrete element method (DEM) model. The full-scale DEM model simulates the movement of all the particles in the system through a series of calculations considering particle deformability and using the explicit time integration algorithm to resolve particle collisions. These features lead to excessive computational requirements, which is unacceptable for the long-time mesoscale numerical simulations of this work. To avoid this, our model just tracks moving particles in the air and gains the post-collision velocity components of particles by solving the momentum equations in connection with Coulomb's law of friction.
Sand particles were assumed as small spheres. Ignoring the influences of particles' rotation, the equation of particle velocity components is simplified as where m p is the particle mass; v is the particle velocity; g is the gravity acceleration; and f f p and f cp stand for the external forces of hydrodynamical and particle collisions, respectively. In aeolian situations, the density ratio between particles and the fluid, i.e., s = ρ p /ρ f , is very large. Thus, the hydrodynamical force f f p here is simply dominated by the drag force from the fluid: where d is the particle diameter; u is the wind velocity located at the particle location; and C d is the drag coefficient. We used the following approximation that where Re p = |u − v|d/ν is the particle Reynolds number; ν is the kinematic viscosity coefficient of the fluid; C ∞ d ∼ = 0.5 is the drag coefficient of the grain in the turbulent limit (Re p → ∞); and Re c p ∼ = 24 is the transitional particle Reynolds number [20].
Collision events take place while the centroid distance between a pair of particles is smaller than the sum of their radii. Comparing with the iteration time step, the collision happens during a very short time. Thus, we did not calculate f cp directly but deduced the post-collision velocity v 1 of a particle by [21] where the variables with subscript 1 or 2 refer to different particles, 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 effective restitution coefficients: where ε and µ are microscopic restitution coefficients for the normal and tangential components, respectively; η = m 1 /m 2 is the mass ratio of two colliding particles. In this work, ε = 0.9 and µ = 0 while calculating particle-bed collisions. For mid-air collisions, the tangential restitution coefficient µ was calculated from [22] where C f = 0.4 [23] is the coefficient of friction, and V n and V t are relative velocities between two colliding particles in normal and tangential directions, respectively. According to Equation (3), Lämmel et al. deduced a splash function after considering all possible impacting situations [24]. If a particle impacts the bed, its movement after collision can be quantified in terms of the mean restitution coefficientē and a distribution of the rebound angle θ 1 . For a given impact angle θ 1 , where Some bed particles are entrained by the residual energy of the impactor, and their kinetic energy E 2 can be drawn from a log-normal distribution where is the minimum energy required for an ejecta considering aerodynamic entrainments [25]. τ f w is the wind shear stress derived from simulation, and τ f t is the aerodynamic threshold deduced from reference [26]. The number of entrained bed grains follows the energy balance approach underlying Equation (10): The ejection angles of all splashed low energy particles were set to be 90 • , and their initial positions were randomly distributed near the impact point.
For the steady and homogeneous fluid field that we considered in this paper, inertia and the horizontal stress gradients of the fluid were neglected. Regarding x as the streamwise direction and z as the vertical direction, the vertical stress gradient of the fluid can be calculated by where τ f is the airborne shear stress; φ p is the volume fraction of particles; is the counter force from N p particles in the cell volume V c ; and . stands for ensemble averaging.
A Prandtl's mixing length model with the kinematic turbulent viscosity ν t = l 2 m |∂u/∂z| was used in this work. τ f then can be expressed as The mixing length scale l m was provided by where κ = 0.42 is the von Karman's constant; u * is the frictional velocity describing the original wall shear stress of the particle-free fluid field. The bedform surface in this model was triangulated and expressed by a series of key points. Elevation variation ∆z at a key point was calculated at every iteration step from the following equation: where N s is the net number of the particles trapped/escaped from the bed surface near this point; φ b = 0.6 is the average volume fraction of the grain bed. ∆ g is the area of one surface grid. The terrain generated in this way was not simply a ladder-shaped surface. Connecting all the key points of the surface generated a realistic micro topography with slopes. During the splash procedure, the coordinate system of the particle movement was rotated towards the direction of the slope surface. Additionally, the angle of repose of the sand particles was taken into consideration. The angle of the slope was not allowed to be larger than 30 • .
Our numerical model implements all the algorithms and equations mentioned above. A more detailed description of the methodology can be found in reference [19].

The Sand Flux Variation During Ripple Formation
Using the numerical method described above, we first simulated the development of the ripple surface and checked the consequent sand flux variation. It is worth mentioning that in this work we used monodisperse particles to avoid the influence of the particle size distribution. The size of the computational domain was x × y × z = 4000d × 40d × 1200d, and the boundary condition was periodic in the x and the y direction. Witĥ g = g 1 − ρ f /ρ p , the Shields number S = ρ f u * 2 / ρ pĝ d in all cases ranged from 0.017 to 0.095. Here, we kept S below 0.1, because according to references [19,27,28], for larger Shields number S, the relationship between sand flux Q and S will no longer be linear. The range of S for which the linear scaling law Q ∝ S becomes invalid is called the Bagnold-like region. In this region, the average velocity of the particles in the air is not only related to the particle's dynamic thresholds (threshold Shields number S d or threshold friction velocity u * d -they correspond to the smallest wind strength that can maintain the particle transport) but also controlled by the velocity of particles flying above the transport layer [29]. Except for the Shields number, particle transport can be influenced by the other two dimensionless numbers, which are the density ratio s = ρ p /ρ f and the Galileo number Ga = sĝd 3 /ν. For the subsequent simulations, we fixed these two parameters as s = 2098 and Ga = 38.
At the beginning of our simulations, the sand particle transport was triggered by randomly distributed inducing particles. These particles arouse other resting grains from the surface causing chain reactions and eventually leading to a saturated sand flux. Before calculating ripple development, we reserved enough time (t * = t/ d/ĝ from 0 to 12,000) for the system to saturate. Equation (15) was unused during this time, and the sand flux rapidly reached to a steady value. We used this procedure to make sure that the following sand flux variation was only caused by the surface deformation. Then, the surface started to deform. As the space-time diagrams show in Figure 2, small bumps emerged from the sand bed regularly and migrated slowly in the wind direction. It is obvious that the wavelength (λ, the horizontal distance from crest to crest) and the amplitude (A, the vertical distance from trough to crest) of ripples grew with time. To have a clearer look at the morphological development of sand ripples, we calculated the autocorrelation C a (∆λ, t) = z(x, t)z(x + ∆λ, t) − z(x, t) 2 of the ripple profiles and averaged the results over the y direction. Then, the average amplitude was calculated from A(t) = 2 2C a (0, t), and the average wavelength λ(t) corresponds to the x-coordinate of the first peak on C a (∆λ, t). As the results show in Figure 3, the ripple wavelength and the amplitude both increased with time. Moreover, the growth rate became more and more insignificant, and it eventually became zero at the fully developed state of the aeolian sand ripples. As the bright yellow background demonstrates, in Figure 3a, there was a stage where the ripple amplitude grows exponentially. It is called the initial stage. Former works [19,30] have proved that before this stage (t * < 35,000 in this work), the bed surface is noticeably three-dimensional. Small bumps emerge from the surface at an arbitrary location. These bumps grow in the transverse direction and link to each other. Eventually, they become the primitive form of aeolian sand ripples. Then, because of the symmetrical ripple form, the bed surface no longer varies in the y direction. More and more low-energy particles tend to climb along the stoss slope of the sand ripple causing an exponential growth on the amplitude. From theoretical analyses, we have already known that this exponential-increase amplitude comes from the first-order contribution of the wavy bed surface [11,12,15]. The most unstable mode of the surface profile caused by this contribution represents the initial wavelength, which corresponds to the first plateau of the wavelength in Figure 3b. Then, after the initial stage (t * > 70,000 in this work), coarsening takes place. It makes adjacent bumps merge to each other. We can see there was a continuous growth of the ripple wavelength in the initial stage and a subsequent step growth in the coarsening stage. The later abrupt wavelength growth was caused by merging ripples.
Next, we checked the sand flux variation during ripple formation. In this work, the horizontal sand flux Q was calculated as the average horizontal particle momentum per unit bed area, i.e., Q = ∑ N m p v x T /∆, where N is the particle number in the air; ∆ is the bottom area of the entire computation domain; and . T represents the time average. We nondimensionalized Q via Q * = Q/ ρ p d sĝd . The temporal evolution of Q * in different wind strengths is shown in Figure 4a. According to former experimental and numerical works [19,[31][32][33], we have already known that there is a linear relationship between Q and the Shields number S. Thus, instead of Q * , here we checked Q * /(S − S d ) to reveal more information. It is worth mentioning that the threshold S d used here was evaluated at t * = 20,000. As the result shows, Q * holds a relatively stable value in the beginning of bed deformation. Then, an obvious decrease happens within the initial stage of the ripple formation. After the initial stage, a new stable value of Q * was reached. Except the decrease in sand flux, we noticed Q * /(S − S d ) derived from different wind strengths no longer coincided to each other after the initial stage. This can be caused by two possibilities: (a) S d varies during the initial stage or (b) the relationship between Q * and S − S d is not linear on the ripple surface. The latter is less likely to be true, because as we know, the linear relationship has already been proved by many experiments, which are always performed on ripple surface.  We also checked the temporal development of the ripple index (RI = λ/A) in Figure 4b and found RI increases until the initial stage begins. One needs to be aware that the value of RI does not represent the ratio between λ and A before the initial stage. Because as mentioned before, the typical ripple form has not emerged at this time. Thus, the increase in RI in this stage represents nothing, and the flux was unrelated to this variation until the initial stage. The sand ripple started to grow in the initial stage. RI decreased in this stage, which is similar to Q. This decrease lasted to the end of the initial stage, when RI and Q all reached to a relatively stable value. However, according to Figure 3, λ and A held their increase after the initial stage. Thus, combining all the results, we can draw a conclusion that Q is controlled by the ripple index RI. Additionally, we noticed that the final RI derived from our model was about 20. This coincides with the results from wind tunnel experiments and field observations [1,7,9,10]. However, they have not revealed the decrease in RI within the initial stage.

The Relationship between the Sand Particle Transport and the Ripple Index
To perform further studies on the relationship between Q and RI, we ran several extra cases with pre-rippled sinusoidal wavy surfaces (triangular wavy surfaces and sinusoidal wavy surfaces showed no differences on the result). For all the cases, the amplitude A was set as 8d, which is observed in Figure 3a as a typical amplitude within the initial stage. The wavelength varies from 50d to 1000d, making the ripple index RI range from 6.25 to 125. The study on the situations with RI < 15 was necessary, because for some specific conditions, the sand bed will develop to mega ripples, which has a much smaller ripple index [34]. For every pre-rippled simulation case, the code module that contains Equation (15) is disabled, and we reserved enough time for all the cases to ensure them to reach the saturated state. This procedure helped us to strictly control the value of RI. To derive authentic results, for all the pre-rippled cases, the density ratio s was set between 1325 and 5300; Ga ranged from 27 to 60; and S varied from 0.01 to 0.1. The other settings were the same as the settings of the ripple formation cases mentioned before.
We first tried to simplify the problem. As Figure 5a shows, reptating particles and saltating particles bounced forward along the wavy bed surface. There are two typical movement modes for them. One is hopping from a stoss slope to another stoss slope or from a lee slope to another lee slope (black lines in Figure 5a demonstrate this kind of particle trajectories); the other one is hopping from a stoss slope to a lee slope or from a lee slope to a stoss slope (red lines in Figure 5a demonstrate this kind of particle trajectories). We defined the former trajectories as monotone trajectories and defined the latter as nonmonotone trajectories. As Figure 5b,c shows, if all the particles have monotone trajectories and if we only study the ensemble average results, the whole system can be estimated by the combination of two simpler sub-systems, which only have a stoss slope or a lee slope. The particle with a non-monotone trajectory can be considered as a link, which transmits energy between the system shown in Figure 5b,c. We checked this energy-transmitting capability by counting the particle number that impacts the surface or take-off from the surface. N in stoss represents the total number of impact particles that hit the stoss slope, and N out stoss represents the total number of the particles take-off from the stoss slope. Similarly, we also had N in lee and N out lee for the lee slope. We checked N lee /N stoss for the impact and take-off particles, because this ratio was useful for the subsequent studies. In Figure 6, we found N in lee /N in stoss decreased linearly with 1/RI, meanwhile, N out lee /N out stoss decreased quadratically. These decreased relationships are all independent to the wind strength and the particle/fluid properties. Moreover, there was a difference between these two ratios, which is demonstrated by the bright yellow shade. This difference was due to the contribution of the non-monotone particle trajectories. Because if all particles have monotone trajectories, we can easily derive N out stoss = N in stoss and N out lee = N in lee . This yields an expression N in lee /N in where we defined this ratio as r a . To continue our study, we adopted this simplification. For simplicity, we assumed r a can be described by the linear fit of all the results in Figure 6, which is demonstrated as the dashed line. Thus, where Γ a = 5 is a constant that will not change with the wind strength or the particle/fluid properties. is the ratio between the impact (take-off) particle number at the lee slope and that at the stoss slope. The black line is a linear fit of N in lee /N in stoss . The red line is a quadratic polynomial fit of N out lee /N out stoss . The dashed line is a linear fit for all the data. The bright yellow shade highlights the difference between the impact particle results and the take-off particle results, reflecting the influence from non-monotone trajectories.
Similar to the definition of Q, we defined the particle transport load M = ∑ N m p T /∆ to represent the average total mass of particles transported above the bed surface per unit bed area. As Figure 5 shows, we separated the whole system to two sub-systems. It is worth mentioning that these two sub-systems are not totally independent to each other. They share the same wind field. For a better understanding, one can imagine that we used the result derived from Figure 5d to estimate the result in Figure 5a. The former only has a few neglectable non-monotone trajectories in the middle. In the following work, for the particles interacting with the stoss slope, we used subscript s to represent their statistical quantities. For the particles interacting with the lee slope, we used subscript l. Then, we could write the stress balance on the bed surface [35] − P zx = τ − τ f w − M s g sin α s + M l g sin α l , where P zx and P zz are horizontal and vertical grain-born stress applied on the bed surface, respectively. They arose from inter-particle or particle-bed contacts. α s and α l are slope angles. We assumed α s = α l = α. τ = ρu * 2 is the air-borne shear stress far away from the bed surface; τ f w is the air-borne shear stress on the bed surface. According to the drag partitioning theory [36], for a saturatde wind-driven particle transport system, τ f w = τ d , where τ d = ρu * 2 d is the shear stress corresponding to the dynamic threshold. Moreover, recall the definition of r a in Equation (16); we find it can also represent the ratio of M, i.e., M l /M s = r a . Use M = M l + M s , and define a parameter µ 0 b = −P zx /P zz . Then we can derive where M * = M/ ρ p d is the dimensionless form of M; µ b = µ 0 b cos α− s sin α(r a − 1)/ (s − 1)/(r a + 1). It shows a linear relationship between M * and S. From the simulation results shown in Figure 7, we proved that M * ∝ S always holds for different RI, different s, and different Ga. Additionally, we found µ b is obviously influenced by RI. The weak effect from s can be observed from the result, as well. To obtain the relationship between µ b and RI, we used α = cot −1 (RI/2) to roughly estimate the slope angle. Then, we have which is derived from Equations (17) and (20).  Figure 7. The relationship between the particle transport load M * and Shields number S. For filled symbols, different colors correspond to different RI. For hollowed symbols in the insets, different colors correspond to different particle/fluid properties. Two insets show the results of two specific RI. All the solid lines represent linear fits.
As Figure 7 shows, not only µ b varies with RI; S d also has a relationship with RI. According to the S d expressions on the slope [28], we similarly scaled S d as Noticing that for the flat surface situation, RI becomes an infinitely great value, which leads to µ b = µ 0 b and S d = S 0 d . Thus, according to Equations (20) and (21) with 1/RI = 0, we can directly derive S 0 d and µ 0 b from our flat surface results. It is worth mentioning that the "flat surface" here refers to the bed surface before the initial stage in the ripple formation simulation cases.
For the cases with s = 2098 and Ga = 38, we display the relationship between µ b and RI in Figure 8a. µ b decreased with RI, and the decrease rate dropped rapidly. µ b became more and more close to µ 0 b as RI increased. The relationship between S d and RI had the similar rule (Figure 8b). What was different was that the S d first dropped to a value smaller than S 0 d at a large RI. Then, as RI continued growing, it gradually increased to the value close to S 0 d . With derived µ 0 b and S 0 d , the results of Equations (21) and (22) are also shown in Figure 8 as the solid lines. Both of them roughly fit the simulation results, but Equation (22) cannot reflect the overly large decrease in S d .  (21) and (22), respectively. All the results are corresponding to the cases with s = 2098 and Ga = 38. RI = 20 is the typical ripple index of the fully developed aeolian sand ripples on earth [2,9,10]. According to our results, what varies most from the flat surface to RI = 20 is not the dynamic threshold S d but the parameter µ b . So, we can say µ b is a parameter that reflects the bed surface geometry. Substituting Equations (21) and (22) into Equation (20), we obtained a brief expression of M * where, S r = S/ 1 − 2sΓ a / µ 0 b RI(s − 1)(Γ a − 2RI) .
Finally, we checked the relationship between M * and Q * . As Figure 9 shows, with specific s and Ga, for all the Shields numbers studied in this work (S < 0.1), their relationship was close to linear. It is interesting that the ripple bed surface barely influenced the result. No obvious variations were observed between the results that corresponded to different RI. Thus, the relationship Q * ∝ S − S d was still valid for the sand particle transport upon the ripple surface. What was different from the result on the flat surface was the dynamic threshold S d , and the proportionality coefficient between Q * and S − S d were controlled by the bed surface. According to Equation (23), we can also scale Q * as Using this expression, we can estimate the relative quantity of the sand flux decrease during aeolian sand ripple formation in wind tunnel experiments. On the flat surface, the dynamic threshold u * d can be estimated as u * d = 0.8u * t [1,37]. u * t is the static threshold, and it can be calculated by u * t = 0.0123 sgd + 3 × 10 −4 / ρ f d [26]. Then, we can derive S 0 d . We used Q 0 to represent the sand flux on flat surface. We easily derived We chose five common RI, which varied from 10 to 30. These RI represent different aeolian sand ripples that emerged in wind tunnel experiments. Figure 10 exhibits the relationship between Q/Q 0 and S for a specific particle diameter d = 250 µm. One can find that for a specific ripple index, the sand flux caused by small wind strength is easier to be influenced by the ripple surface. Q/Q 0 rapidly approached to a constant value that was smaller than 1 when the wind strength increased. To have a clearer look, we extracted the Q/Q 0 values that correspond to S = 0.01 and S = 0.1. As shown in the inset of Figure 10, Q/Q 0 at these two wind strengths all increased with RI. For all the RI, the ripple bed was more influential to the sand flux in the situations with smaller wind strength. This coincides to the experiment result shown in the inset of Figure 4a.

Conclusions
Wind tunnel experiments reported a sand flux decrease during the long-term sand particle transport developments. In this work, we related this decrease to the sand ripple formation at the bed surface. From the numerical simulation and the theoretical analysis, we derived four main conclusions in this work:

1.
We simulated the emergence and development of aeolian sand ripples and found an obvious decrease in the ripple index RI at the initial stage of ripple formation. This RI variation coincides to a sand flux drop, indicating a strong connection between them.

2.
Particle transport load M decreased with RI. Its linear relationship with S always holds; only the proportionality coefficient µ b of this relationship and the dynamic threshold S d varied with RI. We derived their expressions in Equations (21) and (22). 3.
The relationship between M and the sand flux Q was barely influenced by RI. For the S smaller than 0.1, we roughly had Q ∝ M.

4.
We presented the relationship between Q, S, and RI in Equation (24). From this expression, we found the sand flux Q derived at the ripple surface was always smaller than that on the flat surface. A larger RI corresponds to a smaller Q. Cases with small wind strength are easier to be influenced by the ripple bed.
With these conclusions, we finally explained why sand flux decreased during wind tunnel experiments. We hope this work could show some helpful ideas to future experiments. We learned that for some situations, while studying the sand transport, the sand ripple on the bed surface is no longer a neglectable tiny structure. Moreover, Equation (24) can be used for the field observations. By measuring the ripple index and the sand flux at two different wind strengths, one can derive the flat surface dynamic threshold S 0 d , which is a rough estimation of the local sand particle dynamic threshold. We need to do wind tunnel experiments in the future to verify this usage and to conduct a further improvement in the expression of Q.