Nontrivial scaling in supply limited Aeolian sand transport

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 the parametrization of Aeolian processes under low sand availability. Here we show, by means of particle-based numerical simulations, that the Aeolian sand transport rate Q scales with the wind shear velocity u* as Q = a.[1 + b . (u*/u*t - 1)] .[?](d/g) . ρ f. (u* ² - u*t2), 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


Introduction
Aeolian (wind-blown) sand transport produces ripples and dunes and plays a vital role in shaping the Earth's surface.This transport occurs mainly through sand grains hopping along the surface (saltation), thereby transferring to the ground momentum that may set new particles into hopping, rolling or sliding motion (Bagnold, 1941;Shao, 2008;Kok et al., 2012).Furthermore, the particle splash generated by saltating grains provides one main mechanism of dust aerosol emission (Gillette, 1981;Shao et al., 1993), which has major feedbacks with the biosphere, the hydrological cycle and various other components of the Earth system (Mahowald et al., 2014;Schepanski, 2018).The accurate modeling of windblown sand is, thus, important for the development of reliable geomorphodynamic, climate and Earth system models (Shao, 2008).Indeed, previous models of Aeolian transport focused mainly on the transport over either fully erodible beds, such as migrating dunes and ripples (Anderson & Haff, 1988;Shao & Li, 1999;Sauermann et al., 2001;Almeida et al., 2008;Kok & Renno, 2009;Lämmel et al., 2012;Pähtz et al., 2014;Comola et al., 2019), or rigid, fully non-erodible beds, such as consolidated dunes and bare soils (Ho et al., 2011).These studies have shown that wind-blown transport rates follow either a quadratic or a cubic scaling with the wind shear velocity u * -which is proportional to the mean flow velocity gradient in turbulent boundary layer flow -depending upon the bed being fully erodible or fully non-erodible, respectively (Creyssels et al., 2009;Ho et al., 2011).Moreover, a quartic scaling of the sand flux with u * , characterizing a collisional or intense transport regime where the saltation layer is connected to the granular bed through an intermediate granular layer of intense mid-air collisions, has been reported for fully erodible bed conditions when u * exceeds about 4u * t , where u * t stands for the minimal threshold 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 (Shao, 2008;Amir et al., 2014).The characteristics of Aeolian transport over such types of soil, i.e., 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 (Shao, 2008;Kok et al., 2012;Mahowald et al., 2014;Schepanski, 2018).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 (Shao, 2008;Amir et al., 2014;Macpherson et al., 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;Durán et al., 2012;Comola et al., 2019), 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 Discrete-Element-Method 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 (Anderson & Haff, 1988;Almeida et al., 2008;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).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, where ρ f = 1.225 kg/m 3 is the air density and v r i = v i − u(z i ) 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 i = |v r i |, while the drag coefficient , where the Reynolds number Re i = ρ f υ r i d i /µ, 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 (Owen, 1964;Anderson & Haff, 1988), 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 F d x (Z j ) 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 Fig. 1 (see Suppl.Mat. 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 = 0decreases 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, i.e., Θ(x) = 1 if x ≥ 0 and Θ(x) = 0 if x < 0 (see Fig. S3).Therefore, the term C b u * /( √ gd 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 = 15 d 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 u 2 * denoting the mean shear stress of the turbulent wind flow over the surface, and τ t = ρ f u 2 * t corresponding to the minimal threshold τ for transport (Fig. 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 Suppl.Mat., Fig. 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, 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 (Fig. 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 Eqs. ( 5)-( 8).Specifically, the coefficient b in Eq. ( 8) controls the transition from the cubic to the quadratic scaling of Q with u * in Eq. ( 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.Error bars denote the standard deviation from averaging over 5 s within the steady state.
To shed light on the microscopic origin of Eq. ( 5), we note that momentum conservation (Bagnold, 1941;Sørensen, 2004;Ho et al., 2011), 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 Suppl.Mat.) are related to the mean horizontal grain velocity u 0 = (u 0↓ + u 0↑ )/2 (or slip velocity) through the approximate scaling expressions hop ∝ u 2 0 /g and u 0↓ − u 0↑ ∝ u 0 (Ho et al., 2011), which leads to 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).
manuscript submitted to Geophysical Research Letters 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 (Fig. 3c).We find that with C u ≈ 1.68, where the RHS of Eq. ( 9) is the multiplicative factor of [τ − τ t ] in Eq. ( 5), i.e., including the values of u * t , a and b estimated from Fig. 2. Therefore, Eq. ( 9) elucidates the microscopic origin of Eq. ( 5).Since all scenarios (δ 0 , u * ) considered here are associated with saturated transport conditions in the steady state (see Suppl.Mat., Fig. S3), i.e., 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 Eq. ( 9).Our simulations further show that, as sand availability decreases and the transport layer expands, transport can be sustained at increasingly lower u * (Fig. 2a and Eq. ( 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.
In the regime where saltating particles collide onto a sand bed of thickness 2 d m , and in the presence of roughness elements on the hard ground underneath, sand particles are  ejected through splash events mainly backwards, i.e., 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 (Fig. 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 Eq. ( 5), and constitute behavior opposite to the effect of the bed thickness on b and u * t , which contribute to enhancing Q (see Eqs. ( 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 (Fig. 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 (Fig. 4b), the steady-state sand flux Q in this regime is only weakly affected by the amount of mobile grains on the ground (Fig. 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 (Figs. 2 and 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 (Werner, 1995;Kroy et al., 2002;Durán et al., 2010), 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.

Open Research
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).

Figures S1 to S4
Table S1 Introduction In this Supplemental Material, we briefly review the features of the Discrete-Element-Method referred to in the main document, including the complete set of the equations of motion, the details of the numerical integration of these equations, and the models of particle-particle interactions adopted in the simulations of Aeolian sand transport.Furthermore, we present the results of our numerical simulations performed to verify our model, the vertical profiles of the wind velocity and grain-borne shear stress during steady-state transport, and the behavior of the transport layer thickness as a function of the thickness of mobile sand layer, as mentioned in the main document.

S1 Discrete-Element-Method
In the Discrete-Element-Method, the equations of motion are solved for every particle in the system under consideration of the main forces acting on them.These forces are, in the process of non-suspended Aeolian transport of cohesionless particles, the drag force, the inter-particle contact forces and the gravitational force.

S1.1 Equations of motion and contact force model for the sand particles
The equation of translational motion for a particle of mass m i at position r i reads, where F d i is the drag force on particle i, computed with the model described in the main document, g is gravity, N p is the number of particles in the system, j denotes the index of a neighbouring particle that is in contact with particle i, and F c ij denotes the contact force exerted by particle j on i (with Contact between particles j and i occurs with their center-to-center distance is smaller than the sum of their radii, i.e., the contact force acts only if the particles overlap.To model the Corresponding author: Sandesh Kamath, skamath@uni-koeln.decontact force, the following equation is used to define the overlap, where d i and d j are the diameters of particles i and j, respectively, r ij = r i − r j , with r j standing for the position of particle j, and e ij,n = r ij /r ij denotes the normal unit vector pointing from the center of particle j to the center of particle i, with r ij = |r ij |. There are various contact force models for application in DEM simulations, and the modelling of these forces is still an active matter of research (Cundall & Strack, 1979;Schäfer et al., 1996;Brilliantov et al., 1996;Silbert et al., 2001;Di Renzo & Di Maio, 2004;Pöschel & Schwager, 2005;Kruggel-Emden et al., 2007;Luding, 2008;Machado et al., 2012;Parteli et al., 2014;Fan et al., 2017;Schmidt et al., 2020;Santos et al., 2020).In our simulations, we adopt the linear spring-dashpot model, because this model has been employed in previous simulations of wind-blown sand that reproduced the scaling laws associated with Aeolian transport over fully erodible beds (Carneiro et al., 2011(Carneiro et al., , 2013;;Durán et al., 2012;Comola et al., 2019).
Specifically, F c ij can be described as the sum of a normal component, F c ij,n , and a tangential component, F c ij,t .Each of these components encodes an elastic term and a dissipative term, while the magnitude of the tangential force is bounded by the Coulomb friction criterion.
The equations for F c ij,n and F c ij,t read (Cundall & Strack, 1979;Silbert et al., 2001;Santos et al., 2020) where m eff = m i m j /(m i + m j ), with m i and m j denoting the masses of particles i and j, respectively, k n , k t , γ n , γ t and µ s are model parameters, discussed in Section S1.3 below, while the relative normal velocity v ij,n and the relative tangential velocity v ij,t between particles i and j are computed via with v ij = v i − v j denoting the difference between the velocities of particles i and j (v i and v j , respectively), and ω i and ω j standing for their respective rotational velocities.
Moreover, in Eq. ( 4), ξ ij,t is the tangential displacement accumulated as the particles are in contact.The displacement is set as zero at initiation of the contact and is computed in the reference frame of the rotating particle pair to compensate for the effect of rigid body rotations, as described in detail in previous work (Silbert et al., 2001;Santos et al., 2020).
The equation of rotational motion for particle i reads with I i = m i d 2 i /10 and ω i denoting the moment of inertia and the angular velocity of particle i, respectively, and M ij corresponding to the torque on particle i associated with The contact forces between mobile sand particles and the rigid particles constituting the roughness elements of the bed are computed using the same model as in the previous section, but considering that the rigid particles have an infinite mass (Verbücheln et al., 2015).Specifically, the normal and tangential components of the contact force from a rigid particle j on a mobile particle i are computed with Eqs. ( 3) and ( 4), respectively, by setting m eff = m i .Furthermore, contact forces between rigid particles are not considered.

S1.3 Model parameters
Table S1 displays the values of the parameters in Eqs. ( 3) and (4), i.e., the elastic constants k n and k t , the damping coefficients γ n and γ t , and the Coulomb friction coefficient, µ s .The elastic and damping constants are taken from previous models for Aeolian sand transport over fully erodible beds (Carneiro et al., 2011(Carneiro et al., , 2013;;Comola et al., 2019).In particular, the elastic constant for normal contact, k n , is estimated using where d m = 200 µm is the mean particle size adopted in our simulations, while Y = 1 MPa is the Young's modulus adopted in previous work (Carneiro et al., 2011;Comola et al., 2019) and in our computations.Furthermore, for the elastic constant for tangential constant, we use k t = k n /3, while the friction coefficient is consistent with values adopted previously (Comola et al., 2019).

S1.4 Numerical implementation and particle-wind coupling
To solve the equations of motion of the granular phase, we employ LAMMPS (Largescale Atomic/Molecular Massively Parallel Simulator), which is an open source DEM solver based on MPI implementation (Plimpton, 1995).Furthermore, we have extended this solver to incorporate the hydrodynamic description of the turbulent wind flow over the granular surface, developed in previous work (Carneiro et al., 2011;Durán et al., 2012) and briefly reviewed in the main document.To this end, we have included new modules (LAMMPS "fixes") into the granular package of the DEM solver, to set the initial (logarithmic) vertical profile of the mean horizontal wind speed, to compute the drag force, and to update this drag force and the wind profile owing to the process of momentum exchange between the particles and the wind.These modules are available from the corresponding author upon request.
To verify our numerical simulations, we compare our numerical predictions for the heightintegrated mass flux (Q) of wind-blown particles over a fully erodible bed with corresponding wind-tunnel observations of this flux as a function of the wind shear velocity, u * .We compute Q using the following equation, where N is the number of particles in the system, m i and v x i denote the mass and horizontal speed of the i−th particle, respectively, and A = L x • L y is the horizontal area of the simulation domain.To measure this flux, we start the wind-blown transport process as described in the main document and wait until this transport achieves steady state.The typical (physical) time separating the begin of wind-blown transport from the steady state is about 2-3 seconds and independent of u * (Carneiro et al., 2011;Durán et al., 2012;Pähtz et al., 2014;Comola et al., 2019).All results reported in the present work refer to the characteristics of steady-state wind-blown transport, and denote mean quantities obtained from averaging over about 5-10 seconds during steady-state transport.
Furthermore, by suitably normalizing the steady-state flux Q, we obtain the following nondimensional quantity,  We see in Fig. S1 that our numerical predictions for Q(Θ) (circles) agree quantitatively well with observations from wind-tunnel experiments (Creyssels et al., 2009), denoted by the stars.The best fit to our simulation results using Q = aΘ + b yields a ≈ 0.5 and b ≈ 0.0026 (dashed line in Fig. S1), from which obtain the minimal threshold Θ t ≈ 0.0064 below which no transport occurs ( Q = 0).From Eq. ( 10), this value of Θ t leads to the minimal threshold wind shear velocity for sustained transport, u * t ≈ 0.165 m/s.
We note that the value of u * t predicted from our simulations is consistent with the prediction that u * t is about 80% of the minimal threshold wind shear velocity u * ft required to initiate transport, with A ft ≈ 0.1 (Bagnold, 1941;Shao & Lu, 2000).Indeed, by applying the mean particle size d m = 200 µm of our simulations in Eq. ( 11), we obtain u * ft ≈ 0.206 m/s, i.e., our model is consistent with the relation u * t ≈ 0.8 u * ft predicted for wind-blown transport.

S3
The modified wind profiles for varying δ 0 The initial vertical profile of the horizontal downstream wind velocity u x is logarithmic and follows Eq. ( 3) of the main document.However, this wind velocity profile is updated every time-step, since the acceleration of the grains extracts momentum from the air thus creating a negative feedback on the wind.The modification of the wind velocity profile is computed using Eqs.( 4) and ( 5

S4 Computation of the mean hop length and the mean horizontal impact and lift-off velocities
To compute the mean hop length hop and the average horizontal impact and lift-off velocities, u 0↓ and u 0↑ , respectively, we consider only the grains with a minimum vertical lift-off velocity of √ 6gd, i.e., the grains that achieve a minimum height of 3 d m above the bed height.The values of hop , u 0↓ and u 0↑ , are then averaged over a time window of 5 seconds during steady-state sand transport.
The mean horizontal grain velocity or slip velocity u 0 is then computed using Furthermore, to obtain the mean hop length, we start from the mean hop time, which is given by (Ho et al., 2011), where v 0 ↑ is the mean ascending vertical velocity of the grains (also averaged over 5 seconds in the steady state).Furthermore, the horizontal acceleration of the grains is given by, so that the mean hop length is approximated as, S5 Relation between the steady-state bed layer thickness, δ s , and the initial bed layer thickness, δ 0 As mentioned in the main document, once the sand transport process begins, the initial thickness of the mobile sand bed applied in the numerical simulations, δ 0 , decreases toward a smaller value δ s , which is achieved when transport conditions have reached steady state.
As depicted in Fig. S3, δ s and δ 0 are linearly related to each other, and the difference between both values of bed thickness displays a slight increase with u * owing to the effect of wind shear velocity on enhancing erosion.However, we find that the scaling laws reported in the main document are valid whatever value of bed thickness is chosen, while the following relation applies,  In Fig. S3, the filled symbols correspond to numerical simulations in which the wind is carrying the maximum possible number of particles, i.e., the flux is saturated.Starting with δ 0 = 15 d m , for instance, and under a given value of u * − u * t , we observe no change in the average number of particles in the Aeolian layer (or, equivalently, the mass density of dragged particles) upon a decrease in the initial bed thickness δ 0 , as long as the scenarios associated with the filled symbols in Fig. S3 are considered.However, the empty symbols in this figure constitute scenarios where the bed thickness is so small, that the wind flow does not dispose of enough particles on the ground to drive transport toward the saturated flux.These empty symbols are associated with a value of steady-state bed thickness equal to 0 and are referred to as under-saturated.Specifically, for these empty symbols, the wind eroded the entire sand bed and still the amount of sand transport is not enough to saturate the sand flux (the most extreme, non-vanishing flux scenario of such under-saturated regime in our simulations would be, in particular, the case of one single grain hopping downwind).
We have thus not considered these under-saturated scenarios in our analysis of Q(u * ) in the main document, since we are interested here in an expression for the saturated flux that accounts for the bed erodibility.Moreover, it is a straighforward conclusion that, in the under-saturated regime, the mass density of the transport layer decreases upon a reduction of the initial bed thickness δ 0 /d m .Nevertheless, we note that under-saturated transport scenarios constitute an interesting topic to be investigated in future work.For instance, such scenarios have applications to areas that are devoid of any sand availability but subjected to an upwind flux that is under-saturated (for instance in the presence of upwind vegetation or moisture), or over inter-dune bedrock areas within fields of sparsely distributed dunes (Fryberger et al., 1984)).

S6 Transport layer thickness as a function of the bed thickness
As explained in the main document, the transport layer expands gradually as the soil erodibility conditions change from fully erodible to rigid.To quantify this process, we compute the characteristic length-scale l ν associated with the nearly exponential decay of the particle concentration ν(z) with the height z above the ground, i.e., ν(z) = ν 0 exp(−z/l ν ) (17 where ν 0 is the particle concentration extrapolated to the bed (z = 0).Fig. S4 shows the behavior of l ν as a function of the thickness of mobile sand layer on the ground, δ 0 .We see that, when the rigid ground is flat, l ν decreases monotonically as δ 0 increases toward 15 d m (the fully erodible bed scenario).However, when the rigid surface is armoured with nonerodible elements (which in our simulations have the same size as the mobile particles), a minimum in l ν is observed near δ 0 = 2 d m .This minimum can be explained by the prevailing occurrence of backward ejecta in the range δ 0 2 d m (as described in the main document), which leads to lower values of l ν in this range when non-erodible particles cover the ground, compared to the flat ground scenario.As δ 0 becomes much larger than the particle size, the effect of the immobile particles on the Aeolian transport thickness becomes negligible, and l ν approaches asymptotically the value corresponding to the fully erodible bed as shown in 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 × 1000), with d m = 200 µm denoting the mean grain size (Fig.

Figure 1 :
Figure 1: (a) Snapshot of the numerical experiment at t = 0, indicating the dimensions of the simulation domain and the undisturbed wind profile.(b) Side-view of an excerpt of the sediment bed, displaying a layer of mobile particles (blue) of thickness δ 0 on top of the immobile particles constituting the rough ground.

Figure 2 :
Figure 2: (a) Sand flux Q rescaled with the excess shear stress, τ − τ t , plotted as a function of (u * − u * t ) for different values of the initial bed thickness, δ 0 ; inset: the minimal threshold shear velocity for sustained transport, u * t as a function of the steady-state bed thickness, δ s .(b) Circles and squares denote the parameters a and b in Eq. (5), respectively, as obtained from the best fit to the data in (a).The continuous lines in (a) and (b) denote the best fits using Eqs.(6)-(8) in the range δ s /d m ≤ 10 (the continuation of these fits toward larger δ s /d m or fully erodible bed scenario is indicated by the dashed line as a guide to the eye).Error bars denote the standard deviation from averaging over 5 s within the steady state.

Figure 3 :
Figure 3: (a) Mean hop length, hop , and (b) difference between the mean grain horizontal velocities at impact and lift-off, u 0↓ − u 0↑ , as a function of the slip velocity u 0 .The dashed lines in (a) and (b) denote hop ≈ 0.065u 2 0 and u 0↓ − u 0↑ ≈ 0.43u 0 , respectively, obtained from the best fits to the simulation data.In (c), the slip velocity is shown as a function of u * − u * t for different values of δ 0 .The legend in (c) applies as well to both (a) and (b).

Figure 4 :
Figure 4: By means of granular splash numerical experiments with impact angles and velocities characteristic of wind-blown sand transport (a), we find that most ejected grains have negative horizontal lift-off velocity, when the value of the bed layer thickness is 2 d m , and positive otherwise (b).The snapshots correspond to a simulation using a bed layer thickness ≈ 2 d m .Most of the mobile (blue) particles lying on the rigid grains (red) have been rendered transparent for better visualization of the splashed particles.

Figure 5 :
Figure 5: Sand flux Q as a function of δ 0 , obtained with u * = 0.30 m/s.We considered the non-erodible surface consisting of a smooth flat ground (blue) and immobile particles (red).
which we plot in Fig.S1as a function of the Shields number,Θ = u 2 * ρ f (ρ p − ρ f )gd m (10)where ρ p = 2650 kg/m 3 and ρ f = 1.225 kg/m 3 denote the densities of the particles and the air, respectively, while d m = 200 µm is the mean particle diameter and g = 9.81 m/s 2 is gravity.

Figure S1 :
Figure S1: Normalized steady-state flux Q as a function of the Shields number Θ, considering a fully erodible bed (δ 0 ≈ 15 d m ).
) in the main document.The vertical profiles of the modified wind velocity u x and the grain-borne shear stress τ p are shown for different values of the mobile layer thickness δ 0 in Fig. S2.

Figure
Figure S2: (a) The modified wind profiles for different values of δ 0 alongside the initial logarithmic profile; (b) grain-borne shear stress profile as a function of the height above the bed for different δ 0 .The results were obtained with u * = 0.30 m/s.
16) where C b ≈ 0.02 is an empirical parameter and Θ(x) denotes the Heaviside function, i.e., Θ(x) = 1 if x ≥ 0 and Θ(x) = 0 if x < 0. Therefore, the term C b u * /( √ gd m ) denotes the thickness of the total eroded layer, relative to the particle size, from the beginning of transport until steady state (i.e., as the bed thickness evolves from δ 0 toward δ s ).

Figure
Figure S3: (a) Values of the bed layer thickness at steady-state transport, δ s , plotted against the initial values of bed layer thickness, δ 0 , for different values of the wind shear velocity u * .The plot in (b) denotes a zoom into the region of bed layer thickness comparable to the particle size.Filled symbols correspond to saturated transport conditions, while empty symbols denote under-saturated scenarios (the same color code used for the filled symbols in the legend applies to specify u * in these empty symbol scenarios).

Figure S4 :
Figure S4: Thickness l ν of the transport layer as a function of the thickness δ 0 of the mobile sediment layer covering the non-erodible surface, which consists of a flat horizontal surface (blue symbols) and a sheet of immobile particles (red symbols).The results were obtained with u * = 0.30 m/s.

Table S1 :
Parameters of the Discrete-Element-Method.