Turbulent characteristics of saltation and uncertainty of saltation model parameters

It is widely recognised that saltation is a turbulent process, similar to other transport processes in the atmospheric boundary layer. Due to a lack of high-frequency observations, the statistic behaviour of saltation is so far not well understood. In this study, we use the data from the Japan–Australia Dust Experiment (JADE) to investigate the turbulent characteristics of saltation by analysing the probability density function, energy spectrum and intermittency of saltation fluxes. Threshold friction velocity, u∗t , and saltation coefficient, c0, are two important parameters in saltation models often assumed to be deterministic. As saltation is turbulent in nature, we argue that it is more reasonable to consider them as parameters obeying certain probability distributions. We estimate these distributions using the JADE data. The factors contributing to the stochasticity of u∗t and c0 are examined.


Introduction
It is well recognised that saltation, the hopping motion of sand grains near the Earth's surface, is a turbulent process (Bagnold, 1941). However, early studies focused mainly on its mean behaviour. The most well known is the Owen (Owen, 1964) saltation model, which predicts that the vertically integrated saltation flux is proportional to u * cubed, where u * is friction velocity, defined as u * = √ τ/ρ with τ being surface shear stress (N m −2 ) and ρ being air density (kg m −3 ). A dedicated investigation on turbulent saltation was conducted by Butterfield (1991), which revealed the significant variability of saltation fluxes concealed in conventional time-averaged data. Stout and Zobeck (1997) introduced the idea of saltation intermittency and pointed out that even when the averaged u * is below the threshold friction velocity, u * t , saltation can still intermittently occur. The latter authors emphasized saltation intermittency caused by fluctuations of turbulent wind, but stochasticity of u * t can also play a role. Turbulent saltation has attracted much attention in more recent years (e.g. McKenna-Neuman et al., 2000;Davidson-Arnott and Bauer, 2009;Sherman et al., 2017) and large-eddy simulation models have been under development to model the process (e.g. Dupont et al., 2013). However, due to a lack of high-frequency field observations of saltation fluxes, the statistical behaviour of turbulent saltation is to date not well understood.
A related problem is how saltation can be parameterised in wind erosion models. For example, for dust modelling, it is important to quantify saltation, as saltation bombardment is the main mechanism for dust emission. In wind erosion models, u * t is a key parameter which depends on many factors including soil texture, moisture, salt concentration, crust and surface roughness. In models, it is often expressed as 7596 D. Liu et al.: Turbulent characteristics of saltation and uncertainty of saltation model parameters al. (1993) scheme and f θ the Fécan et al. (1999) scheme. The corrections f sl and f cr are so far not well known.
For homogeneous saltation, the saltation flux can be computed using the Kawamura (1964) scheme, multiplied here by the fraction of erodible surface area σ f , where d is the particle diameter in a sand particle size range and g is acceleration due to gravity. The saltation coefficient, c 0 , is usually estimated empirically from field and/or wind tunnel experiments. It falls between 1.8 and 3.1 according to Kawamura (1964), and is commonly set to 2.6 (White, 1979) in wind erosion models. The total (all particle sizes) saltation flux, Q, is a particle-size weighted average of Q(d) where d 1 and d 2 define the upper and lower limits of saltation particle size respectively, and p s (d) is the soil particle size distribution. Observations show, however, c 0 varies considerably from case to case (e.g. Gillette et al., 1997;Leys, 1998), and as the data presented later in this paper show, for a given location, it may vary from day to day and even during a wind erosion event.
While wind erosion modules built into numerical weather and global climate models (e.g. Shao et al., 2011;Kok et al., 2014;Klose et al., 2014) are in general more sophisticated than what is described above and include a dust emission scheme, the estimate of Q is essentially done using Eqs. (1) to (3) or similar. Thus, the estimates of u * t and specification of c 0 are critical to wind erosion and dust modelling.
In most wind erosion models, both u * t and c 0 are treated as being deterministic. As saltation is turbulent, it is more rational to treat u * t and c 0 as parameters that satisfy certain probability distributions. Saltation intermittency also implies that u * t and c 0 depend on the scale of averaging. Shao and Mikami (2005) noticed that u * t for 10 min averaged Q and 1 min averaged Q are quite different. Namikas et al. (2003) and Ellis et al. (2012) have also noticed that averaging intervals of surface shear stress are important to quantifying sediment transport because both shear stress and saltation flux are turbulent.
Between 23 February and 14 March 2006, Ishizuka et al. (2008 carried out the Japan-Australia Dust Experiment (JADE) in Australia. In JADE, both u * and Q, together with a range of atmospheric and soil surface quantities, were measured at relatively high sampling rates. The loamy sand soil surface at the JADE site was very mobile and thus the JADE data are representative of surfaces that are almost ideal for sand drifting. In this study, we analyse some aspects of the turbulent behaviour of saltation using the JADE measurements of saltation fluxes. In light of the analysis, we ask the question of what the most likely values of u * t and c 0 are and how representative they are. We also estimate the probability distribution of the two parameters. .0 E). The size of field is about 1 km in the E-W direction and about 4 km in the N-S direction. A range of atmospheric variables, land surface properties, soil particlesize distributions and size-resolved sand and dust fluxes were measured. During the study period, 12 wind erosion episodes were recorded. The data set is particularly valuable in that particle-size-resolved sand and dust fluxes (Shao et al., 2011) were measured. The details of the experiments and data sets can be found in Ishizuka et al. (2008Ishizuka et al. ( , 2014 and hence only a brief summary is given here. In JADE, three sand particle counters (SPCs) (Yamada et al., 2002) were used to measure saltation at the 0.05, 0.1 and 0.3 m levels with a sampling rate of 1 Hz. A SLD (super luminescent diode) light source is used to detect particles flying through the light beam. The frequency of the input signal is 1-30 kHz, implying that particles moving at a speed of less than 30 m s −1 can be detected. A SPC measures the saltation of particles in the range of 39-654 µm in 32 bins with mean diameters of 39, 54, 69 µm etc. with irregular increments ranging between 15 and 23 µm. At each measurement height, the saltation flux density (M L −2 T −1 ), q, is obtained as the sum of q j (saltation flux for size bin j ) for the 32 size bins, i.e.
The saltation flux, Q, is then estimated by integrating q over height, namely In computing Q, we assume q = q 0 exp(−az) with q 0 and a being fitting parameters from the measurements. Prior to the field experiment, the SPCs were calibrated in the laboratory and during JADE. They were checked in a mobile wind tunnel at the site and compared with other saltation samplers. But as q was measured only at three heights, the vertical resolution of q is relatively poor and inaccuracies in the Q estimates are unavoidable, which we are unable to fully quantify. However, the profiles of q are well behaved and thus the inaccuracies in the absolute values of the Q estimates are not expected to be so large as to affect the conclusions of this study. Q is computed using the SPC data at 1 s intervals. We denote its time series by Q 1 s . From Q 1 s , the 1 min averages, Q 1 min , and 30 min averages of saltation fluxes, Q 30 min , are derived (Ishizuka, 2018). All these quantities are also computed for individual particle size bins as Atmospheric variables, including wind speed, air temperature and humidity at various levels, as well as radiation, precipitation, soil temperature and soil moisture were measured using an automatic weather station (AWS). These quantities were sampled at 5 s intervals and their averages over 1 min intervals were recorded. Two anemometers were mounted at heights of 0.53 and 2.16 m on a mast for measuring wind speed. Also available are the Monin-Obukhov length and sensible heat fluxes. From the wind measurements, surface roughness length z 0 and friction velocity u * are derived, assuming a logarithmic profile (with stability correction) of the mean wind. The roughness length for the experimental site is estimated to be 0.48 mm. Friction velocity is computed with 1 min averaged wind data, denoted by u * 1 min , and 30 min averaged wind data, denoted by u * 30 min (Ishizuka, 2018). In atmospheric boundarylayer studies, there is no standard for the length of time over which one should average wind to correctly estimate u * , but it is common to average it over 10 to 30 min. But the time period depends on the purpose of the averaging. If u * is used as a scaling velocity for the atmospheric boundary layer, e.g. as measure of turbulence intensity, it is necessary to average over a sufficiently large time interval to obtain a constant u * . In this paper, u * is a surrogate of shear stress, the variation of which drives that of saltation. Therefore, short averaging times are preferred, subject to them being larger than the response time of aeolian flux to shear stress. Anderson and Haff (1988) and Butterfield (1991) suggested that this response time is of the order of 1 s.
Observations of surface soil properties, including soil temperature and soil moisture, were made at 1 min intervals. The surface at the JADE site was relatively uniform. A survey of ground cover over an area of 900 × 900 m 2 at the site was made on 11 March 2006. The area was divided into nine tiles and surveyed along one transect of 300 m long in each tile. Photographs were taken every 5 m by looking down vertically to a point on the ground. Surface cover was estimated to be ∼ 0.02 (see Appendix of Shao et al., 2011).
The wind erosion model, as detailed in Shao et al. (2011), is used for computing the saltation fluxes using the JADE atmospheric and surface soil measurements as input. The saltation model component is as described in Section 1, consisting of Eqs. (1)-(3). The fraction of erodible surface area, σ f , used in Eq. (2), is 1 minus the fraction of surface cover. The soil particle size distribution (psd), p s (d), required for Eq. (3), is based on soil samples collected at the JADE site and analysed in the laboratory. The analysis was done using a Microtrac (Microtrac MT3300EX, Nikkiso Co. Ltd.), a particle size analyser based on laser diffraction light-scattering technology. Water was used for sample dispersion. Depending on the methods (pretreatment and ultrasonic vibration) used, the soil texture can be classified as sandy loam (clay 0.3 %, silt 25 % and sand 74.7 %) or loamy sand (clay 11 %, silt 35 % and sand 54 %). The sandy loam psd is used in this study, which has a mode at ∼ 180 µm (see Shao et al., 2011, Fig. 5, Method A).
The default value of c 0 is set to 2.6, as widely cited in the literature (e.g. White, 1979) and the default value of u * t is computed using Eq. (1) with u * t (d) computed using the Shao and Lu (2000) scheme, f λ using the Raupach et al. (1993) scheme, f θ the Fécan et al. (1999) scheme, and f sl and f cr set to one. The frontal area index λ and soil moisture θ are both observed data from JADE.

Method for parameter estimation
Different choices of c 0 and u * t would lead to different model-simulated saltation fluxes which may or may not agree well with the measurements. By fitting the simulated saltation fluxes to the measurements, we determine the optimal estimates of c 0 and u * t and the probability density function (PDF) of these parameters. The method based on the Bayesian theory is used for this purpose.
Suppose X = ( x 1 , x 2 , . . ., x n ) is a measurement vector, with x i being the measured value at time t i , and A is a model with a forcing vector F and model parameter vector β. Let the initial state of the system be i 0 , so the modelled value of the system, X = (x 1 , x 2 , . . ., x n ), can be expressed as The error vector is given by E(β) = X − X, here and is fully attributed to β. Given X, the posterior parameter PDF, p(β X) , can be estimated from the Bayes theorem: where p(β) is the prior parameter PDF and p( X |β) is the likelihood. If p(β) is given, then the problem of finding p(β X) reduces to finding the maximum likelihood. Assuming the error residuals are independent and Gaussian distributed with constant variance, σ 2 , the likelihood can be written as In this case, maximising the likelihood is equivalent to minimising the error, i.e. The solution of Eq. (10) gives an optimal (i.e. with maximum likelihood) estimate of mean β. This is the popular least squares method. A disadvantage of the method is that it assumes a Gaussian posterior parameter PDF and computing the β variance requires prior knowledge of the accuracy of the data.
As an alternative, the approximate Bayesian computation (ABC) method has been proposed (e.g. Vrugt and Sadegh, 2013). It is argued that a parameter value β * should be a sample of p(β X) as long as the distance between the observed and simulated data is less than a small positive value This procedure explicitly provides an estimate of parameter PDF for a given data set. The ABC method is numerically simple: from a prior PDF (e.g. uniform) of β a β * is stochastically generated and the model is run. If Eq. (10) is satisfied, then β * is accepted or otherwise rejected. This procedure is repeated and the a priori PDF of β is mapped to a posterior PDF of β. The ABC method has the disadvantage that it is numerically inefficient. More efficient techniques based on the same principle exist, e.g. the Markov chain Monte Carlo simulation (Sadegh and Vrugt, 2014). In this study, we apply the Differential Evolution Adaptive Metropolis (DREAM) algorithm proposed by Vrugt et al. (2011) for estimation of hydrologic model parameters. The algorithm integrates differential evolution (Storn and Price, 1997) and self-adaptive randomised subspace sampling to accelerate a Markov chain Monte Carlo simulation. A full description of the DREAM algorithm is beyond the scope of our study. Interested readers should refer to the above-cited references for details.
3 Statistical features of saltation

Time series
To provide an overview of the data set used in this study, Fig. 1 shows the time series of Q 1 min and u * 1 min , and Fig. 2, Q 30 min and u * 30 min . During the 20-day period, aeolian sand drift occurred almost every day at the site according to the field logging book, but only 12 events were recorded using the SPCs. Saltation fluxes were not measured on days 55, 58, 59, 64 and then days 66 to 70 due to either instrument maintenance or use of the SPCs for other purposes (e.g. wind tunnel experiments). The figures show that both Q and u * fluctuate significantly and saltation is turbulent. Figure 1b shows an enlarged plot of the Q 1 min and u * 1 min time series for days 61 and 62. At the JADE site, u * t was about 0.2 m s −1 . On day 61, u * was mostly larger than this value and saltation was almost continuous, while on day 62, u * was close to this value and weak saltation occurred frequently also when u * was below 0.2 m s −1 . Figure 2b is as Fig. 1b   time series of Q 30 min and u * 30 min can be directly seen in Fig. 2b. In Fig. 3a, b and c, Q is plotted against u 3 * . Several interesting features can be identified. For the majority of the points, the Q ∼ u 3 * relationship appears to hold, but this relationship can vary significantly even for the same data set from event to event. For example, large differences exist between days 70 and 71 (denoted by D70-71, an event of intensive wind erosion) and day 72 (a day of weak wind erosion), as seen in both Fig. 3a and b. There may be many likely reasons for the differences in the Q ∼ u * relationship but the most conspicuous are differences in atmospheric turbulence (e.g. gustiness) and time-varying surface conditions (e.g. particle sorting and aerodynamic roughness). Figure 3d shows the time series of (u * 1 min − u * 30 min ), a measure of turbulent fluctuations. It is seen that saltation is associated with not only high surface shear stress but also high shear stress fluctuations. The large difference in the Q ∼ u * relationship between D70-71 and D72 (Fig. 3b) is probably attributed to the strong differences in turbulent fluctuations (Fig. 3d): D70-71 was a hot gusty day with top (2 cm) soil temperature reaching 53 • C, while D72 was cooler and less gusty with soil temperature 5 • C lower. Also hysteresis is observed in the Q ∼ u * relationship, as shown in Fig. 3c, using D71 and D72 as examples. Figure 3d shows that for all three events selected (D70-71, D71 and D72), saltation has a relatively short (0.5 to 2 h) strengthening phase, followed by a longer weakening phase. During an erosion event for the same u * , saltation is stronger in the strengthening phase than in the weakening phase. An exam-ination of Fig. 3d suggests that the hysteresis cannot simply be attributed to the intensity of turbulence. We speculate that it is probably more related to flow-saltation feedbacks (e.g. stronger splash entrainment in the strengthening phase) and the modification of surface aerodynamic conditions (e.g. particle sorting and reduced surface roughness Reynolds number).

Probability density function of saltation fluxes
How well the saltation model performs, whether u * t and c 0 are universal and how they are probabilistically distributed must depend on the turbulent properties of saltation. As the JADE saltation fluxes are sampled at 1 Hz, we can use these data to examine (to some degree) the statistical behaviour of saltation. In Fig. 4, the PDFs of the saltation fluxes for different particle size groups are plotted and computed using Q 1 s and Q 1 min . It is seen that the PDFs generally behave as In the case of Q 1 s , there seems to be a distinct change in α at a critical value of Q c ∼ 3 g m −1 s −1 , with α ∼ 1 for Q < Q c and α ∼ 4 for Q > Q c . The PDFs derived from Q 1 min appear to follow the basic functional form of Eq. (11). Again, α is about 1 and tends to be larger for large Q values. Figure 4 shows that the PDFs of Q depend significantly on the interval of time averaging; i.e. after averaging, smaller saltation fluxes become more frequent. This is because the time series of Q 1 s is more intermittent (see also Fig. 6).  The PDFs of Q 1 s and Q 1 min integrated over all particles are shown in Fig. 5b. Again, the PDFs show the general behaviour of p(Q) ∼ Q −1 . In theory, p(Q) can be derived from the PDF of u * , p(u * ). From Eq. (2), we have This can be used to obtain (13) Figure 5a shows the p(u * ) estimated from u * 1 min together with the fitted Weibull distribution. For the fitting, we ensure that p(u * ) for u * > 0.2 m s −1 is best approximated. Figure 5b shows the p(Q) estimated from Q 1 min . We computed p(Q) using Eq. (13) with the fitted p(u * ), assuming u * t = 0.2 m s −1 and c 0 = 2.6. It is seen that the observed and modelled p(Q) have qualitative similarities but using Eqs. (12) and (13) we cannot reproduce the observed p(Q) well. For example, the model fails to predict the low frequent strong saltation fluxes and the mode of saltation fluxes. Tests using several smaller u * t values (0, 0.05 and 0.1) are also carried out. With smaller u * t values, the mode of the predicted saltation fluxes is shifted to smaller values, but the predictions are far from satisfactory.

Saltation intermittency
Following Stout and Zobeck (1997), the intermittency of saltation, γ , is defined as the fraction of time during which saltation occurs at a given point in a given time period.
It should be pointed out that saltation intermittency describes only the behaviour of the process at u * ∼ u * t ; i.e. saltation intermittency is merely a special case of turbulent saltation. Several formulations of γ are possible. Stout and Zobeck (1997) assumed that saltation occurs only in time windows in which u * exceeds u * t . Therefore, if p(u * ) is known, then γ for a given u * t can be estimated as Stout and Zobeck (1997) used the counts per second of sand impacts on a piezoelectric crystal saltation sensor as a measure of saltation activity and found that γ a rarely exceeded 0.5. In Eq. (14a) u * t is fixed and thus saltation intermittency is attributed entirely to the fluctuations of u * . In reality, u * t also fluctuates and satisfies certain PDFs (Raffaele et al., 2016).
In analogy to Eq. (14a), γ for a given u * can be estimated as More generally, we can define saltation intermittency as or Equations (14c) and (14d) The computation of saltation intermittency function γ a (u * t ) is done by integrating p(u * ) (Fig. 5a) to fixed value of u * t . In Fig. 6a, γ a as a function of u * t is plotted. The behaviour of γ a (u * t ) is as expected: it is one at u * t = 0 and decreases to zero at about u * t = 0.5 m s −1 as in the case of JADE, u * rarely exceeded this value. For u * t = 0.2 m s −1 , γ a is 0.35, comparable with the result of Stout and Zobeck (1997) who reported an intermittency of 0.4. As p(u * t ) is not known, Eq. (14b) cannot be used directly, but we can compute γ b (u * ) using the JADE data. First, it is computed using Q 1 min . This is done by selecting a fixed u * say u * c , and counting the time fraction, T u * , which satisfies |u * − u * c | < ε (used is ε = 0.05 m s −1 ) and the time fraction, T Q1 min , which satisfies |u * − u * c | < ε and Q 1 min > 0. By definition, saltation intermittency is T Q1 min /T u * as plotted in Fig. 6a. It is seen that for Q 1 min , γ b (u * ) increases from about 0.6 at u * ∼ 0.1 m s −1 to about one at u * = 0.3 m s −1 . This shows that in JADE a considerable fraction of the saltation fluxes was recorded at u * below the perceived threshold friction velocity (about 0.2 m s −1 ), saltation is more intermit under weak wind conditions and becomes non-intermittent for u * > 0.3 m s −1 . The increase in γ b (u * ) with decreasing u * for u * < 0.1 m s −1 is, however, unexpected. The expected γ b (u * ) for small u * is as depicted using the dashed line. A likely reason for the unexpected behaviour of γ b (u * )is that during a wind erosion event, grains in saltation may continue to hop even when u * is temporarily reduced to small values. The uncertainty in the data also needs to be considered, as the sample size for determining the ratio T Q1 min /T u * becomes smaller. More complete data sets are required to answer these questions. Finally, γ c is computed by using Eq. (14d) and is found to be around 0.73. For the 1 s case, we cannot plot γ b as a function of u * , because u * is not available at such high frequency. We computed γ c for individual particle size groups Figure 6. (a) Saltation intermittency function γ a (u * t ) and γ b (u * ). See text for more details. (b) γ c as a function of particle size for Q 1 s , Q 1 min and Q 30 min . (Fig. 6b) using Q 1 s , Q 1 min and Q 30 min , which is the time fraction of saltation for a given particle size, d, during the saltation event. It is found that γ c (d) decreases with d; i.e. the saltation of larger particles is more intermittent. Also, γ c (d) increases with increased averaging time intervals, implying that the small-scale features of turbulence play an important role in intermittent saltation.

Spectrum of saltation fluxes
Spectral analysis is widely used for characterising the variations in a stochastic process on different scales. Using the JADE data, we computed the power spectrum of saltation fluxes, P Q (f ) at frequency f , and of friction velocity, P u * (f ), using a non-uniform discrete Fourier transform. For comparison, the power spectra are normalised with the respective variances of the signal. In atmospheric boundarylayer studies, the spectra of various turbulence quantities have been thoroughly investigated (Stull, 1988). Examples of spectra from Reynolds shear stress can be found in Mc-Naughton and Laubach (2000). Figure 7 shows P Q (f ) and P u * (f ) (Fig. 7a) as well their co-spectrum (Fig. 7b). P Q (f ) is computed using both Q 1 s and Q 1 min , and P u * (f ) with u * 1 min . It is seen that the power spectra of Q and u * have qualitatively very similar behaviour. Both have a maximum at about 10 −5 Hz, a minimum at about 10 −4 Hz and another peak at about 2 × 10 −3 Hz. The maximum at 10 −5 Hz is related to the diurnal patterns and changing synoptic events, which drive the wind erosion episodes, the minimum at 10 −4 Hz is due to the lack of turbulent winds at the timescale of several hours, while the peak at 2 × 10 −3 Hz is caused by the minute-scale gusty winds/large eddies in turbulent flows. Also the Q − u * co-spectrum shows that Q and u * are most strongly correlated on diurnal/synoptic and gust/large-eddy timescales. P Q (f ) computed using Q 1 s again reveals the peaks at 10 −5 Hz and at 2 × 10 −3 Hz. The power of the Q spectrum then decreases with frequency. As the sampling rate of saltation flux is limited to 1 s in this study, the features of P Q (f ) at frequencies larger than 0.5 Hz are not resolved.

Estimates of saltation model parameters
Given the turbulent nature of saltation, it is rational to treat u * t and c 0 in the saltation model as parameters obeying certain probability distributions. To examine the behaviour of these parameters, we introduce two coefficients r c 0 and r u * t and multiply them by the theoretical values of c 0 and u * t , respectively, in Eq. (2), i.e. u * t = r u * t u * t,theory c 0 = r c 0 c 0,theory .
As introduced in Sect. 1, we assumed c 0,theory = 2.6 and computed u * t,theory using Eq. (1) with observed soil moisture and fraction of cover. The two coefficients r c 0 and r u * t are varied to generate a model estimate of Q using Eqs. (2) and (3) with observed u * . The probability distributions of r c 0 and r u * t are estimated using the following techniques. Let us denote the time series of the modelled saltation flux by Q M,i, (i = 1, N ) and the corresponding measurement by Q D,i .The absolute error, δQ A , and Nash coefficient, I Nash , are used as measures for the goodness of the agreement between the model and the measurement. They are defined as  The prior PDFs of r c 0 and r u * t are assumed to be uniform.
In the numerical experiment, we randomly generate r c 0 and r u * t and seek their values, such that δQ A ≤ ε and I Nash > η. These experiments are repeated for Q 1 min and Q 30 min . The plots of δQ A and I Nash as functions of r c 0 and r u * t show that for certain values of r c 0 and r u * t , the above conditions are satisfied. Figure 8 shows that for Q 1 min , the best simulation is achieved with r c 0 = 1.23 and r u * t = 1.05, while for the Q 30 min , with r c 0 = 0.94 and r u * t = 0.91. This suggests that the optimal estimates of u * t and c 0 are close to the corresponding theoretic values but are dependent on the timeaveraging intervals, with both u * t and c 0 being larger for shorter averaging intervals. The parameter PDFs p(r u * t ) and p(r c 0 ) are estimated with the DREAM algorithm, again using the absolute error and the Nash coefficient as goodness of agreement between the model simulated and measured saltation fluxes. The results are shown in Fig. 9. All PDFs are fitted to a distribution. As seen in Fig. 9a and c, the most frequent r u * t values are 1.12 and 1.04 for Q 1 min and Q 30 min , close to the estimates of 1.05 and 0.91 found in Fig. 8. For Q 1 min , r u * t is ∼ 1.12 ± 0.2 and for Q 30 min ∼ 1.04 ± 0.3. This implies that saltation sometimes occurs when u * is below the theoretical u * t value and sometimes saltation does not occur even when u * is above it, as already seen in Fig. 6a. In the case of p(r c 0 ) ( Fig. 9c and d), the most frequent values of r c 0 for Q 1 min and Q 30 min are respectively 1.04 and 0.92, close to the optimal estimates of 1.23 and 0.94 shown in Fig. 8. But r c 0 varies over a wide range, for instance, for Q 30 min between 0.5 and 5; i.e. c 0 is a rather stochastic parameter.
In nature, many factors influence sediment transport, but the stochasticity of the parameters is determined primarily by the turbulent fluctuations of friction velocity (or surface shear stress), the randomness of threshold friction velocity and soil particle size distribution (representing particle response to forcing). Studies have shown, for instance, that small changes in soil moisture can have a large influence on saltation (Ishizuka et al., 2008) and soil moisture in the very top soil layer can vary significantly over relatively short time periods. Over the period of 18 days during JADE, soil moisture in the top 0.05 m layer varied between 0.02 and 0.04 m 3 m −3 (4 and 8 % in relative soil moisture, assuming a saturation soil moisture of 0.5 m 3 m −3 ). In this study, the influence of soil moisture on saltation is accounted for via Eq. (1) using the soil moisture measurements in the top 0.05 m layer (see also Fig. 4a in Shao et al., 2011). While measured soil moisture is used in the wind erosion model, the randomness associated with its spatial-temporal variations is not, which is most likely reflected in the stochasticity of u * t .
The stochasticity of c 0 arises because saltation fluctuates depending on turbulence and particle size. To demonstrate this, we divided the time series of the saltation fluxes into two subsets, one with Q D,i ≤ 3 g m −1 s −1 representing weak saltation and one with Q D,i > 3 g m −1 s −1 representing significant saltation. This separation is arbitrary but sufficient for making the point that c 0 depends on u * , which is also a measure of turbulence intensity. The parameter PDFs, p(r u * t ) and p(r c 0 ), for the subset Q D,i ≤ 3 g m −1 s −1 is shown in Fig. 10. For Q 1 min and Q 30 min , the most frequent r u * t values are now 0.99 and 0.85, somewhat smaller than the estimated values for the full set (Fig. 9). In comparison, the most frequent r c 0 values are now 0.30 and 0.29, three to four times smaller than for the case when the full set is considered (Fig. 9). This suggests that c 0 has a clear dependency on u * and is smaller for smaller u * . This is because saltation is more intermittent in the case of smaller u * (i.e. smaller excess shear stress) and thus, c 0 , a descriptor of the relation between time-averaged saltation flux and friction velocity, is smaller for more intermittent saltation.
We fit the PDFs, p(r u * t ) and p(r c 0 ), for individual particle size bins and found that the most frequent r u * t values do not differ substantially among the particle sizes, but r c 0 depends systematically on particle size. For example, the most frequent r c 0 values for 101, 151, 203, 315 and 398 µm are 0.5, 1.3, 1.7, 3.1 and 4.0. These values are obtained by first estimating p(r c 0 ) for the individual particle size bins with the measured saltation flux for the corresponding bins and then normalising p(r c 0 ) with the mass fraction of the size bins of the parent soil. A least squares curve fitting shows that the most frequent r c 0 value depends almost perfectly (R 2 = 0.996) linearly on particle size: for the particle size range (100 to 400 µm) we tested, with d being particle size in µm. We have shown that both u * t and c 0 satisfy certain PDFs that depend on the properties of the surface, atmospheric turbulence and soil particle size. Figure 9 shows that for a fixed choice of u * t and c 0 , even if they are optimally chosen, a portion of the measurements cannot be represented by the model. Then, how does the saltation model perform if a single fixed u * t and a single fixed c 0 are used, as is often the case in aeolian models? The p(Q) computed using the model and derived from the JADE measurements are shown for Q 1 min and Q 30 min in Fig. 11. The model is applied to estimate the saltation flux for individual particle size groups using the optimally estimated u * t and c 0 (with r u * t = 1.12 and r c 0 = 1.04 for Q 1 min , and r u * t = 1.04 and r c 0 = 0.92 for Q 30 min ), and the total saltation flux is computed by integration over all particle size groups, i.e. using Eq. (3). Figure 11 shows that for this option, the model overpredicts the probability of large Q but underpredicts the probability of small Q in both cases of Q 1 min and Q 30 min . Obviously, to better reproduce the Q 1 min and Q 30 min PDFs, more values of r u * t and r c 0 sampled from the parameter PDFs are required. We have therefore modelled Q 1 min with other choices of r u * t (1.12 and 0.56) and r c 0 (2.08, 0.01) and plotted the corresponding Q 1 min PDFs as well as the averaged Q 1 min PDF of the three simulations. Similarly, we performed Q 30 min model simulations with other r u * t (1.04) and r c 0 (1.84) values and examined the Q 30 min PDFs. With the additional choices of the r u * t and r c 0 values, the Q 1 min and Q 30 min PDFs can be better reproduced.

Summary
In this paper, we used the JADE data on saltation fluxes (resolution 1 s) and frictional velocity (resolution 1 min) to analyse the statistical behaviour of turbulent saltation and estimate the probability distribution of two important parameters in a saltation model, namely the threshold friction velocity, u * t , and saltation coefficient, c 0 .
Saltation fluxes show rich variations on different scales. It is found that while the widely used Q ∼ u 3 * relationship holds in general, it can vary significantly between different wind erosion events. In several wind erosion events observed in JADE, saltation hysteresis occurred. We examined the probability density function of the saltation fluxes, p(Q), and found that it generally behaves like Q −α with α ∼ 1. For Q 1 s , there is a distinct change in α at Q = 3 ∼ 4 g m −1 s −1 with α ∼ 1 for smaller Q and α ∼ 4.0 larger Q. It is shown that p(Q) is dependent on the averaging time intervals as a consequence of saltation intermittency.
We introduced the saltation intermittency functions γ a (u * t ), γ b (u * ) and redefined saltation intermittency γ c as the fraction of time during which saltation occurs at a given point in a given time period, and computed these saltation