A Numerical Study of Aeolian Sand Particle Flow Incorporating Granular Pseudoﬂuid Optimization and Large Eddy Simulation

: A numerical investigation of aeolian sand particle ﬂow in atmospheric boundary layer is performed with a Eulerian–Eulerian granular pseudoﬂuid model. In this model, the air turbulence is modelled with a large eddy simulation, and a kinetic–frictional constitutive model incorporating frictional stress and the kinetic theory of granular ﬂow is applied to describe the interparticle movement. The simulated proﬁles of streamwise sand velocity and sand mass ﬂux agree well with the reported experiments. The quantitative discrepancy between them occurs near the sand bed surface, which is due to the di ﬀ erence in sand sample, but also highlights the potential of the present model in addressing near-surface mass transport. The simulated proﬁles of turbulent root mean square (RMS) particle velocity suggest that the interparticle collision mainly account for the ﬂuctuation of sand particle movement.


Introduction
Aeolian sand transport is a leading factor in the geophysical evolution of arid and semiarid regions on earth [1][2][3]. It also causes environmental problems such as air quality deterioration and land degradation [4,5]; the sand transport in atmospheric boundary layer is critical to the latter phenomenon. The interpretation of such gas-particle two-phase flow relies on the understanding of the momentum exchange between the two phases, which is usually specified as the distributions of velocities and mass flux of the sands [6,7].
Intrusive equipment such as ruggedized probes has been used to quantify these distributions at satisfactory temporal resolutions [8,9], which nevertheless inevitably disturb the measured flow fields. The use of nonintrusive measurements, such as particle image velocimetry and high-speed photography [7,[10][11][12], avoids this problem. However, they have difficulties in addressing "dense" flow under strong wind or near the sand bed surface.
Numerical simulation has been an indispensable method for accompanying the measurements [2,13]. There are generally two methods for simulating gas-particle flow. One method tracks the information from individual particles in the Lagrangian point of view [14,15]. Such a discrete model constructs the interparticle and gas-particle interactions by virtue of the Newtonian laws; therefore, it would become inefficient in cases of an extremely large numbers of particles. The other method treats moving ∂ ∂t α g ρ g + ∂ ∂x j α g ρ g u g,j = 0 ∂ ∂t α g ρ g u g,i + ∂ ∂x j α g ρ g u g,i u g,j = −α g ∂p g ∂x i + ∂ ∂x j α g µ gl + µ gt ∂u g,i ∂x j + ∂u g, j ∂x i − 2 3 ∂u g,i ∂x j δ ij + α g ρ g g i + β u s,i − u g,i (1) where ρ g = 1.225 kg/m 3 is the gas density; α g , u g,j and p g are the volume fraction, the velocity component in the jth direction and the pressure of the gas phase, respectively; δ ij is the Kronecker delta, g i is the component of gravity in the ith direction, u s,i is the velocity component of the sand phase in the i th direction, β is the drag coefficient between the two phases and µ gl = 1.8 × 10 −5 kg/(m·s) is the gas viscosity. µ gt = ρ g (C g ∆) 2 |S g | is the gas subgrid-scale eddy viscosity with |S g | = (2 S g,ij S g,ij ) 1/2 , where C g = 1 is the Smagorinsky model constant, ∆ = (∆x∆y∆z) 1/3 is the filter width and S g,ij = (∂u g,j /∂x i + ∂u g,i /∂x j )/2 is the filtered rate of strain of the gas phase [23,24]. It is undeniable that the Smagorinsky mode does not well model the process of energy backscatter. However, such kind of Smagorinsky LES turbulence model has still successfully and reliably been used in the simulations of aeolian sand transport in equilibrium state [22]. Therefore, considering the emphasis on sand phase analysis and the difficulty of model implementation, the classical Smagorinsky model is adopted.
The filtered equations of motion for the "continuous" sand phase are given by Mathiesen et al. [17] and shown in Equation (2): where a smooth function is written as ϕ 2 (α s ) = arctan[96(α s − α s,min )]/π + 0.5, d s is the average diameter of the sand particles, Ω = 31 • is the angle of internal friction of the sand particles and I 2D is the second invariant of the strain rate tensor determined by the rate of the strain tensor of the sand phase flow S s,ij = (∂u s,i /∂x j + ∂u s,j /∂x i )/2. Under the Eulerian-Eulerian framework, an eddy viscosity model of the sand phase µ st is introduced to take into account the subgrid-scale motion of the sand particles [34]. µ st is modelled using an eddy viscosity assumption: µ st = α s ρ s (C s ∆) 2 (2 S s,ij S s,ij ) 1/2 , where C s = 0.14 [23]. The bulk viscosity is described as ξ s = 4α s 2 ρ s d s (1 + e)g 0 (θ/π) 1/2 /3. Then, the transport equation for turbulent kinetic energy is shown in Equation (6): Specifically, k s is the conductivity of granular pseudotemperature, as shown in Equation (7): γ s is the dissipation due to inelastic collisions, as shown in Equation (8): And Φ = −3βθ is the fluctuating energy exchange between the two phases, for which the details are shown in Equation (9): With Re s = α g ρ g d s | u g − u s |/µ gl as the particle Reynolds number [35,36]. Discretization of the governing equations is performed on a collocate grid using the finite control volume technique. The pressure-velocity coupling is achieved by the SIMPLE method with momentum interpolation [37]. The QUICK scheme is adopted for the convective terms [38]. The second-order central difference scheme is adopted for the diffusion terms. A hybrid of the Adams-Bashforth and Crank-Nicolson schemes is adopted for time advancement of the gas-phase integration [39]. A fully implicit scheme is adopted for the integration of the volume-fraction equation of the sand phase. Due to the strong coupling of drag forces between the two phases, the two-phase partial elimination algorithm is used to decouple the drag [40].

Conditions for Simulation
The simulation is conducted in a rectangular domain, as shown in Figure 1a. The coordinates X, Y and Z denote the streamwise, vertical and spanwise directions, respectively. In the X and Z directions, the mesh is uniformly distributed with mesh cell sizes of dx = Lx/Nx = 0.0042 m and dz = Lz/Nz = 0.001563 m, respectively. The periodicity condition is assumed hereby. In the Y direction, the mesh is refined near the bottom, with the mesh cell size of dy 1 = 2.3606 × 10 −4 m adopted. The symmetrical boundary conditions are set for both the gas and sand phases at the top of the computational domain: the vertical component of the wind and sand velocities as well as the normal gradients for other variables, including the concentration of gas and sand phases, are assumed to be zero in the Y direction at the top. The nonslip condition and the Werner-Wengle wall function are assumed for the bottom surface [41]. As shown in Equation (10), the boundary conditions are applied to the wall-tangential velocities and the granular pseudotemperature of the sand particles: where Ψ = 0.01 and e w = 0.45 are the specularity coefficient and the particle-wall collision restitution coefficient, respectively [29]. A logarithmic wind profile is applied for the initial condition, as shown in Equation (11): where κ = 0.41 is the von Karman constant, Y 0 = d s /30 the roughness length, and u * is the friction velocity. With regard to the LES simulation, a random fluctuating part of the velocities (u g ',v g ',w g ') is superimposed on the logarithmic profile to generate turbulent flow in the boundary layer. The initial vertical and span wise components of wind velocity are set to v g ' and w g ', respectively. At the initial time t = 0, both the gas pressure and the granular pseudotemperature of the sand are set to zero. An initial sand phase field is adopted by using a random distribution of sand concentration, the total concentration of sand is reasonably set as 1 × 10 −4 in the model. The environmental temperature is maintained at a value of 15 • C. When using the Courant-Friedrichs-Lewy (CFL) to determine the limited range of time step, the maximum Courant number can be estimated as CFL max = ∆tmax u g /dx + v g /dy + w g /dz ≤ 0.3.
Then, ∆t ≤ 0.3/max u g /dx + v g /dy + w g /dz ≈ 8.7 × 10 −5 s. In order to guarantee the precision of the LES turbulence simulation, the iteration time step can be reasonably set to a value smaller than 8.7 × 10 −5 s. On the premise of satisfying CFL criterion, after comprehensive consideration, the final decision to set the time step ∆t for both phases as 1.0 × 10 −5 s. . In order to guarantee the precision of the LES turbulence simulation, the iteration time step can be reasonably set to a value smaller than 8.7 × 10 −5 s. On the premise of satisfying CFL criterion, after comprehensive consideration, the final decision to set the time step ∆t for both phases as 1.0 × 10 −5 s. In the Eulerian-Eulerian granular pseudofluid model, the drag coefficient is a key parameter for simulating the momentum transfer between two phases [42]. For a given slip velocity of 1.50 m/s and a particle diameter of 110 μm, the distribution of the drag coefficient β as a function of the volume fraction of the sand phase αs is given in Figure 1b, which shows that the drag coefficients for three classic models are similar to each other. Therefore, the expression of the drag coefficient by Di Felice is selected for the present simulation [35]. To compare with the simulation, three wind tunnel experiments are introduced with friction velocity u as the unknown [43][44][45]. Equation (11) is used to estimate the values for u in experiments based on half the size of the wind tunnel height and its corresponding freestream velocity. Then, the freestream winds in the simulation are determined so that values of u are produced that are as close to those of in experiments as possible. One of the values for u at ds ≈ 250 μm is shown in Table 1.   In the Eulerian-Eulerian granular pseudofluid model, the drag coefficient is a key parameter for simulating the momentum transfer between two phases [42]. For a given slip velocity of 1.50 m/s and a particle diameter of 110 µm, the distribution of the drag coefficient β as a function of the volume fraction of the sand phase α s is given in Figure 1b, which shows that the drag coefficients for three classic models are similar to each other. Therefore, the expression of the drag coefficient by Di Felice is selected for the present simulation [35].
To compare with the simulation, three wind tunnel experiments are introduced with friction velocity u * as the unknown [43][44][45]. Equation (11) is used to estimate the values for u * in experiments based on half the size of the wind tunnel height and its corresponding freestream velocity. Then, the freestream winds in the simulation are determined so that values of u * are produced that are as close to those of in experiments as possible. One of the values for u * at d s ≈ 250 µm is shown in Table 1.    Figure 2 shows a typical view of the distribution of the sand concentration α s , which is averaged spatially in the Z direction at each time step and then averaged temporally over the entire computation period. It turns out that the sands are inhomogeneously distributed and more than half of their mass is concentrated within a height taking up 20% of the entire thickness of the near-surface sand flow, which is similar to earlier simulations and observation [22,36,46]. The maximum sand concentration is on the order of 10 −5 , which is consistent with Liu and slightly lower than that observed by Zhang et al. [47,48].  Figure 2 shows a typical view of the distribution of the sand concentration αs, which is averaged spatially in the Z direction at each time step and then averaged temporally over the entire computation period. It turns out that the sands are inhomogeneously distributed and more than half of their mass is concentrated within a height taking up 20% of the entire thickness of the near-surface sand flow, which is similar to earlier simulations and observation [22,36,46]. The maximum sand concentration is on the order of 10 −5 , which is consistent with Liu and slightly lower than that observed by Zhang et al. [47,48].

Sand Velocity
The streamwise sand velocity us and mass flux q are averaged spatially along two periodic directions at each time step and then averaged temporally over the entire computation period. Figure  3 shows the simulated height profiles of the streamwise sand velocity. To optimize the comparison, the ensemble average value us and Y are normalized by the freestream wind velocity U0 and the thickness of the boundary layer Yb, respectively. In the modeling, the thickness of the boundary layer Yb is set as 0.10m. This height is basically a representative of the maximum saltating height of blown sand particles, which is consistent with the height measured by the wind tunnel experiments [44,47]. First, as shown in Figure 3a, sand velocity varies logarithmically above a certain height (Y/Yb = 0.2), which is qualitatively consistent with the two experiments and simulations but different from the linear pattern obtained by using a particle image velocimetry (PIV) measurement [10,22,49]. At a given height, the sand velocity (or the interphase momentum transfer) increases with friction velocity. As the height increases, the variation of the sand velocity becomes less sensitive to friction velocity, which can be attributed to the presence of fewer obstacles to the full acceleration of the sand particles in the more "dilute" region. On this issue, the simulation agrees qualitatively well with the two experiments. Figure 3b shows a specific comparison between the height profiles of the present simulation and those of the two experiments that are chosen from the results depicted in Figure 3a. The pattern obtained by the simulation agrees very well with the power function-like result of Dong, which was achieved with a PIV measurement [44]. This is because the PIV technique calculates the Eulerian

Sand Velocity
The streamwise sand velocity u s and mass flux q are averaged spatially along two periodic directions at each time step and then averaged temporally over the entire computation period. Figure 3 shows the simulated height profiles of the streamwise sand velocity. To optimize the comparison, the ensemble average value u s and Y are normalized by the freestream wind velocity U 0 and the thickness of the boundary layer Y b , respectively. In the modeling, the thickness of the boundary layer Y b is set as 0.10 m. This height is basically a representative of the maximum saltating height of blown sand particles, which is consistent with the height measured by the wind tunnel experiments [44,47]. First, as shown in Figure 3a, sand velocity varies logarithmically above a certain height (Y/Y b = 0.2), which is qualitatively consistent with the two experiments and simulations but different from the linear pattern obtained by using a particle image velocimetry (PIV) measurement [10,22,49]. At a given height, the sand velocity (or the interphase momentum transfer) increases with friction velocity. As the height increases, the variation of the sand velocity becomes less sensitive to friction velocity, which can be attributed to the presence of fewer obstacles to the full acceleration of the sand particles in the more "dilute" region. On this issue, the simulation agrees qualitatively well with the two experiments. Figure 3b shows a specific comparison between the height profiles of the present simulation and those of the two experiments that are chosen from the results depicted in Figure 3a. The pattern obtained by the simulation agrees very well with the power function-like result of Dong, which was achieved with a PIV measurement [44]. This is because the PIV technique calculates the Eulerian velocity field by correlating interframe windows that contain continuous grayscale distribution, which is similar to the Eulerian-Eulerian simulation that regards the sand phase as a continuum. The difference between their magnitudes is predictable because the values for the friction velocity of the experiment are rather crudely estimated. Only above a certain height, the logarithmic pattern is shared by the simulation and the result obtained from the phase Doppler particle anemometer (PDPA) measurement by Kang et al. [45]. When approaching the sand bed surface, however, the result from the PDPA technique deviates from the results of the simulation. This is probably because as a point measurement, the PDPA technique is sensitive to the interferences of interparticle motions in the 'dense' region.
sand samples in experiment are mixing sands with nonspherical shapes, which inevitably lead to a vertical shift of the turning point from that in simulation where particles are treated as ideally spherical and with a single diameter. Moreover, as the height increases, the difference between profiles becomes greater, indicating that the variation of the sand velocity becomes more sensitive to particle size. This feature is shared by Dong et al. [44] in the upper layer. Generally, it is not easy for the PIV technique to provide satisfactory near-wall velocity information either because the strong light reflection on the boundary deteriorates the image quality necessary for interframe correlation, or because the size of the interrogation window has to be large enough to guarantee the accuracy of the correlation, which in turn sacrifices the spatial resolution.  Figure 3c shows the height profile of the ensemble average streamwise sand velocity varying with particle size. There is a turning point (Y/Y b = 0.2 for the simulation) that divides the profile into two parts. In the upper part, the normalized sand velocity increases with decreasing particle size, which is in contrast to the lower part. This result can be explained as follows: the finer the sand particle, the more easily it is accelerated by wind in the "dilute" region where the momentum transfer from wind to sand is dominant; in contrast, when approaching the sand bed surface, a larger particle becomes superior in retaining momentum after interparticle collisions in the "dense" region. Figure 3c indicates that Dong el al (2006) did not approach down enough to obtain their own turning point, below which we would share the same pattern that particles with bigger diameter suspended lower than the smaller ones. The difference in the height of the turning point is probably due to the fact that sand samples in experiment are mixing sands with nonspherical shapes, which inevitably lead to a vertical shift of the turning point from that in simulation where particles are treated as ideally spherical and with a single diameter. Moreover, as the height increases, the difference between profiles becomes greater, indicating that the variation of the sand velocity becomes more sensitive to particle size. This feature is shared by Dong et al. [44] in the upper layer. Generally, it is not easy for the PIV technique to provide satisfactory near-wall velocity information either because the strong light reflection on the boundary deteriorates the image quality necessary for interframe correlation, or because the size of the interrogation window has to be large enough to guarantee the accuracy of the correlation, which in turn sacrifices the spatial resolution.

Sand Flux
The height profile of the sand mass flux can be defined in light of the sand concentration and streamwise sand velocity [10,44,50], as shown in Equation (12): Similar to the definition in Section 4.2, the angle bracket in the above equation is defined as the ensemble average treatment for the quantity. Figure 4 shows the simulated height profile of q(Y) normalized by the reference value q(Y 0 ) at Y 0 = 2 × 10 −3 m, as previously obtained by Dong et al. [43,44]. One can conclude that (1) all the flux profiles are made up of three parts. The first part is a near-surface layer in which the flux increases with height. The second part refers to the upper layer in which the flux decays rapidly. The third part is a peak of flux around Y/Y 0.1 = 0.2 that joins the two layers. Part (2) of the conclusion is that at a given height, the flux increases with the wind velocity shown in Figure 4a and with particle size shown in Figure 4b. These two features can be found in Dong as measured with a segmented sand sampler and Dong as measured by the PIV technique and calculated with Equation (12) [43,44].

Sand Flux
The height profile of the sand mass flux can be defined in light of the sand concentration and streamwise sand velocity [10,44,50], as shown in Equation (12): Similar to the definition in Section 4.2, the angle bracket in the above equation is defined as the ensemble average treatment for the quantity. Figure 4 shows the simulated height profile of q(Y) normalized by the reference value q(Y0) at Y0 = 2 × 10 −3 m, as previously obtained by Dong et al. [43,44]. One can conclude that (1) all the flux profiles are made up of three parts. The first part is a nearsurface layer in which the flux increases with height. The second part refers to the upper layer in which the flux decays rapidly. The third part is a peak of flux around Y/Y0.1 = 0.2 that joins the two layers. Part (2) of the conclusion is that at a given height, the flux increases with the wind velocity shown in Figure 4a and with particle size shown in Figure 4b. These two features can be found in Dong as measured with a segmented sand sampler and Dong as measured by the PIV technique and calculated with Equation (12) [43,44].
One can also conclude from Figure 4 that the results from Dong and the simulation agree well in the upper layer, but they deviate from each other towards the sand bed surface [43]. This is probably because the segmented sand sampler used in Dong is an intrusive tool that has a drag effect on the sand flow [43]. Such an effect becomes especially evident near the surface where a wind flow with limited strength propels a large number of sand particles. As a result, certain particles obtain momenta and flying heights smaller than they are expected to be, thus resulting in a drop in q above the surface (Y/Y0.1 = 0.4 in this case) and an overestimation of q in the near-surface layer. A nonintrusive measurement technique such as PIV, on the other hand, has to finish two steps to obtain q according to Equation (12). Moreover, it obtains only the relative size of q, because the absolute value of the sand concentration cannot be obtained directly from grayscale images [44].

Fluctuation of the Sand Phase
In the present numerical model, the pseudotemperature θ directly reflects the characteristics of fluctuation of sand movement. This "temperature" is irrelevant to the thermal one, but is an index of the complex interaction between the viscous and inertial forces asserted on sand particle. Aiming to solve the pseudofluid modelling of granular kinematics, Mathiesen proposed the relation between the turbulent RMS (root mean square) velocity of particle and pseudotemperature: = √< 2 >= √3 , in which cs 2 = us' 2 + vs' 2 + ws' 2 ), that is, the sum of the squares of components of particle fluctuating velocity in streamwise, vertical and spanwise directions [17]. Based on the isotropic assumption us' 2 = vs' 2 = ws' 2 that is reasonable when particles do not conglomerate, one can obtain cs 2  = 3us' 2 . In all, the turbulent RMS particle velocity is given as = √< ′ 2 >= √ . One can also conclude from Figure 4 that the results from Dong and the simulation agree well in the upper layer, but they deviate from each other towards the sand bed surface [43]. This is probably because the segmented sand sampler used in Dong is an intrusive tool that has a drag effect on the sand flow [43]. Such an effect becomes especially evident near the surface where a wind flow with limited strength propels a large number of sand particles. As a result, certain particles obtain Atmosphere 2020, 11, 448 9 of 14 momenta and flying heights smaller than they are expected to be, thus resulting in a drop in q above the surface (Y/Y 0.1 = 0.4 in this case) and an overestimation of q in the near-surface layer. A nonintrusive measurement technique such as PIV, on the other hand, has to finish two steps to obtain q according to Equation (12). Moreover, it obtains only the relative size of q, because the absolute value of the sand concentration cannot be obtained directly from grayscale images [44].

Fluctuation of the Sand Phase
In the present numerical model, the pseudotemperature θ directly reflects the characteristics of fluctuation of sand movement. This "temperature" is irrelevant to the thermal one, but is an index of the complex interaction between the viscous and inertial forces asserted on sand particle. Aiming to solve the pseudofluid modelling of granular kinematics, Mathiesen proposed the relation between the turbulent RMS (root mean square) velocity of particle and pseudotemperature: , that is, the sum of the squares of components of particle fluctuating velocity in streamwise, vertical and spanwise directions [17]. Based on the isotropic assumption u s ' 2 = v s ' 2 = w s ' 2 that is reasonable when particles do not conglomerate, one can obtain c s 2 = 3 u s ' 2 .
In all, the turbulent RMS particle velocity is given as u rms = u 2 s = √ θ. Figure 5 shows the variation of the spatially and temporally averaged u rms with height by the present simulation and by the PDPA measurement of Kang et al. [45], the latter of which is by far the most detailed experimental data about the turbulent RMS velocity of aeolian sand particles ever found. The friction velocity in experiment is deduced by collectively using the logarithmic law and the mainstream wind. It can be seen that, within the saltation layer, the fluctuation of particle movement increases along with the height and reaches the peak at Y = 0.005 m. As Y > 0.005 m, it drops with the increase of height. Near the sand bed surface, the sand concentration is relatively high, and the above-mentioned fluctuation is supposed to be mainly determined by the intensive interparticle collision. Far away from the saltation layer, on the contrary, the concentration of saltation is relatively low, and the fluctuation can only be determined by the drag of air flow. Now the distribution patterns by simulation and experiment qualitatively agree with each other that, the fluctuation of sand movement near the surface is much higher than that high up in the air. Therefore, the interparticle collision must be the dominating factor accounting for the fluctuation of sand particle movement. It is noteworthy that, close to surface (Y = 0 m) the fluctuation drastically reduces to zero, indicating that the sand movement is inhibited strongly by the sand bed. As shown in Figure 5a, the magnitude of fluctuation increases along with the friction velocity. This is because higher wind energy generates more intensive interparticle colliding activity, thus strengthening the fluctuation of sand particle movement. As shown in Figure 5b, the magnitude of fluctuation increases along with the particle diameter. This is because larger particles result in a higher probability of interparticle collision (or a greater chance of "overlapping in space"). However, on top of the saltation layer, sand particles are so sparse that the above-mentioned probability is quite low. Therefore, either the wind strength or the particle size has very limited influence on the fluctuation.
The reasons for the quantitative discrepancy between simulation and experiment in Figure 5 can be concluded as follows. (1) In simulation the particle diameter is uniform, while for experiment there is a distribution for particle diameter (usually not as symmetric as some classical models like the normal distribution), so the mean diameter is different between simulation and experiment. Furthermore, as concluded above, particle diameter is key to the magnitude of fluctuation. (2) In simulation there are some simplifications, such as the isotropic assumption for the fluctuating components and the assumption to exclude large-scale fluctuation. In experiment the measuring result contains both the large-scale and the small-scale fluctuations, although the small-scale components take the majority [17,51]. (3) The simulation is conducted in space while the PIV measurement is realized on a plane, the latter of which may be a projection of the original 3D information. The reasons for the quantitative discrepancy between simulation and experiment in Figure 5 can be concluded as follows. (1) In simulation the particle diameter is uniform, while for experiment there is a distribution for particle diameter (usually not as symmetric as some classical models like the normal distribution), so the mean diameter is different between simulation and experiment. Furthermore, as concluded above, particle diameter is key to the magnitude of fluctuation. (2) In simulation there are some simplifications, such as the isotropic assumption for the fluctuating components and the assumption to exclude large-scale fluctuation. In experiment the measuring result contains both the large-scale and the small-scale fluctuations, although the small-scale components take the majority [17,51]. (3) The simulation is conducted in space while the PIV measurement is realized on a plane, the latter of which may be a projection of the original 3D information.

Drag between the Two Phases
The drag is the essential connection of energy transfer between the two phases, affecting significantly the characteristics of sand movement. The average drag that the air phase passed on the sand phase is Fdrag = β∆U∆x∆y∆z, which represents the overall drag in a unit volume. fdrag = β∆U is called the drag density, where ∆U = √( − ) 2 + ( − ) 2 +( − ) 2 is defined as the slip velocity between the two phases. Figure 6 shows the variation of spatially and temporally averaged drag density with height by the present simulation, for which an experimental comparison is very difficult to find. However, according to the above results and discussions, since the velocity, flux and fluctuation of sand phase have been proven to be comparable with those of the experiment, it is sufficient to believe that the simulation result of drag density is instructive (at least qualitatively). It can be seen that, as the height goes up, the drag density increases drastically in the start and reaches the maximum value at Y = 0.005 m. Then it decreases gradually, and finally reduces to zero at the top surface of saltation layer (Y = 0.1 m). Near the sand bed surface, sand particles are propelled by the wind and start to move. Due to the inertia, however, the velocity difference between the two phases is so huge that fdrag has a great potential to increase and finally reach the peak. As the height goes up, there is a continuous momentum transfer between the two phases, so the velocity difference shrinks, leading to the decrease of fdrag. At the sand bed surface, fdrag = 0 due to the nonslip effect. As shown in Figure 1b, with the slip velocity and particle diameter fixed, the drag coefficient β increases along with the sand concentration. Therefore, distribution of fdrag shown in Figure 6 is consistent with that of the sand concentration. As shown in Figure 6a, the slip velocity increases along with the friction velocity, thus lifting the value of fdrag. As shown in Figure 6b, smaller particle leaves a larger room for the drag of

Drag between the Two Phases
The drag is the essential connection of energy transfer between the two phases, affecting significantly the characteristics of sand movement. The average drag that the air phase passed on the sand phase is F drag = β∆U∆x∆y∆z, which represents the overall drag in a unit volume. f drag = β∆U is called the drag density, where ∆U = u g − u s is defined as the slip velocity between the two phases. Figure 6 shows the variation of spatially and temporally averaged drag density with height by the present simulation, for which an experimental comparison is very difficult to find. However, according to the above results and discussions, since the velocity, flux and fluctuation of sand phase have been proven to be comparable with those of the experiment, it is sufficient to believe that the simulation result of drag density is instructive (at least qualitatively). It can be seen that, as the height goes up, the drag density increases drastically in the start and reaches the maximum value at Y = 0.005 m. Then it decreases gradually, and finally reduces to zero at the top surface of saltation layer (Y = 0.1 m). Near the sand bed surface, sand particles are propelled by the wind and start to move. Due to the inertia, however, the velocity difference between the two phases is so huge that f drag has a great potential to increase and finally reach the peak. As the height goes up, there is a continuous momentum transfer between the two phases, so the velocity difference shrinks, leading to the decrease of f drag . At the sand bed surface, f drag = 0 due to the nonslip effect. As shown in Figure 1b, with the slip velocity and particle diameter fixed, the drag coefficient β increases along with the sand concentration. Therefore, distribution of f drag shown in Figure 6 is consistent with that of the sand concentration. As shown in Figure 6a, the slip velocity increases along with the friction velocity, thus lifting the value of f drag . As shown in Figure 6b, smaller particle leaves a larger room for the drag of air flow, so f drag increases along with the decrease of particle size. This is especially obvious at the lower part of saltation layer where the sand concentration is higher.

Interparticle Collision
The interparticle collision restitution coefficient e is a representative of the elasticity of interparticle collision and defined as the ratio of kinetic energies after and before collision. Du et al. suggest that the pseudotemperature of particle θ is closely related to the interparticle collision coefficient [52]. Therefore, e is critical to the quantitative simulation of fluid driven particle movement. As a matter of fact, due to the shape irregularity of natural quartz sand particles, it is ineffective to use the spherical particles as the substitute in the experiments to measure e. As a result, it is impossible to have a general solution for e among different experiments and simulations. Wang et al. suggest e ∈ [0.428, 0.455] by conducting a wind tunnel experiment using natural sands [11], while Haff and Anderson agree on e ∈ [0.6, 0.8] [14]. Kang and Guo adopted e = 0.87 and e = 0.89 in their DEM simulation of aeolian sand flow [36], while Kang and Liu used e = 0.7 in their simulation [49]. In all, the present work is not looking for a universally effective value (or range) for e, but to investigate the influence of e on the characteristics of fluctuation of particle movement.
Atmosphere 2020, 11, x FOR PEER REVIEW 11 of 15 air flow, so fdrag increases along with the decrease of particle size. This is especially obvious at the lower part of saltation layer where the sand concentration is higher.

Interparticle Collision
The interparticle collision restitution coefficient e is a representative of the elasticity of interparticle collision and defined as the ratio of kinetic energies after and before collision. Du et al. suggest that the pseudotemperature of particle θ is closely related to the interparticle collision coefficient [52]. Therefore, e is critical to the quantitative simulation of fluid driven particle movement. As a matter of fact, due to the shape irregularity of natural quartz sand particles, it is ineffective to use the spherical particles as the substitute in the experiments to measure e. As a result, it is impossible to have a general solution for e among different experiments and simulations. Wang et al. suggest e  [0.428, 0.455] by conducting a wind tunnel experiment using natural sands [11], while Haff and Anderson agree on e  [0.6, 0.8] [14]. Kang and Guo adopted e = 0.87 and e = 0.89 in their DEM simulation of aeolian sand flow [36], while Kang and Liu used e = 0.7 in their simulation [49]. In all, the present work is not looking for a universally effective value (or range) for e, but to investigate the influence of e on the characteristics of fluctuation of particle movement. Figure 7 shows the variation of the profiles of urms and interparticle collision frequency fc with e. Here, the interparticle collision frequency is based on a kinetic theory of granular flow, which can be calculated by combining the pseudotemperature and sand phase concentration [26]: This is because as the particles save more energy after the collision, there is a larger opportunity for the subsequent interparticle collision to happen, thus leading to a stronger fluctuation of sand particle movement. The influence of e on the two profiles is more obvious near the sand bed surface than up in the air, which is due to the fact that with a higher sand concentration, the effect of interparticle collision will be magnified to a greater extent. From Figure 7a one can also notice that with e changing from 0.6 to 0.8, the peak RMS velocity varies from 1.124 to 1.555 m/s, corresponding to an increase ratio of 38.3%. Therefore, the fluctuation of particle movement is very sensitive to the interparticle collision restitution coefficient. It is noted that since the experimental verification for the simulation results in Figure 7 cannot be obtained at present and the e value used in some modeling studies based on DEM is not within the range of [0.6, 0.8] [36]. Therefore, it is speculated that the classical statement such as e  [0.6, 0.8] mentioned above may not be an optimal condition for the quantitative simulation of turbulent aeolian sand transport.  Here, the interparticle collision frequency is based on a kinetic theory of granular flow, which can be calculated by combining the pseudotemperature and sand phase concentration [26]: f c = 6.77α s θ 1/2 / d s 1 − (α s /α s,max ) 1/3 . It can be seen that, u rms and f c increases along with e. This is because as the particles save more energy after the collision, there is a larger opportunity for the subsequent interparticle collision to happen, thus leading to a stronger fluctuation of sand particle movement. The influence of e on the two profiles is more obvious near the sand bed surface than up in the air, which is due to the fact that with a higher sand concentration, the effect of interparticle collision will be magnified to a greater extent. From Figure 7a one can also notice that with e changing from 0.6 to 0.8, the peak RMS velocity varies from 1.124 to 1.555 m/s, corresponding to an increase ratio of 38.3%. Therefore, the fluctuation of particle movement is very sensitive to the interparticle collision restitution coefficient. It is noted that since the experimental verification for the simulation results in Figure 7 cannot be obtained at present and the e value used in some modeling studies based on DEM is not within the range of [0.6, 0.8] [36]. Therefore, it is speculated that the classical statement such as e ∈ [0.6, 0.8] mentioned above may not be an optimal condition for the quantitative simulation of turbulent aeolian sand transport.

Conclusions
In this study, the aeolian sand transport within a turbulent atmospheric boundary layer was simulated using a Eulerian-Eulerian model coupling the granular pseudofluid optimization and the gas-phase large eddy simulation. This continuum model is efficient in simulating a dense particle flow and makes it possible to avoid determining the specific modes of near-surface particle motion,

Conclusions
In this study, the aeolian sand transport within a turbulent atmospheric boundary layer was simulated using a Eulerian-Eulerian model coupling the granular pseudofluid optimization and the gas-phase large eddy simulation. This continuum model is efficient in simulating a dense particle flow and makes it possible to avoid determining the specific modes of near-surface particle motion, i.e., saltation or creep, for which distinct laws of motion must be applied [1,2].
The simulated profiles of streamwise mean sand velocity and sand mass flux qualitatively agree with the reported experiments. The discrepancy between them occurs near the sand bed surface. This highlights the potential of the present model in addressing the near-surface sand transport layer, which is problematic for nonintrusive measurement techniques with limited spatial resolution. The simulated profiles of turbulent RMS particle velocity indicate that the interparticle collision is mainly responsible for the fluctuation of sand particle movement. Therefore, the interparticle collision restitution coefficient should be carefully chosen in simulation. The above conclusions may be helpful as reference in developing the Euler-Euler fluid driving particle transport. To fully achieve geophysical significance, the next step is to simulate the sand flow in nature rather than in a wind tunnel. For this effort, it is necessary to develop more reasonable constitutive relations to deal with a wider range of particle sizes, and it is also necessary to reproduce the fluid turbulence as appropriate for the typical spectra of wild gusts.