Scaling Laws in Aeolian Sand Transport Under Low Sand Availability

Previous studies of wind‐blown sand have considered either fully erodible or non‐erodible soils, but the transport over sparsely sand‐covered soils is still poorly understood. The quantitative modeling of this transport is important for parameterizing Aeolian processes under low sand availability. Here we show, by means of numerical simulations, that the sand transport rate Q scales with the wind shear velocity u∗ as Q=a⋅1+b⋅u∗/u∗t−1⋅d/g⋅ρf⋅u∗2−u∗t2 $Q=a\cdot \left[1+b\cdot \left({u}_{\ast }/{u}_{\ast \mathrm{t}}-1\right)\right]\cdot \sqrt{d/g}\cdot {\rho }_{\mathrm{f}}\cdot \left({u}_{\ast }^{2}-{u}_{\ast \mathrm{t}}^{2}\right)$ , where u∗t is the minimal threshold u∗ for sustained transport, d is particle size, g is gravity and ρf is air density, while u∗t and the empirical parameters a and b depend on the sand cover thickness. Our model explains the transition from the quadratic to cubic scaling of Q with u∗ as soil conditions change from fully erodible to rigid and provides constraints for modeling Aeolian transport under low sand availability.

KAMATH ET AL. 10.1029/2022GL097767 2 of 9 for sustained transport (Pähtz & Durán, 2020;Ralaiarisoa et al., 2020). However, natural Aeolian systems encompass a broad range of soil types characterized by low sand availability on the ground, including bare and crusted soils sparsely covered with mobile sediments (Amir et al., 2014;Shao, 2008). The characteristics of Aeolian transport over such types of soil, that is, when the thickness of the mobile sand layer on the rigid ground is comparable to a few grain diameters, are poorly understood.
The quantitative understanding of these characteristics is important for various fields, in particular to improving wind-blown sand and dust schemes in climate models. Once in the atmospheric circulation, dust substantially affects the planet's climate and biosphere, atmospheric geochemistry, the hydrological cycle, and various other components of the Earth system, yet estimates of vertical dust flux and atmospheric dust budget are counted amongst the largest uncertainty sources in climate simulations (Kok et al., 2012;Mahowald et al., 2014;Schepanski, 2018;Shao, 2008). Since dust is rarely entrained directly by wind but is, instead, emitted mainly by the impacts of wind-blown sand grains onto the ground (Shao et al., 1993), an accurate model for the Aeolian sand transport rates over various types of soil, from fully erodible to fully non-erodible, is required. However, it is difficult to derive such a model from analytical computations alone, given the broad range of natural soil erodibility conditions associated with sparsely covered bare, gravel and crusted soils (Amir et al., 2014;Macpherson et al., 2008;Shao, 2008;Wang et al., 2011).
Therefore, here we perform the direct computation of grain trajectories during Aeolian sand transport by means of particle-based simulations, or Discrete-Element-Method (DEM). This type of simulation has been applied previously to investigate Aeolian transport over fully erodible beds (Carneiro et al., 2011;Comola et al., 2019;Durán et al., 2012), thereby introducing a helpful means to elucidate processes that are difficult to assess in wind tunnel or field experiments, such as the mechanisms of sediment transport very close to the bed. Indeed, using DEM simulations, it is possible to resolve these mechanisms, as well as their impact on the resulting sand flux, without any need for assumptions about the splash process, the rebound dynamics or the modification of the wind profile in the transport layer-which are rather directly computed. As we discuss in the subsequent sections, our DEM simulations show that the scaling of the sand flux with u * displays considerable and yet unreported dependence on the availability of sand on the ground-characterized here through the thickness of the mobile sediment layer covering the non-erodible surface.

Numerical Experiments
The DEM consists of solving Newton's equations of motion for all particles in the system under consideration of the main forces acting on them (Cundall & Strack, 1979). In contrast to other types of numerical models of soil erosion (Almeida et al., 2008;Anderson & Haff, 1988;Kok & Renno, 2009), DEM models of Aeolian transport do not rely, thus, on a splash function to represent the ejection of particles from the soil owing to grain-bed collisions. Rather, the lift-off velocities of the rebound and ejected particles are obtained by directly solving their equations of motion under consideration of particle-particle interactions (Lämmel et al., 2017;Yin et al., 2021). We start our simulations by pouring sand-sized spherical particles of diameter d uniformly distributed in the range 160 ≤ d/μm ≤ 240 onto a flat horizontal rigid bed at the bottom of the simulation domain-which has dimensions (L x × L y × L z )/d m = (200 × 8 × 1,000), with d m = 200 μm denoting the mean grain size (Figure 1). In doing so, we generate a thin bed of N p randomly poured particles on the ground, where the bed thickness δ 0 is determined by N p . For instance, N p = 30, 000 for δ 0 ≈ 15 d m .
Furthermore, we adopt periodic boundary conditions in the along-wind (x) and cross-wind (y) directions and impose a reflective horizontal wall at the top of the simulation domain, to avoid that particles escape through crossing the upper boundary at z = L z . However, we find that removing this reflective wall would allow only few particles for escaping, thus leading to a negligible change in the results of our simulations.
Once the particles come to rest and the bed has been formed, a few particles are injected into the simulation domain to impact on the ground, thus producing a splash and ejecting grains into air. The Aeolian drag force on the particles is computed with the expression, KAMATH ET AL.

