Bursts in discontinuous Aeolian saltation

Close to the onset of Aeolian particle transport through saltation we find in wind tunnel experiments a regime of discontinuous flux characterized by bursts of activity. Scaling laws are observed in the time delay between each burst and in the measurements of the wind fluctuations at the fluid threshold Shields number θc. The time delay between each burst decreases on average with the increase of the Shields number until sand flux becomes continuous. A numerical model for saltation including the wind-entrainment from the turbulent fluctuations can reproduce these observations and gives insight about their origin. We present here also for the first time measurements showing that with feeding it becomes possible to sustain discontinuous flux even below the fluid threshold.

. Bagnold defined also a dynamic threshold. Once started, saltation will cease at a somewhat lower shear velocity due to momentum transfer from the impacting grains. The ratio between the dynamic and the fluid threshold is 0.8. However, wind shear velocities above the fluid threshold can show temporarily no sand transport. Also Rasmussen and Sorensen observed saltation to be discontinuous below the fluid threshold in the field 4 . Intermittent fluctuations of the turbulent wind produce sand bursts, which can be observed individually if their density is low so that they do not superpose. Induced by surface-layer turbulent structures, the sand bursts evolve downstream into coherent streamers 5,6 . The burst effects of turbulence were first identified in sediment transport in water 7,8 . In aeolian saltation, Stout defined an intermittency function and found that bursts of saltation were approximately normally distributed with respect to the relative wind strength 1 . Schönfeldt statistically correlated the wind gusts in the velocity time series with the saltation events to propose varying dynamic and fluid thresholds. All studies mentioned here were performed in the field, i.e. not under reproducible or fully controlled conditions. Numerical work of Dupont et al. first to demonstrated a correlation between the flow turbulence and the sand streamers using Large Eddy Simulations (LES) 9 . The numerical observation of sand streamers requires large systems and, consequently several assumptions concerning particle interactions to facilitate the computational simulations. No previous works have provided yet an insight of sand bursts at the particle scale.
Here, we report on a systematic investigation of sand bursts not only experimentally observed in a wind tunnel but also numerically reproduced using a stochastic process to model the temporal fluctuations in a turbulent channel flow. The Discrete Element Model (DEM) is a flexible technique to include the main forces acting on a saltating particle and correctly reproduce the particle-particle and fluid-particle iterations. The first one is found in the particle splash and mid air collisions, and the second one is in the particle lift and particle drag 10,11 .

