Modelling aeolian sand transport using a dynamic mass balancing approach

Knowledgeofthechangingrateofsediment ﬂ uxinspaceandtimeisessentialforquantifyingsurfaceerosionand deposition in desert landscapes. Whilst many aeolian studies have relied on time-averaged parameters such as wind velocity ( U ) and wind shear velocity ( u ⁎ ) to determine sediment ﬂ ux, there is increasing ﬁ eld evidence that high-frequency turbulence is an important driving force behind the entrainment and transport of sand. At thisscale ofanalysis, inertiain thesaltation system causeschanges in sediment transport to lag behindde/accel-erationsin ﬂ ow.However,saltationinertiahasyettobeincorporatedintoafunctionalsandtransportmodelthat can be used for predictive purposes. In this study, we present a new transport model that dynamically balances the sand mass being transported in the wind ﬂ ow. The ‘ dynamic mass balance ’ (DMB) model we present accounts for high-frequency variations in the horizontal ( u ) component of wind ﬂ ow, as saltation is most strongly associatedwiththepositive u componentofthewind.TheperformanceoftheDMBmodelistestedby ﬁ ttingitto two ﬁ eld-derived (Namibia's Skeleton Coast) datasets of wind velocity and sediment transport: (i) a 10-min (10 Hz measurement resolution) dataset; (ii) a 2-h (1 Hz measurement resolution) dataset. The DMB model is shown to outperform two existing models that rely on time-averaged wind velocity data (e.g. Radok, 1977; Dong etal.,2003),when predicting sandtransport over the two experiments.Forall measurement averaging intervals presented in this study (10 Hz – 10 min), the DMB model predicted total saltation count to within at least 0.48%, whereas the Radok and Dong models over- or underestimated total count by up to 5.50% and 20.53% re-spectively. The DMB model also produced more realistic (less ‘ peaky ’ ) time series of sand ﬂ ux than the other two models, and a more accurate distribution of sand ﬂ ux data. The best predictions of total sand transport are achieved using our DMB model at a temporal resolution of 4 s in cases where the temporal scale of investigation is relatively short (on the order of minutes), and at a resolution of 1 min for longer wind and transport datasets (ontheorderofhours).Theproposednewsandtransportmodelcouldprovetobesigni ﬁ cantforintegratingtur-bulence-scale transport processes into longer-term, macro-scale landscape modelling of drylands. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
A central aspect of understanding desert landscape dynamics relies on the quantification of surface deposition and erosion. Numerous formulae have been developed to predict sediment flux (q) and these generally rely on temporally-and spatially-averaged measures of erosivity such as wind velocity (U) and shear velocity (u ⁎ ) (e.g. Bagnold, 1941;Kawamura, 1951;Zingg, 1953;Owen, 1964;Radok, 1977;Lettau and Lettau, 1978;Dong et al., 2003). No single sediment flux model has been found to be broadly applicable under the diversity of natural environmental conditions that exist, as most are derived theoretically or from wind-tunnel observations without rigorous field-testing (Leenders et al., 2011). Large discrepancies arise between field-derived and modelled rates of sand flux (Sherman et al., 2013) because existing models employ mean parameter values to determine the aggregate behaviour of a landscape, when aeolian transport often appears to be dominated by extremes in erosivity and erodibility parameters (Walker and Nickling, 2002;Okin and Gillette, 2004;Wiggs et al., 2004;King et al., 2005;Okin, 2005;Raupach et al., 2006;Walter et al., 2012;Li et al., 2013).
In particular, the use of u ⁎ for estimating sand transport has been contested due to its dependence on a log-linear relationship between wind velocity and height, which collapses at high measurement frequencies (Bauer et al., 1992;Namikas et al., 2003;Wiggs et al., 2004). Some wind tunnel and field experiments have identified relationships between sediment transport and u ⁎ at timescales on the order of 10 0 -10 1 s (Butterfield, 1991(Butterfield, , 1993Spies et al., 2000;Namikas et al., 2003). Whilst this is incompatible with the traditional u ⁎ approachwhich requires wind shear to be calculated over periods of minutes to hourssuppressing the short-term variability of the transport system by using time-averaged series can provide reasonable relationships between wind velocity and sediment flux (Wiggs et al., 2004;Baas and Sherman, 2005b;Davidson-Arnott et al., 2012).
There is increasing evidence that high frequency turbulence is an important driving force behind aeolian sediment entrainment and transport in the boundary layer (Butterfield, 1991;Sterk et al., 1998;Namikas et al., 2003;Schönfeldt and von Löwis, 2003;Baas and Sherman, 2005a;Leenders et al., 2005;Baddock et al., 2011;Weaver and Wiggs, 2011;Wiggs and Weaver, 2012;Chapman et al., 2013). This has followed the relatively recent field application of sonic anemometers and high-frequency impact sensors, which allow quasi-instantaneous fluctuations from the mean wind velocity in the horizontal (u′), lateral (v′) and vertical (w′) directions, and therefore the Reynolds shear stress (−ρu0w0, where ρ is air density), to be measured in the field. Findings from physical experiments (e.g. Kline et al., 1967;Grass, 1971;Lu and Willmarth, 1973;Roy et al., 1999) and numerical simulations (e.g. Kim et al., 1987) suggest that certain characteristic fluid excursions are of greater importance to general turbulent flow dynamics than are spectra of fluctuating motions or eddies. Quadrant analysis (Lu and Willmarth, 1973) defines four discrete categories of momentum exchange in turbulent airflow, based on the relative signs of u′ and w′: ejections (Q2 events; u′ b 0, w′ N 0), sweeps (Q4; u′N 0, w′ b 0); outward interactions (Q1; u′N0, w′ N 0) and inward interactions (Q3; u′ b 0, − w′ b 0). Research in the fluvial domain suggests that ejection and sweep events in combination (the 'bursting process') are a key determinant of sediment transport (e.g. Best and Kostaschuk, 2002). However, several aeolian studies (e.g. Sterk et al., 1998;Schönfeldt and von Löwis, 2003;van Boxel et al., 2004;Leenders et al., 2005;Baddock et al., 2011;Wiggs and Weaver, 2012;Chapman et al., 2013) have shown that wind-driven saltation flux is better associated with turbulent structures characterised by a positive streamwise fluctuating velocity (u′N 0, sweeps and outward interactions). This difference between turbulent processes operating in the aeolian and fluvial domains is likely to be due to the greater disparity in density ratio between sand particles and wind compared with water (Wiggs and Weaver, 2012).
Despite improvements in high-frequency anemometry and saltation impact sensors that allow for many of the relevant transport processes to be measured in the field, the relative importance of turbulence for aeolian sediment transport remains unclear. This is partly due to the spatially discrete nature of transport events such as sand streamers (Butterfield, 1991;Baas and Sherman, 2005a), and inertia in the saltation system causing changes in sediment transport to lag behind de/accelerations in flow (Lee, 1987;Anderson and Haff, 1988;Butterfield, 1991Butterfield, , 1993Butterfield, , 1999Spies et al., 2000;Wiggs and Weaver, 2012). This lag has been notably explored through the concept of the saturation length L sat , a single length characterising the distance over which sediment flux adapts to changes in the wind shear (Sauermann et al., 2001;Andreotti et al., 2010). The saturation time T sat can be analytically derived from this, and is proportional to the reciprocal of the sum of localised erosion and deposition rates in a system (Gao et al., 2014). Evidently this makes it difficult for T sat to be accurately quantified in field settings over relatively short timescales. Whilst L sat has been shown to control the length of the smallest dunes evolving from flat sand beds (Claudin and Andreotti, 2006;Claudin et al., 2013) and to be a key parameter for the realistic modelling of barchan dunes (Narteau et al., 2009), it remains difficult to isolate the mechanism responsible for limiting saturation in aeolian systems (see Andreotti et al., 2010;Pähtz et al., 2015). As a result of these complexities, no empirical, predictive model has explicitly incorporated the dynamic balance between the mass of sand moving into (i.e. erosion) and dropping out of (i.e. deposition) the wind flow.
The aim of this study is to derive and test a new sediment transport model that accounts for the lag in flux that emerges when wind flow is analysed at sufficiently high frequencies. The model is fitted to empirical field data, and its performance at a range of measurement resolutions and averaging intervals (10 Hz-10 min) is compared to two existing transport models (Radok, 1977;Dong et al., 2003), using two separate wind and sand flux datasets.

New dynamic mass balance (DMB) model for aeolian sand transport
A new aeolian sediment transport model is proposed that incorporates the effects of a lag between high frequency flow turbulence and sediment flux. In the generalised model, referred to in this study as the 'dynamic mass balance (DMB) model', an increase in mass flux depends on above-threshold horizontal wind velocity. Whilst research in the fluvial domain (e.g. Best and Kostaschuk, 2002) has emphasized the importance of w for sediment transport, evidence from aeolian literature (e.g. Sterk et al., 1998;Schönfeldt and von Löwis, 2003;van Boxel et al., 2004;Leenders et al., 2005;Baddock et al., 2011;Wiggs and Weaver, 2012;Chapman et al., 2013;Mayaud et al., 2016) suggests that the w component tends to be an order of magnitude lower than u in natural wind flow, and that wind-driven flux is more simply associated with a positive u component. This is supported by the findings presented in Section 4.1.1, which show that the w component is poorly correlated with flux data. Therefore, vertical wind velocity is not incorporated in our DMB model.
The settling rate of saltating grains is explicitly accounted for in the DMB model. The model is formulated as a differential (i.e. time-rateof-change) equation that incorporates inputs to and outputs from a notional body of transported sand, such that the transported mass asymptotically approaches equilibrium under constant conditions, meaning a finite time is needed for both entraining and dropping the transported material. This behaviour reflects the inherent inertia in the saltation system that causes sediment transport to lag behind flow de/acceleration (Lee, 1987;Anderson and Haff, 1988;Butterfield, 1991Butterfield, , 1993Butterfield, , 1999Spies et al., 2000;Wiggs and Weaver, 2012), and is conceptually equivalent to the T sat parameter (Andreotti et al., 2010;Gao et al., 2014;Pähtz et al., 2015). The DMB model is given as: where Q iP is the predicted mass flux for the given time interval (kg m −1 s − 1 ), u is the mean horizontal wind velocity over the given time interval (m s − 1 ), u t is the horizontal wind velocity threshold (m s −1 ), d is the grain diameter (mm), D is a reference grain diameter (mm), and a 1 (kg m −3 s), b 1 (dimensionless) and c 1 (dimensionless) are empirically fitted constants. In this way, the first square bracket in the model describes instantaneous sediment entrainment. This is switched 'on' or 'off' depending on whether the wind velocity is above threshold, using the threshold criterion H. Debate exists within the literature about the power to which U or u ⁎ should be raised (Sarre, 1988;Ho et al., 2011;Sherman and Li, 2011). Here, sand flux is proportional to the square of above-threshold wind velocity, as is the case for many existing transport models. The second square bracket in the model describes depositional processes. Parameter b 1 controls the vertical settling of grains due to gravity (represented by the ratio d/D), a fundamental attribute of sediment transport that is primarily controlled by the grain size (Farrell and Sherman, 2015). Parameter c 1 controls the vertical settling of grains due to horizontal wind velocity reduction (represented by the ratio u/u t ), which reflects the fact that the saturation length L sat is strongly controlled by the relation between the gravitational constant and the average sand particle speed above the surface (Pähtz et al., 2015).

models
Existing sand transport models rely on time-averaged parameters such as mean wind velocity (U) and wind shear velocity (u ⁎ ) to determine sediment flux, whereas our DMB model, by virtue of its differential formulation, is well-suited to being employed with higher-frequency wind velocity and sand transport data. Two commonly used sand transport models were chosen to assess the comparative performance of our new model. The Radok (1977) and Dong et al. (2003) models formulate sediment flux as a function of U, and are therefore theoretically similar to other models employing u ⁎ due to their usage of mean wind characteristics. The Radok and Dong models have been shown to predict sediment transport more accurately than other U-based transport equations (Leenders et al., 2011), so they are considered good models against which to test the performance of our DMB model. The Radok model was initially developed to describe the movement of snow. It is fully empirical but has the limitation of not incorporating a threshold term, and therefore predicts transport even when wind velocity is below threshold. Therefore, it is theoretically only applicable in cases where wind velocity exceeds the threshold (Leenders et al., 2011). The Radok model is given by: where a 2 and b 2 are empirically fitted constants (kg m −1 s −1 and s m −1 respectively).
The semi-empirical model of Dong et al. (2003) has been widely applied at a range of temporal and spatial scales (e.g. Leenders et al., 2011;Bailey and Thomas, 2014). It incorporates a threshold wind velocity that should allow the sediment transport formulation to be applied to complete wind datasets. The Dong model is given by: where ρ is the density of air (1.25 kg m −3 ), g is acceleration due to gravity (9.81 m s −2 ), and a 3 is an empirically fitted constant. It was noted in preliminary tests that despite the inclusion of a threshold velocity term, the Dong model predicted some mass flux when wind velocities were slightly below threshold, but no mass flux at the actual threshold velocity. This arises from the fact that for Dong's model, ð1− ut u Þ 2 →0 as u → u t .
Therefore, in this study a threshold condition was imposed on the Dong formulation, such that no mass flux was predicted below the threshold wind velocity.

Exploring model parameters
By inserting a variety of parameter values into the three model formulations, it is possible to investigate the sensitivity of model behaviour to each parameter. Fig. 1 displays the effect of varying different parameters on the shapes of sediment flux time series predicted by the three models. The models were initially run with estimates based on extensive preliminary testing, and with values ±30% of this initial estimate. The theoretical threshold wind velocity was set at 6.4 m s −1 , to reflect the threshold velocity derived from the dataset of Experiment 1 (see Section 4.1).
As shown in Fig. 1, the Radok and Dong models respond much more rapidly (i.e. with no lag) to sudden increases and decreases in wind velocity compared to the DMB model. The Radok (Fig. 1e,f) and Dong (Fig. 1g) parameters effectively alter the magnitude of simulated transport but do not alter the onset or rate of decline in saltation activity. In contrast, the DMB model clearly displays a lag in its sediment transport response, producing exponential increases or decreases in mass flux as opposed to the sudden step-changes evident with the comparison models.
In cases of constant above-threshold wind velocity, mass flux predicted by the DMB model tends to an asymptote more rapidly at nearthreshold wind velocities than at very high wind velocities (Fig. 1bd). Mass flux approaches equilibrium within 5-8 s (cf. Butterfield, 1991Butterfield, , 1993Butterfield, , 1999Spies et al., 2000) when wind velocity is comparatively close to the threshold (e.g. at t = 5-25 s, t = 100-110 s), but takes longer (10-20 s) to approach equilibration when the constant wind velocity is higher (e.g. at t = 30-40 s). Using the example parameter values shown in Fig. 1, the relationship between approximate equilibration time (T E ) and wind velocity (u) over a realistic range of abovethreshold wind velocities (7-15 m s − 1 ) for the DMB model is given by: T E ≈ 1.5u − 5. This result is theoretically only applicable when flux is tending towards equilibrium, but equilibrium between wind flow and saltation at shorter timescales is not easily recognised in field data.
In the DMB model, the empirical scaling parameter a 1 controls to a large extent the magnitude of the mass flux towards which the model equilibrates for a given wind velocity (Fig. 1b). Parameters b 1 (Fig. 1c) and c 1 (Fig. 1d) are an important determinant of sediment transport near equilibrium and also in cases where wind velocity decreases gradually; this is because they represent the physical settling of grains due to gravity and horizontal wind velocity reduction, respectively. Higher values of b 1 and c 1 produce a longer lag in the exponential reduction of flux following an abrupt drop in wind velocity. Following Farrell and Sherman (2015), the universal settling velocity of a sand grain (w 0 , in m s −1 ) can be related to a reference grain size (D, in mm) by: The grain size value typically assigned for aeolian sand transport is D = 0.25 mm (e.g. Dong et al., 2003;Sherman et al., 2013), which produces a settling velocity w 0 = 1.79 m s −1 . Except in exceptional circumstances, aeolian saltation occurs within the bottom 1 m of the nearsurface boundary layer (McKenna Neuman and Scott, 1998;Nield et al., 2013). In the absence of any wind, Eq. (4) predicts that it takes 0.6 s for a grain at 1 m to fall to the ground. If all velocities in Eq. (1) are set to zero (in order to mimic settling due only to gravity) and a reference grain size is assumed, then dQ iP dt ¼ −Q iP =b 1 , which gives Q iP = Q 0 e − t/b1 as the solution over time (t). From an arbitrary initial mass flux Q 0 , Q iP is reduced by 86% (i.e. the majority of grains in the saltation cloud fall out) within the first 0.6 s when b 1 = 1/3, which is similar to the fallout rate predicted by the empirically derived Eq. (4). As a result, b 1 is set as a constant of 1/3 m s −1 in the DMB model.

Methods
Two field experiments were undertaken in order to provide collocated wind and sand transport data to fit and compare the performance of the three models described in Section 2. In this study, we differentiate the temporal resolution at which the data were measured (the 'measurement resolution') from the interval over which the data are subsequently averaged for analysis purposes (the 'averaging interval'). Experiment 1 (10 Hz measurement resolution) was conducted over a 10-min period to explore model behaviour in detail, and averaging intervals for this experiment did not exceed 10 s in order to provide enough data points to maintain statistical reliability. Experiment 2 (1 Hz measurement resolution) was designed to provide longer wind and sediment transport datasets that allowed for longer averaging intervals to be calculated. A sonic anemometer (Campbell CSAT-3) was deployed at the field site to obtain high-frequency time series measurements of wind velocity data at 0.30 m height (Fig. 2b). This measurement height accounts for the turbulent structures that are closest to the surface (and therefore most relevant for sand transport), whilst minimising interference of saltating sand grains with the sonic signal (Weaver and Wiggs, 2011). For Experiment 1, the sonic anemometer sampled three-dimensional wind velocity (uhorizontal, vspanwise and wvertical) at a frequency of 10 Hz over a period of 30 min, from which a representative 10-min period was selected. This conforms to the sampling period length recommended by van Boxel et al. (2004) for taking account of the largest turbulent eddies forming within the boundary layer. For Experiment 2, wind velocity was measured at a frequency of 10 Hz and sub-sampled to 1 Hz, over a period of 120 min. In both experiments, the anemometer was aligned parallel with the prevailing wind direction and level in the horizontal plane.

Study areas and data collection
In order to capture the high-frequency dynamics of saltation (Bauer et al., 1998), 10 Hz measurements of grain impact counts were collected using a Saltation Flux Impact Responder (Safire) (see Baas, 2004, for an extensive description of the instrument). The Safire was deployed 1.2 m away from the base of the sonic anemometer array, with its sensing ring flush with the sand surface. Four vertically separated, wedge-shaped sand traps (sampling efficiency N 90%; see Nickling and McKenna Neuman, 1997;Weaver and Wiggs, 2011) were also deployed, to calibrate saltation count with equivalent mass sand flux. The four sand traps were spaced 1.25 m apart to the side of and level with the sonic anemometer array, perpendicular to the dominant wind direction. The mass flux data collected over the experimental periods were averaged across the four traps. For Experiment 1, the sand traps were installed over the full 10-min experimental period, whilst for Experiment 2, the sand traps were installed for a 40-min period.

Field data analysis
The sonic anemometer data were post-processed to correct for any minor errors in aligning the anemometer heads into the mean approach wind flow. The streamwise frame of reference was rotated in accordance with the local streamline angle using yaw rotation (cf. van Boxel et al., 2004;Walker, 2005;Weaver and Wiggs, 2011;Chapman et al., 2012). Roll and pitch rotation were omitted from the correction procedure because the streamlines were considered to be horizontal over the flat desert surface. Subsequent to yaw rotation correction the average spanwise velocity was zero (V = 0), so no further analysis of this component was undertaken. The remaining velocity data were decomposed into their average (U; V and W ) and fluctuating (u′, v′ and w′) components, and the kinematic shear stress derived from the covariance of the instantaneous horizontal and vertical turbulent components (−u ′ w′).
For Experiments 1 and 2, the Safire datasets were used in conjunction with their respective sand trap datasets to interpolate high-frequency variations in sediment flux. Any attempt to calibrate saltation count to mass flux requires impact responder data to be collected in the same location and under the same bed conditions as sand trap data (Arens, 1996;Baas, 2004;Davidson-Arnott et al., 2009). In this study, the mass flux observed over a given time interval, Q iM (kg m −1 s −1 ), was obtained by summing the mass collected in each vertically-separated sand trap compartment (averaged across the four sand traps) over the entire experimental period (cf. Dong et al., 2003), multiplying by the proportional occurrence of saltation during the time interval (gained from the Safire data), and adjusting for sampling period and sand trap inlet size: where M c = measured mass of sediment in compartment (kg), S i = saltation count over the time interval, S tot = total saltation count over the entire experimental period, A = area of sampler's inlet (0.0075 m 2 ), t = length of time interval (s) and c represents compartment numbers. Eq. (5) assumes a linear relationship between saltation count and sand transport rate, which has been shown to be a reasonable assumption (e.g. Weaver, 2008;Davidson-Arnott et al., 2009).

Results
4.1. Fit of transport models (Experiment 1) 4.1.1. Wind and sediment transport characteristics Fig. 3a displays the time series for 10 Hz horizontal wind velocity and mass flux over the duration of Experiment 1, and Table 1 provides a summary of descriptive characteristics. The mean values for wind velocity (7.43 m s −1 ) and mass flux (0.00447 kg m −1 s −1 ) clearly mask the true turbulent nature of the transport system, with wind oscillating at a range of frequencies (10 −1 -10 1 Hz), and saltation occurring in a typically intermittent manner. Peaks in above-threshold wind velocity were often accompanied by rapid responses in sediment flux, and mass flux maxima (~0.07 kg m −1 s −1 ) were of the same order of magnitude as those recorded in previous similar studies (e.g. Dong et al., 2003, in a wind tunnel; Raats, 1996, Leenders et al., 2011, in the field).
As shown in Fig. 3b, the relationship between the instantaneous (10 Hz) horizontal component of wind flow (u′) and mass flux was moderately strong (R 2 = 0.33). Note that all coefficients of determination presented in this study are adjusted R 2 values. By sub-sampling the data to a resolution of 1 Hz, the coefficient of determination increased to R 2 = 0.48, similar to coefficients for 1 Hz data (0.34-0.71) reported in previous studies (Sterk et al., 1998;Schönfeldt and von Löwis, 2003;Leenders et al., 2005). Above-average (N0 m s −1 ) fluctuations in u ′ were primarily responsible for the majority of the larger transport events. Whilst some sand transport did occur when instantaneous horizontal wind velocities were lower than average, the mass fluxes associated with a negative u′ ( x = 0.17 × 10 −2 kg m − 1 s − 1 ; s = 0.28 × 10 −2 kg m −1 s −1 ) were significantly lower (p b 0.0001) than those associated with a positive u′ (x = 0.74 × 10 −2 kg m −1 s −1 ; s = 0.91 × 10 −2 kg m −1 s −1 ).
The relationship between the instantaneous vertical component (w ′) and mass flux (Fig. 3c) was much weaker (R 2 = 0.01) than for u′. Consequently the relationship between − u ′ w′ (a parameter that gives equal weighting to variations in both u′ and w′) and mass flux was very low (R 2 b 0.01). This has previously been noted in other studies, because sweeps (Q4) and outward ejections (Q1) offset each other in −u ′ w′ calculations due to their opposing contributions to shear stress (Leenders et al., 2005;Wiggs and Weaver, 2012;Chapman et al., 2013). Therefore, the horizontal flow component is likely to be the principal driver of entrainment and deposition, supporting the exclusion of the vertical flow component from the formulation of the DMB model.

Fitting method
The three transport models were used to fit the observed mass fluxes (Q iM ) to wind characteristics over the 10 min sampling interval of Experiment 1, using the original 10 Hz data and sub-sampled 1 Hz data. The optimisation for the fitting was performed using Sequential Least Squares Programming (SLSQP; Python® programming language), a routine based on sequential quadratic programming that allows a function of several variables to be minimised (Kraft, 1988). In order to  Table 1 Descriptive statistics of 10 Hz wind parameters at 0.30 m height and mass flux during the 10-min study period: u is wind velocity; u′ and w′ are the horizontal and vertical turbulent velocity fluctuations; −u′w′ is kinematic shear stress; Q is the mass flux over the study period. identify areas within the multidimensional parameter spaces where optimal fitting would occur, 1000-run Monte Carlo simulations were performed for each model with randomised initial conditions. The best-fit parameters for the Dong and Radok models were insensitive to initial conditions. The two-dimensional parameter space of the DMB model (with coefficient b 1 set as a constant) was somewhat sensitive to initial conditions, with some of the parameters producing more than one optimal fit, which suggests that the numerical fitting routine may be finding local minima. Therefore, the optimisation routine for the DMB model was run by setting realistic constraints on each parameter's range (based on a first round of Monte Carlo simulations) and allowing all parameters to simultaneously vary within a confined solution space. This fitting methodology was also applied to the saltation count dataset of Experiment 2. Following Ravi et al. (2012), the wind velocity threshold (u t ) for Experiment 1 was taken to be the minimum wind velocity producing saltation for N1% of the 10 min sampling interval (i.e. at least 6 s), which was determined as 6.4 m s −1 . This corresponds well to the threshold of 6.7 m s − 1 derived using the modified Time Fraction Equivalence Method (TFEM) (Stout andZobeck, 1996, 1997;Wiggs et al., 2004). Fig. 4 shows the relationships between horizontal wind velocity and observed and predicted mass flux for Experiment 1. At 10 Hz, the DMB model ( Fig. 4a) produced the best correlation between observed and predicted flux (R 2 = 0.53), compared to the Radok (R 2 = 0.42) and Dong (R 2 = 0.43) models (Fig. 4b). The spread of the data produced by the DMB model closely matched the shape of the observed data (Fig. 4a), whereas the Radok and Dong models simulated completely deterministic transport responses following exponential curves (Fig.  4b). The predictive performance of all three models improved at the 1 Hz averaging interval, but again the DMB model produced the most realistic data spread (Fig. 4c) and the best correlation (R 2 = 0.70) compared to the Radok (R 2 = 0.63) and Dong (R 2 = 0.62) models (Fig. 4d). All the coefficients of determination presented were statistically significant at the 95% confidence level. Whilst the DMB model clearly replicates the spread of data and absolute values of mass flux better than the Radok and Dong models, it is also useful to compare changes in predicted mass flux over time (Fig. 5) The three models produced time series that matched well with observed variations at both resolutions, although all models underestimated a peak in sand flux at t = 5 mins and overestimated flux by up to 0.02 kg m −1 s −1 at t = 2.5 min. The DMB model (Fig. 5c, d) tended to predict a less 'peaky' time series that more closely resembled the observed data than the Radok and Dong models, especially at the 1 Hz averaging interval. The lag that emerges naturally from the formulation of the DMB model allows it to reproduce the absence of sediment transport at lower velocities (e.g. at t = 6-7 min).

Model performance
Comparing the total modelled and observed saltation count over a given period allows for a different method of assessing model performance than simple correlation analysis. Under the assumption that the relationship between saltation count and mass flux is linearas shown by Weaver (2008) and Davidson-Arnott et al. (2009), and as is implicit in Eq. (5)the grain impact data were used to directly fit the three transport models. Saltation counts were used instead of flux rate because they can be directly summed to produce a cumulative measure of sand transport intensity. Fig. 6 displays the difference between total modelled and observed saltation count over the duration of Experiment 1, over a range of time-averaging intervals using the original 10 Hz data (1 s, 2 s, 3 s, 4 s, 5 s, 10 s). The DMB model performed well at all averaging intervals (Fig. 6a), overestimating total count by a maximum of 0.48% (at 3 s) and underestimating by a maximum of 0.37% (at 4 s). This represents a significant improvement on the performance of existing models, with the Radok model consistently overestimating total count by up to 5.50% (at 10 Hz) and the Dong model consistently underestimating count by up to 20.53% (at 10 Hz). Correlation between total modelled and observed count (Fig. 6b) generally improved as the interval increased, with the averaging out of peaks and lags in the high-frequency saltation series suiting the interval-independent formulations of the Radok and Dong models in particular. Overall, the three models performed best at the 4 s averaging interval, predicting total saltation count most accurately with the highest coefficients of determination.

Model bias
For a complete assessment of model performance, it is necessary to identify any systematic biases the models produce. The frequency distribution of residuals (Q iM -Q iP ) for the fitted model curves in Experiment 1 are shown in Fig. 7. At 10 Hz and 1 Hz (Fig. 7a,b), all three model fits produced residuals with an approximately symmetrical unimodal distribution, but strong positive skewness and high kurtosis values imply the existence of bias in each model. The residuals for the Radok and Dong models were positive for mass fluxes N0.035 kg m −1 s −1 (Fig.  7c,d), which suggests that both models systematically underestimated mass flux above this value. The DMB model systematically underestimated mass fluxes at higher levels in the 10 Hz data, around  (1 Hz), the Radok model underestimated by 73% (10 Hz) and 48% (1 Hz), and the Dong model underestimated by 69% (10 Hz) and 48% (1 Hz). These high-flux events occurred b 1% of the total experimental period during Experiment 1, but were responsible for 11.9% of the total transported sand.
Given the high-end biases identified in Fig. 7, it is important to separate the ability of the three models to simulate increasing and decreasing mass flux, as these two transport behaviours arise from differing entrainment and depositional processes. Specifically, above-threshold horizontal wind velocity likely controls the saltation system during periods of increasing mass flux, whilst the importance of vertical settling rates is likely greater during periods of decreasing mass flux. A robust local regression method (MATLAB® 'rloess' method; span 1%), which uses weighted linear least squares and a second-degree polynomial model, was employed to smooth the 10 Hz and 1 Hz curves and categorize them in discrete periods of rising or falling mass flux. At 10 Hz (Fig.  8a), the DMB model simulated mass flux equally well during periods of decreasing and increasing transport. The Radok and Dong models (Fig.  8b, c) simulated mass flux slightly more accurately during periods of decreasing transport. At 1 Hz (Fig. 8d, e, f), all three models predicted transport rates better during periods of decreasing mass flux than during periods of increasing mass flux, which suggests that depositional processes are slightly better parameterised than erosional processes. In all cases, the linear fits in Fig. 8 were found not to deviate from unity to any significant extent (using Student's t-test at the 95% confidence level). The three models are therefore statistically robust, even over the relatively short time period of Experiment 1.
Results from Experiment 1 demonstrate that, during a relatively short (10 min) experimental period characterised by intermittent saltation, the DMB model predicted sand transport better than the Radok and Dong models, both in terms of mass flux and total saltation count. This is evident from: (a) the more realistic spread of transport data produced by the DMB model at these resolutions; (b) the less 'peaky' time series due to the lag naturally emerging in the DMB model; and (c) the ability of the DMB model to predict total saltation count with far less error than the other two models, albeit with slightly lower R 2 values for averaging intervals longer than 4 s.

Testing model behaviour
The three fitted models were tested against the data of Experiment 2, in order to examine their transferability over a longer time periodwhich allows for coarser time averaging intervals to be investigatedand in a different wind/sediment transport regime to Experiment 1. The same fits derived for each model from the 1 Hz Experiment 1 mass flux dataset were used for predicting mass flux over Experiment 2 (the a 1 coefficient of the DMB model was scaled to the length of the temporal interval). The wind velocity threshold (u t ) was taken to be 6.4 m s −1 , based on the data of Experiment 1. Mass flux equivalents were only available for 40 min of the total experimental period of Experiment 2, so as for Experiment 1, the full 2-h saltation count dataset from the Safires was used to fit the three models directly.
For the 40 min of grain impact data that were calibrated to provide mass flux equivalents, the relationships between horizontal wind velocity and observed and predicted mass flux are shown in Fig. 9. At 1 Hz measurement resolution, the DMB model (Fig. 9a) produced the best correlation between observed and predicted flux (DMB model: R 2 = 0.48; Radok model: R 2 = 0.25; Dong model: R 2 = 0.36) and the most representative spread of data. Coefficients of determination improved at the 1 min averaging interval (DMB model: R 2 = 0.67; Radok model: R 2 = 0.67; Dong model: R 2 = 0.64), although the Dong and DMB models underestimated flux at wind velocities b7.5 m s −1 . These results support findings from Experiment 1 that the DMB model predicts sand flux in a more realistic manner than the other two models, although its relative performance declines over a longer averaging interval. The observed wind velocity/saltation count time series for the complete (2-h) duration of Experiment 2 is presented in Fig. 10a, alongside the model predictions. Wind velocity exceeded the threshold velocity (6.4 m s −1 ) for 98.4% of the total experimental period, and, unlike in Experiment 1, almost continuous sediment transport was recorded. The saltation count series predicted by the three models is shown in Fig.  10b, c, d. The Radok and Dong models produced similar saltation counts throughout the study period, sometimes peaking significantly at times when there was no corresponding peak in the observed data (e.g. at t = 800 s, t = 3200 s and t = 5500 s). In contrast, the DMB model predicted a less peaky time series that more closely matched the observed data. The coefficient of determination between observed and modelled saltation count was highest for the DMB model (R 2 = 0.64) compared to the Radok model (R 2 = 0.50) and the Dong model (R 2 = 0.51).   Fig. 11 compares the total observed and predicted saltation count from each model over the duration of Experiment 2, at different averaging intervals. As for Experiment 1, the DMB model performed well across all averaging intervals (Fig. 11a), overestimating total count by a maximum of 0.84% (at 10 Hz), whilst the Radok and Dong models generally did not predict total count as closely. At an averaging interval of 10 min (not shown in the figure), the DMB model overestimated total count by 0.26%, the Radok model by 0.42% and the Dong model by 4.29%. The coefficients of determination were similar for all three models (Fig. 11b), even at greater averaging periods (R 2 = 0.97 for all at the 10 min averaging interval). Whilst the DMB model provided the most accurate overall predictions of total saltation count at temporal resolutions ranging from 10 Hz-10 min, over this significant (2-h) experimental period of well-developed saltation, the difference in performance between the three models was starkest at averaging intervals of 10 s or shorter. For intervals longer than 10 s, the Radok and Dong models performed almost as well as the DMB model.

Discussion
We have shown that a sediment transport model based on dynamic mass balancing, which takes account of the turbulence-scale lag between wind flow and sediment flux, improves our capacity to predict sediment transport over a wide range of measurement resolutions and averaging intervals, from 10 Hz to 10 min. This is achieved by formulating transport in more physically realistic terms, with two main differences compared to existing models: (i) the natural emergence of a temporal lag owing to its differential structure; (ii) the explicit representation of both erosional and depositional components of sediment transport.
First, by maintaining a fraction of the mass flux from one time interval to the next, the DMB model naturally leads to lags in response of the transported sand mass to flow de/acceleration (i.e. an apparent resistance to change). This reflects the inherent inertia in the saltation system that has been observed in numerous wind tunnel and field studies (Lee, 1987;Anderson and Haff, 1988;Butterfield, 1991Butterfield, , 1993Butterfield, , 1999Spies et al., 2000;Wiggs and Weaver, 2012). Indeed, sand in a natural transport system is not immediately deposited as soon as the wind velocity decreases, because of the momentum contained in individual saltating grains. The presence of a lag resulted in stronger coefficients of determination between observed and modelled sediment flux, a much more realistic spread of mass flux data (Figs. 4 and 9),   11. (a) Proportional differences between total modelled and observed saltation count over the period of Experiment 2, at different time averaging intervals ranging from 1 Hz to 1 min. Positive values represent an overestimation of the count, and negative values represent an underestimation; (b) R 2 coefficient of determination for the relationship between total modelled and observed saltation count at time averaging intervals ranging from 1 Hz to 1 min. and less 'peaky' time series of flux (Fig. 5) and saltation count (Fig. 10), over short (10 min) and longer (2 h) experimental period lengths. Moreover, our DMB model predicted total saltation count to within 0.84% or closer for all experiments, whereas the Radok and Dong models over-or underestimated total count by up to 5.50% and 20.53% respectively. The DMB model tended to predict transport rates better during periods of decreasing mass flux (Fig. 8), which suggests that deposition processes are slightly better represented than entrainment processes in its parameterisations. Given that the saturation time T sat is the reciprocal of the sum of localised erosion and deposition rates in a system (Gao et al., 2014), further investigation of the model's behaviour is needed in order to quantify saturation limits using real field data.
Second, our aim was to formulate a model that more explicitly accounts for turbulence-scale processes, as this has been a much soughtafter tool for aeolian research. By accounting in some physical sense for both erosion and deposition in the saltation system, the DMB model provides a method for directly representing the movement of a saltating mass of sand. Crucially, the experiments presented here did not identify a temporal resolution that was too coarse for the DMB model to resolve the entrainment of grains to, and the settling of grains from, the saltation cloud. The DMB model performed well in a high-energy environment with both well-developed saltation conditions (Experiment 2) and partially developed saltation conditions (Experiment 1). This is a strong advantage given that the relationship between saltation intensity and wind speed often breaks down in cases where transport is only partially developed (see Davidson-Arnott et al., 2009).
Nevertheless, there was evidence in Experiment 1 that the DMB model systematically underestimated high mass fluxes (Fig. 7), although it did not suffer from this effect as strongly as the Radok and Dong models. Systematic underestimation may be a consequence of not explicitly considering splash entrainment, which arises in well-developed saltation conditions (Anderson and Haff, 1991;Shao and Raupach, 1992). At higher wind velocities, the additional energy contained in saltating particles could result in a power-law entrainment function that is inappropriate for lower wind velocities. Instrumentation bias at low wind velocities may have also artificially truncated the transport datasets, due to the momentum threshold in Safires for registering a sand grain impact (Baas, 2004;Stout, 2004). The introduction of a lowend bias such as this would evidently complicate the fitting of the models across the full range of observed wind velocities. However, since the distribution of wind velocities in many dryland areas broadly follows a Weibull distribution (Merzouk, 2000;Zhang et al., 2002), the occurrence of large sand transport events remains proportionally low in most natural aeolian transport systems. The three models tested here should therefore be considered statistically robust predictors of sediment transport in the majority of wind and saltation conditions.
Previous studies have proposed that suppressing the short-term variability of the transport system by using time-averaged series can provide relatively good relationships between wind velocity data and sediment flux (Wiggs et al., 2004;Baas and Sherman, 2005b;Davidson-Arnott et al., 2012). Over the short duration of Experiment 1, optimal model predictions (in terms of combining prediction accuracy for total saltation count and high correlation between observed and predicted data) were achieved when wind velocity and saltation data were averaged over a 4 s interval. This value is similar to previously reported response times of mass sand flux (e.g. Butterfield, 1991;Spies et al., 2000;Weaver, 2008). Spies et al. (2000) showed that flux responds to positive changes in wind velocity on the order of 2-3 s, whilst the response to decelerating wind is about 1 s longer. Using cross-spectral analysis, Weaver (2008) revealed a periodicity of 1.5-5 s in the horizontal wind velocity that drives the majority of sediment transport. Therefore, it seems sensible that an averaging interval on this order should allow the models to capture the essential dynamics of the transport system without being affected by sub-second variability. Our DMB model significantly outperforms the Radok and Dong models at this scale of analysis.
However, over the longer dataset of Experiment 2, the best results were achieved at 1-min averaging intervals. The 'averaging out' of peaks and lags in the transport system explains why the Radok and Dong models, which are interval-independent, were more suited to this resolution of analysis than at shorter resolutions. The DMB model still slightly outperformed the Radok and Dong models at the 1 min interval, possibly because longer periodicities in the sand transport system are accounted for by its differential formulation. On the basis of these results, we propose that in cases where the temporal scale of investigation is short (on the order of minutes), sand transport should be predicted using our DMB model run at a 4-s averaging interval. However, when the lengths of wind/transport datasets are longer (on the order of hours), sand transport can be accurately predicted at 1-10 min averaging intervals, with only slight differences in model performance between the DMB, Radok and Dong models.
This has ramifications for how wind and sand transport data are collected in the field, because it informs the appropriate resolution at which data should be measured and therefore the type of anemometry that is best suited to a given scale of investigation. Our recommendations also imply that turbulence-scale processes linked to the saltation lag, which are more explicitly accounted for in our DMB model compared with existing models, can be integrated into longer-term simulations of sand transport without substantially increasing computational expense. The approach we propose could prove to be significant for integrating turbulent transport processes into macro-scale landscape modelling of drylands.

Conclusions
The new dynamic mass balance (DMB) model presented in this study represents an improved method for predicting both high-frequency sand transport variability and coarser-scale total mass transport. The temporal lag emerging from the DMB model accounts for the momentum contained in individual saltating sand grains. This results in better predictions of sand transport than the models of Radok (1977) and Dong et al. (2003) in terms of accurately simulating mass flux distributions and time series of flux and saltation count. This has been shown to be the case over a wide range of measurement resolutions/averaging intervals (10 Hz-10 min). For all experiments and averaging intervals presented in this study, the DMB model predicted total saltation count to within 0.48% or closer, whereas the Radok and Dong models over-or underestimated total count by up to 5.50% and 20.53% respectively. The best correspondence between observed and modelled sand transport, in terms of accurately predicting total saltation count whilst maintaining high correlation with the observed data, was achieved when the DMB model was run using data averaged over 4 s (for Experiment 1) or 1 min (for Experiment 2). We therefore propose that shortterm, local-scale aeolian transport research should be conducted at a 4-s measurement resolution if the appropriate data collection methods are available, and employ our DMB model for predictive purposes. Over longer timescales, running our DMB model with averaging intervals of 1-10 min also yields slightly more accurate predictions of sand transport than the existing models analysed here. This could prove to be a revolutionary way of incorporating fine-scale transport dynamics into macroscale studies of dryland landscape change. methodological discussions and for developing the fitting optimisation routine. Two anonymous reviewers are thanked for their constructive comments that greatly helped to improve this paper.