10.1029/2022GL097767
3 of 9 where ρ f = 1.225 kg/m 3 is the air density and r = − ( ) is the difference between the velocity v i of particle i and the wind velocity u(z i ) at the height z i of the particle's center of mass. Furthermore, r = | r | , while the drag coefficient d is computed trough (Cheng, 1997) , where the Reynolds number Re = f r ∕ , with μ = 1.8702 ×10 −5 kg m −1 s −1 denoting the dynamic viscosity of the air.
The wind velocity profile is constant along x and y throughout the simulations, while the initial vertical profile of the horizontal (downstream) wind velocity, u x (z), reads, where u * is the wind shear velocity, κ = 0.4 the von Kármán constant, z 0 ≈ d m /30 is the roughness of the quiescent bed, and h 0 is the bed height, which is set as the uppermost height within the granular surface where the particles move with velocity smaller than 0.1 u * (Carneiro et al., 2011). However, the acceleration of the particles owing to the action of the drag force extracts momentum from the air (Anderson & Haff, 1988;Owen, 1964), thus leading to a modification of the wind velocity profile. The modified velocity profile is obtained by numerical integration of (Carneiro et al., 2011) where τ p (z) is the grain-borne shear stress and is given by with  d ( ) denoting the horizontal component of the total drag force on the particles with center of mass at Z j , while A = L x ⋅ L y (Carneiro et al., 2011).
Furthermore, in order to obtain a rough rigid bed underneath the mobile sand cover, we deposit the mobile particles on top of a sheet of "frozen" immobile particles as displayed in Figure 1 (see Supporting Information S1 for the set of DEM particle-particle contact force equations, including the presence of the frozen particles). In doing so, the rigid bed provides a model for a fully consolidated dune surface or bare granular surface, where the constituent immobile particles have the same diameter as the mobile grain size.