Methods
At Aarhus university, we used a wind tunnel illustrated in Fig. 1. The wind tunnel is 20 m long with 10 m of working section in the downstream end. The rectangular cross section is 0.60 m wide and 0.90 m high. A set of turbulence spires at the entry and a 3 m long array of roughness blocks downwind of the spires have been designed to produce a boundary layer which has approximately the same characteristics as that forming over an infinitely long sand bed near the fluid threshold 12 . Downwind of the array, the bottom of the tunnel is covered by a 15 mm thick bed of pre-sieved sand samples of diameter D mean = 180μm. The sand bed at low transport rate allows fairly long experiments without sand depletion upwind. A small container at the end of working section collects the transported sand from which we obtain the average flux.
Sand feeding in the wind tunnel has two main purposes. First, it guarantees that no depletion of sand occurs at the upwind end of the bed. Sand is fed into the tunnel via 5 tubes that distribute it fairly evenly across the bed. The feeding rate is approximately the same as that by which sand is transported downstream the tunnel.
Second, feeding is also useful for studying saltation when the shear stress is above the dynamic threshold as it simulates an infinite upwind fetch. The input of particle through feeding at the upwind end of the tunnel triggers saltation by adding momentum in the system.
Saltation is recorded on a standard video camera (2 megapixels) outside the tunnel. An observation point, consisting of a video camera and a pitot tube for wind measurements, is located at the downwind end of the sand bed before the sand collector, as indicated in Fig. 1. This focuses on the central part of the tunnel where saltating particles are illuminated when they cross a vertical laser sheet. Color frames are extracted from the videos (see Fig. 1), converted to gray-scale [0,1], and then into binary black and white images using a pixel threshold of 0.06 13 . Some particles may be lost during this selection. Saltation activity is observed with the increase of white pixels in the time series. From the video records, we measure the number and the length of the bursts. The video frames do not give access to the individual dynamics of the grains. The saturated flux is not quantified by the number of white pixels. We identify the start and end of a burst using a calibrated threshold to discriminate the burst activity from the noise, which is when the number of white pixels from one frame to the next changes by more than 30. We measure a time delay between two sand bursts if sand transport stops for more than 3 seconds, which was arbitrary chosen short enough to include long bursts with short delays between bursts at high wind velocities.
Initially the duration of experiments varied from 2-3 minutes (short-runs), but this produced very poor statistics. Therefore two series of new experiments were made where data were taken over three consecutive periods of 20 minutes (long-runs). The tunnel settings and the environmental conditions are described in the Supplemental Material.
In the numerical simulations, a three dimensional wind channel of dimensions D 700 50 7 5 mean 3 ( × × . ) with periodic boundary conditions in the direction of the wind contains a poly disperse 3D quiescent packing composed by 12 layers of hard spheres with Gaussian distributed diameters with D mean = 200 μm, size dispersion σ D = 0.15D mean , and density ρ s = 2650 kg/m 3 that are subjected to a gravitational field in vertical direction (y-direction) and a logarithmic wind velocity profile u(y) imposed in horizontal direction (x-direction), In Eq. (1), y 0 = D mean /30 is the roughness of the bed, h 0 the bed height, κ = 0.4 the von Kármán constant, and u * the wind shear velocity. The air density is ρ a = 1.174 kg/m 3 . We express the wind velocity through the dimensionless Shields number, defined as Similarly, we defined the fluid threshold Shields number as the minimum wind shear velocity for sand transport. The sand samples in the wind tunnel and the numerical spheres are within the typical range of fine wind blown particles. We normalize the Shields number in our results with θ c as it strongly depends on the weather conditions 14 . The turbulent flow velocity u splits into the mean stream velocity u(y) and a stochastic part u s : For u s , several measurements [15][16][17][18] have shown a highly non-Gaussian behavior of the Lagrangian acceleration distribution of fluid particles. Reynolds 19 formulated a system of stochastic differential equations (SDE) for the logarithm of the dissipation rate, χ = log(ε/〈 ε〉 ), with 〈 ε〉 as the mean dissipation rate that reproduces the Gaussian distribution of the velocities u s and the highly non-Gaussian distribution for the acceleration a s of the particles in fully developed turbulence. Within reasonable bounds, the probability distribution of the flow perturbation has negligible effect on the results in steady state. Every particle is attached to a tracer that generates random perturbations. The dissipation rate 〈 ε〉 controls the width of the probability distribution of these perturbations and consequently the rate of particle lift. Wider distributions produce stronger turbulent accelerations and lift more particles. Details of the Discrete Element implementation, wind profile integration considering the momentum exchange between the wind and the grains, and the turbulence model are discussed in the Supplemental Material.
The sand flux in the direction of the wind is defined as, is the area of the bottom of the channel, v i x and m i are, respectively, the particle velocity in the x-direction and the mass of the particle i. The saturated flux is the average flux in the stationary state. The reinjection of particles at one side of the system domain after they cross the periodic boundaries might act effectively as a sand feed and increase the numerical predictions for the sand flux.

