Anomalous Scaling of Aeolian Sand Transport Reveals Coupling to Bed Rheology

Predicting transport rates of windblown sand is a central problem in aeolian research, with implications for climate, environmental, and planetary sciences. Though studied since the 1930s, the underlying many-body dynamics is still incompletely understood, as underscored by the recent empirical discovery of an unexpected third-root scaling in the particle-fluid density ratio. Here, by means of grain-scale simulations and analytical modeling, we elucidate how a complex coupling between grain-bed collisions and granular creep within the sand bed yields a dilatancy-enhanced bed erodibility. Our minimal saltation model robustly predicts both the observed scaling and a new undersaturated steady transport state that we confirm by simulations for rarefied atmospheres.

Figure 1

Data from laboratory measurements [7, 23, 35] and previous [3] as well as our original sand-transport simulations, based on the discrete element method (DEM) [36, 37] (see Supplemental Material [38] for details), obey the transport-rate scaling in Eqs. (1a) and (1b) (solid line). The DEM simulations allow us to toggle between a complex boundary-layer wind velocity profile (dots and circles) and simplified “fully rough” flow conditions (squares) based on Prandtl’s turbulent closure [42], cf. Eq. (2c), and to study a wide range of particle-fluid density ratios $s$, terminal grain settling velocities $v_s$, and shear stresses $\tau$ in excess of the transport threshold ${\tau }_t$. The dashed line amounts to neglecting midair grain collisions.

Figure 2

Sand transport rate scaling predicted by our minimal two-species saltation model without midair collisions [$Q=(\tau -{\tau }_t)V$] and with a static-bed splash function [9] for terminal grain settling velocities $v_s=\{{10}^{3/2},{10}^2,{10}^{5/2},{10}^3\}$ (circles, squares, diamonds, stars) and particle-fluid density ratios $s=\{4^0,\dots ,4^6\}{v}_s^2/10$ (colors). Small (large) $s$ tend to be on the right (left).

Figure 3

Quasi-two-dimensional DEM-based sand-transport simulations (as in Ref. [3]) of the creep and dilatancy regime of aeolian transport. (a) Granular creep visualized by the average height-resolved horizontal grain velocity $⟨v_x⟩(z)$ in the sediment bed ($z<0$, below the elevation at which high-energy grain-bed collisions occur [4, 38]). Its increase with height $z$ and imposed wind shear stress $\tau$ (solid lines) reveals a characteristic skin depth $\lambda \approx 0.72$ (dashed lines). (b) Because of progressive smoothing, the granular volume fraction profiles $\phi (z)$ (solid lines) around $z=0$ deviate considerably from the limiting form for saltation on a quiescent bed—roughly a step from $\phi \approx 0.58$ to the exponential extrapolations of $\phi (z>0)$ (dashed lines) [23, 35]. They exhibit a focal point ${\phi }_f=\phi (z\approx -\lambda )\approx 0.1$. (c) The constitutive relation $\mu (\mathrm{\Delta },I)$ (solid line) for aeolian creep at subyield conditions ($\mu \lesssim 0.3$) is similar to that of other sheared granular flows [46, 47]. It interconnects the local friction coefficient $\mu =-{\sigma }_{xz}^c/{\sigma }_{zz}^c$, local normalized streamwise velocity fluctuations $\mathrm{\Delta }\equiv (-T_{xx}/{\sigma }_{zz}^c{)}^{1/2}$ with $T_{xx}\equiv ⟨v_x^2⟩-⟨v_x{⟩}^2$, and local inertial number $I\equiv ⟨d{v}_x/dz⟩/\sqrt{-{\sigma }_{zz}^c}$. Here, ${\sigma }_{xz}^c$ (${\sigma }_{zz}^c$) is the shear (normal) component of the structural granular stress associated with grain-grain contacts in the bed [38]. (d) The value $\phi (z=0)$ is taken as a proxy for the number of grains available for splash ejection in Eq. (3), and its ${\tau }_g(0)$ dependence motivates Eq. (4) with ${\tau }_Y\approx 0.17$.

Figure 4

Laboratory measurements, DEM-based sand-transport simulations (cf. Fig. 1), and predictions by our minimal saltation model with cooperative splash according to Eqs. (3) and (4) (inset) collapse on a master curve defined by (a) Eqs. (1a) and (5) (approximately $V\propto s^{1/3}$), corresponding to saturated transport conditions, or (b) an undersaturated steady state [Eqs. (1a) and (6), approximately $V\propto s^{1/6}$, upper inset]. Depending on the initial condition, this state can also be reached and sustained in DEM simulations, based on the code of Ref. [3] (open black circles) or Ref. [37] (open green squares) for $s\gtrsim {10}^5$, regardless of the driving flow velocity profile (cf. Fig. 1). Lower inset: exemplary transition between the steady states, as occasionally spotted in the simulations. Solid (dashed) lines correspond to Eqs. (1a) and (1b) with (without) the term representing midair collisions, which are neglected in our minimal saltation model. Filled symbols as in Figs. 1 and 2.