Results and Discussion
Once transport begins, some of the grains composing the initial bed layer are entrained into flow, so that the bed layer thickness-which has initial value δ 0 at time t = 0-decreases over time until transport eventually achieves steady state. At steady state, the bed layer thickness amounts to , where C b ≈ 0.02 is an empirical parameter and Θ(x) denotes the Heaviside function, that is, Θ(x) = 1 if x ≥ 0 and Θ(x) = 0 if x < 0 (see Figure S3 in Supporting Information S1). Therefore, the term b * ∕ √ m (≲0.5 for all scenarios) denotes the thickness of the total eroded layer, relative to the particle size, from the beginning of transport until steady state.
We note that periodic boundary conditions are applied in our simulations (see Section 2), so that the number of particles in the system is constant over time (Carneiro et al., 2011;Durán et al., 2012;Pähtz & Durán, 2020). Indeed, the domain of our simulations may be interpreted as a small stretch of soil over which the sediment flux is in the steady state. Due to fluctuations associated with the transport dynamics, the difference between the particle mass outflux from and influx into this soil stretch varies over time, but on average, the total number of particles within the associated volume is constant over time.
We begin our discussion by considering an initial bed thickness δ 0 = 15d m , for which we observe steady-state transport conditions (δ s ≈ 14.8 d m ) consistent with the fully erodible bed scenario reported in previous studies. Specifically, our simulations reproduce quantitatively the height-integrated, non-suspended mass flux of transported particles, Q, as a function of u * over fully erodible beds, and the observation that, for moderate wind conditions (u * /u * t ≲ 4), Q is approximately proportional to τ − τ t , with = f 2 * denoting the mean shear stress of the turbulent wind flow over the surface, and t = f 2 * t corresponding to the minimal threshold τ for transport ( Figure 2). Furthermore, our numerical predictions match the experimental observations of the nearly exponential decay of the vertical particle concentration with the height above the ground and the value of u * t ≈ 0.165 m/s predicted for the mean particle size in our simulations (see Figure S1 in Supporting Information S1).
However, as we decrease the initial bed layer thickness δ 0 substantially, we observe a change in the scaling of the steady-state sediment flux with u * . More precisely, our simulation results follow, approximately, the model. 10.1029/2022GL097767 5 of 9 where u * t,∞ ≈ 0.165 m/s, a ∞ ≈ 22.15 and b ∞ ≈ 5.28 denote the values of u * t and the empirical constants a and b, respectively, associated with fully erodible bed scenario (δ s /d m → ∞), while the best fits to the simulation data in the range δ s /d m ≤ 10 yield C t ≈ 0.14, c t ≈ 0.83, C a ≈ 0.47, c a ≈ 0.76 and c b ≈ 2.61 ( Figure 2).
Wind tunnel experiments (Ho et al., 2011) revealed a cubic scaling of Q with u * on fully rigid beds. Here, we find that sediment transport rates over a soil that is not fully rigid but contains, instead, a thin layer of mobile sediment, further depends on this layer's thickness according to Equations 5-8. Specifically, the coefficient b in Equation 8 controls the transition from the cubic to the quadratic scaling of Q with u * in Equation 5 as bed conditions change from fully rigid (δ s = 0) to fully erodible (δ s ≫ d). Moreover, while the coefficient a provides an attenuating factor for Q near the rigid bed scenario, a decrease in bed thickness reduces the minimal threshold shear velocity, u * t , as we elucidate next.
To shed light on the microscopic origin of Equation 5, we note that momentum conservation yields (Bagnold, 1941;Ho et al., 2011;Sørensen, 2004), where ℓ hop denotes the mean hop length of the saltating particles, while u 0↓ and u 0↑ are their mean horizontal impact and lift-off velocities, respectively. Furthermore, ℓ hop and u 0↓ − u 0↑ (computed as explained in Section 4 of the Supporting Information S1) are related to the mean horizontal grain velocity u 0 = (u 0↓ + u 0↑ )/2 (or slip velocity) through the approximate scaling expressions hop ∝ 2 0 ∕ and u 0↓ − u 0↑ ∝ u 0 (Ho et al., 2011), which leads to ≈ ⋅ ( 0∕ ) ⋅ [ − t ] , where C u is an empirical parameter.
An increase in u * over a fully erodible bed leads to an enhancement of the particle concentration in the transport layer without significantly affecting u 0 , so that Q scales quadratically with u * in the fully erodible bed regime (Ho et al., 2011). By contrast, the transport layer over the hard surface is, for a given saltation flux, much thicker than over an erodible bed because of the non-saturated feedback which keeps a larger wind velocity in the saltation layer (Ho et al., 2011). The weak coupling between the particles and the wind in the transport layer over a fully non-erodible surface results in a linear scaling of u 0 with u * , thus yielding a cubic scaling of Q with u * in the fully rigid bed regime (Ho et al., 2011).
Here we find that, in the presence of a thin layer of mobile sand on the hard ground, the scaling of u 0 with u * further depends on δ s (Figure 3c). We find that with C u ≈ 1.68, where the RHS of Equation 9 is the multiplicative factor of [τ − τ t ] in Equation 5, that is, including the values of u * t , a and b estimated from Figure 2. Therefore, Equation 9 elucidates the microscopic origin of Equation 5. Since all scenarios (δ 0 , u * ) considered here are associated with saturated transport conditions in the steady state (see Figure S3 in Supporting Information S1), that is, since the total mass of particles in the transport layer under given u * − u * t is the same for all values of δ 0 considered, the effect of sand availability on the scaling of Q(u * ) is attributed entirely to the dependence of u 0 on this availability, encoded in the parameters on the RHS of Equation 9. Our simulations further show that, as sand availability decreases and the transport layer expands, transport can be sustained at increasingly lower u * (Figure 2a and Equation 6). This finding is further consistent with the wind-tunnel observation that u * t over fully rigid beds is lower than over fully erodible beds (Ho et al., 2011).
To the best of our knowledge, our study is the first one to estimate sediment transport rates from direct numerical simulations of particle trajectories under intermediate soil erodibility conditions between fully erodible and fully non-erodible. We find that our results remain approximately valid when the rigid bed underneath the mobile sediment layer is a smooth flat surface. However, the immobile roughness elements on the hard ground have a crucial effect on the value of the Aeolian sand flux.