Results
Without the sand feeding, the experiments in the wind tunnel for θ c < 0.0072 produced no grain at the end of the wind tunnel. The sand transport shoots up very steeply at θ c ≈ 0.0072 ± 2 × 10 −5 . Figure 2a shows the transported sand (red circles) collected at the container at the downwind end during the long runs. Within the error bars, saltation activity is very rare and the sand flux is virtually zero for θ/θ c = 0.93 as shown in the flux series on Fig. 2b. The gradual increase of θ/θ c shows the build up of the continuous flux (Fig. 2d) going through the discontinuous state (Fig. 2c). As the wind velocity increases, bursts superpose and the average burst size seen by the camera increases until the sand flux becomes continuous. Therefore, the interval for a discontinuous saturated flux is The numerical results from simulations (black squares) are discussed at the end of the manuscript. Figure 3a shows the time series of the sand transport in the experiments at θ c . An example of a sand burst is highlighted in red. The video camera focuses on a narrow window of 75 mm length, which does not allow to observe multiple bursts. The burst size in yellow and the time delay between bursts in blue display scaling laws. Figure 3b shows that the burst size distribution follows a power law with exponent − 0.7. The time delay (or waiting time) between bursts has also a power law distribution with exponent − 1.2, as shown in Fig. 3c. Scaling laws with comparable exponents can be also found in other complex systems with bursting activity, such as neuronal avalanches in rat cortex 20 . At θ/θ c = 1.24, the short and few intervals between bursts require much longer runs for getting good statistics. The high frequencies of the distribution are bounded by the noise stemming from the video camera. The low frequencies are bounded by the finite size of the wind tunnel which limits the formation of large low frequency eddies in the turbulent structure. Consequently, the probability distributions are restricted to a few orders of magnitude. The bursting activity was not prominent enough to observe a power law for θ < θ c . Similar power laws fitted reasonably the time delays of the bursting activity for θ/θ c < 1.2 (see Supplemental Material). However, for θ/θ c > 1.2, it is more difficult to obtain a clear distribution for the reducing time delays.
Whether a burst is generated locally from strong gusts at the area nearby or from the gradual development of the saltation initiated far upstream and propagating downstream cannot be determined unless we measure the velocity at a series of upstream observation points. A sand burst observed in the videos can be the superposition of local bursts generated at different locations of the wind tunnel. However, two different bursts do not have enough space to move alongside in the 0.6 m wide wind tunnel.
The aerodynamic lift initiates the sand bursts. Figure 4 shows the wind velocity at the observation point for a run at θ c . In red we show the intervals where the wind velocity was larger than the average wind speed, i.e., the speed corresponding to θ c . The time delay/waiting time between wind fluctuations above the average velocity also follows a power law of exponent − 1.6, as shown in Fig. 4b. Similarities in the probability distributions suggest a correlation between the occurrence of the bursts and the fluctuations in the wind speed. A definitive conclusion requires a complete data synchronisation between the wind measurements and the video records, which was unreachable with the state of art equipment and the collected material. Figure 5 gives insight about the triggering mechanism. It shows, in (a), that bursts disappear, when the Shields number decreases from θ c to θ/θ c = 0.94. Alternatively, the bursting appears after the Shields number increases from θ/θ c = 0.94 to θ c , shown in (b). Figure 6 shows two different experiments, with (red squares) and without feeding (black circles). With feeding it becomes possible to sustain a discontinuous flux even below the fluid threshold θ c . The occurrence of sand bursts is strongly dependent on the wind fluctuations around the fluid threshold. The momentum input from the sand feeding decreases the dependence of the sand bursts on the wind fluctuations and thus transport is sustained as long as the wind shear velocity is above the dynamic threshold. Thus sand flux can occur for Shields numbers that did not display any transport without feeding. The feeding significantly shifts the fluid threshold Shields number to θ/θ c = 0.73. The area in red shows the hysteresis zone that ends at θ/θ c = 1.12 where the transported sand is the same for both experiments within the error bars. The contribution of the feeding to the sand transport decreases as θ/ θ c → 1.12. We note that the ratio between threshold Shields number with and without feeding is around 0.70 (or u D /u t = 0.84).
In the numerical simulations, the random turbulent perturbations mimic the wind pick of sand grains and reproduce the sand bursts. Figure 7a shows the bursting activity at the sand flux at θ/θ c = 1.02 and 〈 ε〉 = 0.02 kg/(m.s 3 ) at the numerical simulations. From a completely settled particle bed, winds with dissipation rate 0.01 kg/(m.s 3 ) < 〈 ε〉 < 0.06 kg/(m.s 3 ) detach particles from the bed surface starting saltation. If other particles are neither ejected nor lifted by the stochastic force meanwhile, saltation stops. Saltation restarts if any other saltating particle is detached from the particle bed by the stochastic turbulent forces. Increasing 〈 ε〉 , the lift rate increases and the time delay between sand burst decreases until sand flux becomes continuous. Figure 7b shows the decrease of the time delays between bursts with the

Conclusions
In summary, the turbulent gusts create discontinuous saltation through bursting. A discontinuous saturated flux occurs at Shields numbers between θ θ < / . We can not expect a perfect match between our findings in the wind tunnel and Nature. In fact, no wind tunnel is able to reproduce the low frequency of big eddies normally observed in the field. Therefore, the contribution of the bigger eddies can only be measured in field experiments. However, close to the onset of saltation bursts are still quite small and thus dominated typically by smaller turbulent eddies.
Additionally, the same wind shear velocity could have different rates of mass transport in the field and in the wind tunnel. The wind shear velocity is a mean quantity which does not account the large  Aeolian sand bursts evolve downstream into streamers, which are so far, poorly understood. Future studies in the field should focus on these interesting coherent 3D structures.