Aeolian Creep Transport: Theory and Experiment

Aeolian sand transport drives geophysical phenomena, such as bedform evolution and desertification. Creep plays a crucial, yet poorly understood, role in this process. We present a model for aeolian creep, making quantitative predictions for creep fluxes, which we verify experimentally. We discover that the creep transport rate scales like the Shields number to the power 5/2, clearly different from the laws known for saltation. We derive this 5/2 power scaling law from our theory and confirm it with meticulous wind tunnel experiments. We calculate the creep flux and layer thickness in steady state exactly and for the first time study the relaxation of the flux toward saturation, obtaining an analytic expression for the relaxation time.


Introduction
Aeolian sand transport is responsible for desertification processes and bed topography evolution. Therefore, the prediction of the sand transport rate under turbulent wind conditions is fundamental for the management of land resources and sustainable environmental development (Durán et al., 2011;Kok et al., 2012;Pähtz et al., 2020;Shao, 2008). From a physical point of view, sand transport is the dynamic response of a granular bed under the shearing of turbulent wind. The two main types of aeolian sand transport are saltation (grains move by successive jumps over the bed) and creep (rolling and sliding over the bed surface) (Bagnold, 1941;Pähtz et al., 2015). Saltation has been extensively investigated in the past decades (Andreotti, 2004;Ho et al., 2011;Iversen & Rasmussen, 1999;Martin & Kok, 2017;Sørensen, 2004). In contrast to saltation, creep still remains virtually unexplored.
In previous studies, mainly bed traps were used to measure creep transport rates (Bagnold, 1941;, and only very limited theoretical approaches for creep have been reported up to the present (Wang & Zheng, 2004). All proposed creep transport laws are either empirical or semiempirical, and lack a solid physical basis. Therefore, a solid theoretical framework describing creep precisely is badly required. Here we will propose for the first time a theoretical model for aeolian creep, finding a 5/2 power scaling with the Shields number for the creep transport rate. This new law was validated through wind tunnel experiments. Furthermore, we will explore the relaxation of creep transport to steady state. saltation layer is a dilute rapid granular flow, in which the saltating grains are accelerated by the wind before they impact the creep layer. In the creep layer, a dense shear flow of sand is driven by the shear stress τ 0 applied on the surface of the creep layer.

Saltation Layer
To describe flow in the saltation layer, one can use a continuum saltation model based on mass and momentum conservation (Sauermann et al., 2001) where ρ sal and u sal are the density and velocity of the sand, in the saltation layer, respectively. Г is the erosion rate, and f d and f b denote the drag force and friction force acting on a volume element of the saltation layer. The erosion rate Г, which is the difference between the vertical flux of impacting and splashed grains, can be written as (Parteli et al., 2007;Sauermann et al., 2001) where l is the average saltation length, γ is a model parameter which determines how fast the system reaches its steady state. τ g0 is the grain-borne shear stress at the ground, and τ t is the threshold shear stress. A detailed derivation of the drag and friction forces f d and f b can be found in Parteli et al. (2007) or Sauermann et al. (2001). Finally, combining mass and momentum conservation Equations 1 and 2, the erosion rate Equation 3, and the drag and friction forces f d and f b , the closed model for sand transport in the saltation layer can be established (Durán & Herrmann, 2006;Parteli et al., 2007;Sauermann et al., 2001), and the following expression can be obtained for the steady state saltation flux Q s at saturation.
where z m and z 0 denote the mean saltation height and the roughness height, respectively. z 1 is a reference height, between z 0 and z 1 , α is the splash parameter, and C d is the drag coefficient.
If modifications of the total shear stress by suspended grains are neglected, τ can be assumed to be constant at any height above the ground. According to Owen (1964), the total shear stress τ can be split into the airborne shear stress τ a and the grain-borne shear stress τ g , and we obtain τ = τ a + τ g = const.

Creep Layer
The creep layer is a thin dense sheared sheet of flowing grains driven by the surface shear stress τ 0 . Because the creep layer is dense, the air flow will not directly affect the movement of grains inside the creep layer. At the upper surface of the creep layer ( Figure 1) one has τ a (z = s) ≡ τ a0 and τ g (z = s) ≡ τ g0 . τ a0 is applied directly to the upper surface of the creep layer, and τ g0 is the momentum transfer from the saltating grains to the creep layer during the impacts. Since we are formulating a continuum theory, the discrete impacts from saltating grains are here described by a momentum transfer. Therefore, they together drive the creep layer with τ a0 þ τ g0 ¼ τ 0 ¼ ρ a u 2 * . We define a Cartesian coordinate system Oxz, with the x axis pointing in the streamwise and the z axis in the vertical direction, respectively ( Figure 1). Let s(x, t) be the upper surface of the creep layer and b(x, t) its base. Hence, the thickness of the creep layer is given by h(x, t) = s(x, t) − b(x, t). If variations in the bulk density are neglected, the creep can be considered as incompressible. The continuity equation then becomes The momentum balance equation is where u and w are the velocity components in the x and z directions, respectively. The transverse velocity in y direction is not considered here. σ and τ are the normal stress and internal shear stress, respectively. ρ = ϕρ s is the density where ϕ is the sand volume fraction of the creep layer and ρ s the sand density.
Based on Sayed and Savage (1983), the constitutive relation for the creep layer can be expressed by a nonlinear viscous law which is related to pressure and strain rate In Equation 7, τ y is the Coulombic yield stress which is proportional to normal stress σ, and τ v denotes the viscous stress, which is the non-Newtonian stress caused by particle collisions and is proportional to the square of the shear rate _ γ ¼ ∂u=∂z. d is the grain diameter, χ is a model parameter, g denotes gravity, and μ 0 is the coefficient of dynamic friction. Using kinematic boundary conditions at the upper surface and at the base of the creep layer (Iverson, 2012), we obtain E s and E b are the entrainment rate at the upper surface and the erosion rate at the base of the creep layer, respectively. Additionally, if there is no slip at the base of the creep layer, the dynamic boundary conditions along the upper surface and the base of the creep layer become where the superscripts "s" and "b" denote the upper surface and the base, respectively. It is important to mention that these equations can equally well describe rippled surfaces or rugged terrain just by setting the apropriate boundary condition on b(x, t).

Steady State Solution
In the steady state of creep motion, that is, under uniform and steady conditions, the time and space derivatives in Equations 5 and 6 vanish. In this case, combining the constitutive Equation 7 and the boundary 10.1029/2020GL088644

Geophysical Research Letters
conditions, the exact solution for the equilibrium creep motion can be derived. The velocity profile in the creep layer is given as u where h eq = τ 0 /ρgμ 0 is the depth of the creep layer in steady state. Then, the creep flux in steady state can be obtained as This means that the creep flux increases with the shear stress to the power 5/2, and defining the Shields This is the first time to our knowledge that such a law is proposed for creep transport rate.

General Case
Even not being in steady state, the continuity Equation 5 and the momentum balance Equation 6 can be integrated over z from the base to the surface of the creep layer following Leibniz' rule. The depth-integrated continuity and momentum equations are then obtained using the kinematic and no-slip boundary conditions at the base and the upper surface of Equations 8 and 9 where overbars denote depth-averaged quantities defined as φ ¼ 1=h ð Þ∫ s b φdz for any variable φ. E s is the entrainment rate at the upper surface of the creep layer which is determined by the erosion rate Г of the saltation layer through E s = Γ/ρ s . In order to calculate the three unknown fields, namely, the depth h(x, t), the surface elevation s(x, t), and the average velocity u x; t ð Þ, we need to supplement an independent equation. For this purpose we use the depth-integrated energy flux equation, which was obtained from a weighted residual method (Capart et al., 2015;Richard & Gavrilyuk, 2013;Steffler & Jin, 1993). The weight function selects the velocity u(x, t) using Leibniz' rule and the boundary conditions. After the integration of the momentum equation with this weighting function, the following equation can be obtained ∂ ∂t Assuming that the velocity profile away from steady state has the same form as the velocity profile at steady state, we can obtain the average of the square and the cube of the velocity as u 2 ¼ 25=16 ð Þ u 2 and u 3 ¼ Under these assumptions, the governing Equations 14 and 15 for the creep layer become With the energy flux Equation 15, together with the continuity and momentum Equations 13 and 14, one can therefore obtain the depth h(x, t), the surface elevation s(x, t), and the average velocity u x; t ð Þ for aeolian creep.

Geophysical Research Letters
WANG ET AL.
In principle, the time-dependent equations for the saltation and creep layers could be solved numerically to calculate the transient behavior of aeolian sand transport. In order to simplify Equations 16 and 17, we will assume that saltation has reached steady state, and thus only the creep motion remains in a transient state. In this case, E s = 0, and ∂h=∂x ¼ ∂ u=∂x ¼ 0. Then the governing equations in the creep layer may be written as We have calculated the depth h eq and averaged velocity u eq of steady uniform creep motion in the previous section. In order to reduce the number of parameters and simplify the analysis, we first rewrite Equations 18 and 19 in dimensionless variables given by e h ¼ h=h eq , e u ¼ u= u eq , and e t ¼ t=t Assuming that slope effects are small, that is, The creep model of Equations 20 and 21 is a system of two coupled nonlinear ordinary differential equations. The relaxation process toward the steady state of creep transport can be described by the flow diagram of Figure 2 obtained numerically from Equations 20 and 21. Under different initial conditions, the system trajectories evolve along curved paths away from unstable initial points until they reach the stable Steady State A (1, 1). In steady state, the relationship between thickness and average velocity of the creep layer is given by Equation 10. After normalization, the dimensionless relation for the steady state becomes e u ¼ e h 3=2 , which is the black line shown in Figure 2.
In order to further understand the creep dynamics when the system is not in steady state, we will now consider the transient from one steady state to another. Starting from creeping motion at a lower velocity in the Steady State B (0.5, 0.63), the shear stress τ 0 at the surface of the creep layer is suddenly increased, and consequently, the thickness and average velocity of the creep layer increase to a new Steady State A e h; e u e t→∞ ¼ 1; 1 ð Þ, illustrated by a thick solid blue trajectory in Figure 2. Similarly, the creep flux e q ¼ e u  Figure S1 (see the supporting information). It can be seen that the transient of e h, e u, and e q first increases and then slowly relaxes toward steady values. In addition, we see that e q relaxes slower toward saturation than e u and e h. This means that the height and grain velocity first reach saturation, and the flux is the dominant mechanism controlling the transient to creep saturation, because the relaxation is limited by the slowest process. Note that the corresponding time scales reflect complicated processes involving the whole system dynamics. To find analytically the time scales involved in the transient behavior of creep transport, we search for the time needed to relax to the Steady State Q c given by Equation 11. Using Equations 18 and 19, we derive the relaxation equation for the creep flux We see that the creep flux q approaches the steady state exponentially with characteristic time T c According to Equation 23, T c scales as which characterizes the time scale of the response. T c increases with wind shear velocity as ∼ u 3 * and decreases with grain diameter as ∼ d −1 . To the best of our knowledge, this is the first time that the transient behavior of creep transport has been studied. Further work is required to experimentally confirm this transient behavior.

Experiments
Wind tunnel experiments were carried out to validate our theoretical results. The experiments were conducted in a multifunctional wind tunnel located at the Key Laboratory of Mechanics on Western Disaster

Geophysical Research Letters
and Environment, at Lanzhou University, with an experimental setup similar to that of Zhu et al. (2019). For a schematic design see Figure 3a. The experimental section of the wind tunnel was 22 m long, and its cross section was 1.3 m in width and 1.45 m in height. The incoming flow speed could be continuously adjusted in the range between 4 and 40 m/s. Roughness elements and a wedge were placed in front of the working section to accelerate the development of the boundary layer. A bed of about 10 m in length, 0.6 m in width, and 6 cm in thickness was filled with sand from the Tengger Desert (ρ s = 2,650kg/m 3 ) in China, which was sieved into five different categories (Zhu et al., 2019). In our experiment, we used sand with a median particle diameter of 251 μm. The wind velocity was measured at four different heights above the sand bed (4.3, 8.4, 13.3, and 20.5 cm) with I-type hot-wire probes (DANTEC 55P11) located near the end of the working section, which were connected to a constant-temperature hot-wire anemometer. The wind velocity profile satisfied the classical logarithmic law for turbulent boundary layers: u x z ð Þ ¼ u * κ ln z z 0 , and the wind shear velocity u * was obtained after fitting the logarithmic law to the averaged data (κ is the von Kármán constant of about 0.41). To measure the saltation and creep mass fluxes, a special sand trap located at the downwind end of the working section was implemented. The sand trap consisted of an upper and a lower part (Figure 3b). The upper part was 60 cm in height and 2 cm in width, divided into 30 bins, each bin having a vertical opening of 2 cm, and it was used to collect saltating sand grains. The lower part was used to collect creeping sand grains. To avoid that saltating particles enter the trap for the creeping ones, we designed a special device (Figure 3b). The principle of this device was based on the fact that over a flat surface saltating grains typically impact at angles between 9°and 16° (Bagnold, 1936;Nalpanis et al., 1993;Swann & Sherman, 2013;Willetts & Rice, 1989). To catch these grains, the height of the lower edge of the saltation trap was chosen to be 2 mm above the opening of the creep trap. The creep trap was positioned 15 mm downwind from the saltation trap having its opening at the height of the sand bed (see Figure 3b). Therefore, saltating particles could not enter the creep trap.

Experimental Results and Discussion
As shown in Figure 4, the dimensionless creep transport rate grows with the Shields number, meaning that the creep flux increases rapidly with the wind velocity. The solid line in Figure 4 of Q c = ρ s d ffiffiffiffiffi gd p À Á versus the Shields number Θ = τ 0 /(ρ s gd) was obtained from Equation 11. As shown in Figure 4a, theoretical values from Equation 11 are in quantitative agreement with the experimental data. And the creep transport rate Q c is found to scale as Θ 5/2 , confirming Equation 12. Moreover, in Figure 4b, we compare our theory with previous experimental results . The symbols show the dimensionless creep transport rate obtained from experiments  for d = 0.152 mm and d = 0.382 mm under different Shields numbers, and the colored lines are the theoretical prediction of the form Q c ∝ Θ 5/2 . We observe very good agreement between the experimental data and our theory. Next, we will present also our results for the saltation layer (see Figure S2 in the supporting information). There is excellent agreement between the theory and experiment, indicating that Equation 4 accurately predicts the saltation transport rate. Typically, the sand transport rate can be expressed as a 1.5th order polynomial in the Shields number as seen in Equation 4. Recently also, a linear asymptotic relation between the sand transport rate Q s and the Shields number Θ (or u 2 * ) was suggested (Durán et al., 2012;Ho et al., 2011). As shown in Figure S2, a linear fit of the form Q s ∝ (Θ − Θ t ) is also compatible with the experimental data.
In Figure 5, we see both our experimental and theoretical results for the contribution of creep to the total transport as functions of the Shields number. The theoretical values come from Equations 4 and 11. Figure 5 also compares our theoretical prediction to field measurement (Sherman et al., 2019) and previous wind tunnel data . Sherman et al. (2019) used extensive field observations to find evidence that the proportion of creep transport compared to total transport weakly increases with shear velocity, in agreement with our theory and wind tunnel experiments.

Conclusions
The present study on aeolian creep transport quantifies an important geophysical phenomenon that has consequences on both bedform evolution and desertification. We proposed a theory for creep transport subjected to wind shear with the aid of a constitutive relation and a depth-averaged integration, which remarkably agrees with the results of wind tunnel experiments. Therefore, this theory might both provide a new way to fully and quantitatively validate aeolian creep transport and deliver a complete description of all relevant fields, including those that might be difficult to measure experimentally at present, such as the relaxation process of creep. These results can be useful in numerous applications ranging from sediment transport or aeolian erosive mechanisms to the understanding of the dynamics of the morphology of Earth's surface and on other planets.

Data Availability Statement
All data are available in Zenodo Digital Repository at https://zenodo.org/record/3884458#.Xt3VnoUmRhc.