of 9
In the regime where saltating particles collide onto a sand bed of thickness ≲2d m , and in the presence of roughness elements on the hard ground underneath, sand particles are ejected through splash events mainly backwards, that is, the majority of ejecta displays negative horizontal lift-off velocity component. This result can be understood by noting that, as downwind hopping grains impact obliquely upon the thin sand layer covering the rough ground, they mobilize soil grains forward, which, however, collide with the roughness elements located in their front. Upon such collisions, the trajectories of the bed particles mobilized by grain-bed impacts are reflected backwards, as elucidated through our granular splash experiments (Figure 4, where N ej is the number of ejected grains per impact).
These dynamics, which act by attenuating Q upon exposure of the bed roughness elements, are encoded in the coefficient a in Equation 5, and constitute behavior opposite to the effect of the bed thickness on b and u * t , which  contribute to enhancing Q (see Equations 5-8). These competing effects lead to an anomaly in the dependence of Q on the bed thickness, with the emergence of a minimum around δ 0 ≈ 2 d m (or δ s ≈ 1.8 d m ). This anomaly is not observed when the ground is a smooth flat surface ( Figure 5). The bed thickness associated with the minimum Q is independent of u * , thus indicating that the anomaly reported here is purely a signature of the bed roughness and is not affected by the flow properties.
We note that, notwithstanding the strong decrease of N ej with the bed thickness in the regime δ 0 /d m ≳ 2 (Figure 4b), the steady-state sand flux Q in this regime is only weakly affected by the amount of mobile grains on the ground (Figure 2a). Therefore, our simulation results are providing evidence in support of the hypothesis that the magnitude of Q is controlled by the rebound dynamics of sand grains during transport-as assumed, for instance, in a recent purely rebound-based model (Pähtz et al., 2021)-rather than by the splash process. Our results further help to elucidate the observation that cohesion, which affects mainly the splash process by enhancing particle-particle attractive interaction forces within the bed, has little impact on Q and the threshold for Aeolian transport cessation, as these are mainly controlled by rebound dynamics (Comola et al., 2019).
Our model reproduces the scaling laws of Q with u * observed experimentally over fully erodible and rigid beds (Figures 2 and S1 in Supporting Information S1). However, various ingredients that are essential to improve the quantitative assessment of Aeolian sand flux, such as complex particle geometric shapes and aerodynamic entrainment (Li et al., 2020), should be incorporated in future work. Furthermore, we have employed sand-sized non-erodible roughness elements, but natural soils encompass much broader particle size distributions, including gravels, pebbles and rocks. From our results, we expect that such coarser non-erodible elements have even larger impact on the sand flux scaling. Our model is paving the way toward a quantitative representation of sand availability conditions in larger scale models, such as regional Earth system simulations, by explicitly incorporating the information of local steady-state bed thickness in the parameterization of Aeolian sand transport rates.
Previous work developed continuum models for Aeolian flux that explicitly account for sand supply and spatio-temporal variations in bed surface properties, including moisture, shells, non-erodible elements and vegetation (De Vries et al., 2014;Hoonhout & Vries, 2016). Furthermore, the particle-based simulations adopted in the present work provide a means to improve our understanding of the (microscopic) particle-scale mechanisms controlling the response of Aeolian transport processes to different types of soil and particle-bed interactions. Future research combining insights from both types of model could thus help to achieve improved numerical simulations of Aeolian soil morphodynamic processes at different scales (Durán et al., 2010;Kroy et al., 2002;Werner, 1995), by incorporating the effect of sediment availability on sediment flux and erosion/deposition rates.

Conclusions
In conclusion, we have presented the first numerical model for wind-blown sand flux under low sand availability, by characterizing this flux as a function of the thickness of the mobile sediment layer available for transport on the ground. Specifically, we showed that the Aeolian sand flux scales with the excess shear stress multiplied by a coefficient that decreases with the mobile layer thickness covering the non-erodible ground, thereby yielding a model for Aeolian transport rates under intermediate bed erodibility conditions between the fully erodible and fully non-erodible scenarios. Our model elucidates how the scaling of the Aeolian sand flux Q with the wind shear velocity u * changes from quadratic to cubic as bed conditions change from fully erodible to fully non-erodible, respectively (Ho et al., 2011).
We also found that the roughness elements on the rigid bed affect the sediment flux upon rigid bed exposure, by causing an anomaly in the behavior of Q with the bed layer thickness, with the occurrence of a minimum which is independent on the flow conditions. These findings will have an implication for the representation of non-erodible elements associated with different types of soil in future experimental and theoretical studies.

Data Availability Statement
All data included in this work are generated from our numerical model and is available online (https://doi. org/10.6084/m9.figshare.19469501). The data for validation with experiments is available from (Creyssels et al., 2009).