Lagrangian modeling of bedload movement via the impulse entrainment method

At low shear stresses merely above incipient conditions, the characteristic travel patterns of bedload particles remain difficult to predict due to intermittent movement. At these conditions, particles rest, are entrained into the flow when sufficient momentum transfer from near-bed turbulence occurs, and then come to rest again. The overarching goal of this research is to close the critical gaps related to intermittent movement in order to allow for Lagrangian modeling of bedload at near incipient conditions. Thereby, the specific objectives of the present work are to predict the statistics of (1) the particle resting time, tR, and (2) the magnitude of hydrodynamic momentum transfer (or impulse) during entrainment, Ient. To predict these statistics, we employed the conceptual framework of the impulse entrainment method and predicted impulse statistics by simulating turbulent time series realizations with a generic, regime-based streamwise velocity spectrum. Model validation was carried out by directly comparing simulation results with published experimental impulse and particle entrainment statistics. Model predictions showed that an increase in stress was correlated with a sharp decrease in the average tR and an increase in I ent .


Introduction
The characteristic travel patterns of bedload particles have been shown to influence mean bedload fluxes and also the fluctuations in bedload transport. Prediction of the former magnitude has provided one motivation for investigating the characteristics of particle travel [1][2]. Understanding of the latter temporal fluctuations has also advanced recently through detailed experimental investigations focused on particle travel and particle-toparticle interactions [3]. The characteristics of particle travel follow a complex relationship with bed shear stress, however, because of the gradual regime transition that occurs from low intensity bedload transport at stresses near incipient motion to high intensity transport and general motion at higher shear stresses [4]. At the particle scale, the transition from low to high intensity bedload transport has been shown to correspond with transitions from intermittent to continuous motion [5] and from primarily rolling or sliding transport (i.e., travel in regular contact with the bed) to primarily saltation (i.e., particle hopping with less frequent collisions with the bed) [6].
In the regime of high intensity transport, where bedload particles primarily travel in saltation, numerical Lagrangian particle tracking models have been used to successfully predict the statistics of geometric and kinematic particle travel characteristics [1][2]. In these models, numerous particle hops are typically simulated by estimating the governing hydrodynamic forces with drag and lift force equations and then tracking particle trajectory with the change in linear momentum equation. While this type of framework may apply similarly to the transport phase of a moving particle in the low intensity transport regime, the inability of these models to account for particle resting precludes their usage for predicting intermittent transport characteristics without modification. Therefore, a new simulation component is needed in order to predict the characteristics of both particle resting and entrainment. The problem of predicting particle entrainment has been studied extensively, though exact prediction remains challenging even in controlled laboratory settings. Some semitheoretical entrainment models [7][8] consider that there is a critical hydrodynamic force above which entrainment will occur, then relate instantaneous velocity to hydrodynamic force, and finally use the probability density function (PDF) of turbulent velocity fluctuations to predict the probability of entrainment. These models have generally performed well with some fitting of coefficients to predict the probability of entrainment. However, in a pivotal paper highlighting the direct relation between instantaneous velocity and entrainment events, [9] showed that not only the magnitude but also the duration of strong flow events is important for predicting the moments of particle entrainment. The instantaneous impulse magnitude was therefore used to predict entrainment, which is defined as the integral of the streamwise hydrodynamic force, Fx, over a given time interval when Fx continuously exceeds the critical force required to first move the particle, Fcr. By considering Fx to be proportional to the square of the instantaneous streamwise velocity, u 2 , via the drag force equation, the magnitude of an individual impulse event, Ii, can be defined as: where u 2 cr is the critical square velocity associated with Fcr, the condition u 2 > u 2 cr is satisfied continuously between t1 and t2, <u 2 >i is the time-averaged square velocity magnitude (or force magnitude) of the impulse event, and Ti (= t1-t2) is the duration of the impulse event [10]. According to the impulse entrainment method, individual impulse events are consider to result in particle entrainment if Ii is greater than a critical impulse value, Icr. Using this impulse entrainment method, experimental studies have been conducted to investigate the statistics of Ii and its sub-components as well as determine the value of Icr for certain conditions flow and sediment conditions [10][11].
Despite the established correlation of high magnitude impulse events with particle entrainment, limited work has been done to predict the statistics of impulse events for wider ranges of flow and sediment conditions. From experimental data, [11] investigated the PDFs of Ii and suggested that it is lognormally distributed. Their work, however, spanned a limited range of flow conditions and was only carried out for one pair of mobile particle and bed material sizes. Therefore, the primary goal of this research is to develop a method for predicting the PDF of turbulent impulse events in order to predict the statistics of (1) the particle resting time, tR, and (2) the magnitude of hydrodynamic momentum transfer (or impulse) during entrainment, Ient.

Model Development
In order to predict the impulse distribution, we consider that it will be governed, in general, by two probability groups: (1) the probability that u 2 will exceed u 2 cr; and, (2) the conditional probabilities that successive velocities will obtain u 2 values that exceed u 2 cr given the u 2 values of the previous time steps that also exceeded u 2 cr. In simple terms, the first probability group may be expected to control the statistics of u 2 and therefore the percentage of total time occupied by impulse events. The second probability group may be expected to control the distribution of impulse event durations. The notion expressed in the first probability group has been incorporated into previous semi-theoretical probability of entrainment models by using turbulent flow PDFs [7][8]. However, the second probability group is more complex due to the many scales of motion in a turbulent flow and to our knowledge has not been considered before in a semi-theoretical model.
In order to appropriately account for the second probability group, we recognize the correspondence of the described conditional probabilities with the correlational statistics commonly used in time series analysis, namely, the autocorrelation function. Furthermore, we recognize that the autocorrelation function is directly related to the velocity spectrum via the Fourier transform. With this reasoning, we hypothesize herein that by simulating a turbulent streamwise velocity time series realization using the streamwise velocity statistics (mean and turbulent intensity) and the streamwise velocity spectrum, we will be able to appropriately predict the impulse distribution and the distributions of its sub-components.
To generate the streamwise velocity spectrum, Suu, we rely in general upon the work of [12], wherein a generic, four-range spectrum model was delineated as a function of the flow depth, vertical elevation and roughness scales. Of the four ranges considered in [12], the inertial subrange (i.e., "-5/3" power law range) was ignored herein based on findings from experimental investigations on near-bed velocity spectra [13]. The three regimes considered herein correspond to the largest eddies (i.e., the "0" power law range), the intermediate eddies (i.e., the "-1" power law range) and the viscous subrange, for which a "-7" power law was assumed [12][13]. Locations of the breakovers from the "0" to "-1" power law ranges and from the "-1" to "-7" power law ranges were defined to occur at wavenumbers kw1 = (a1H) -1 and kw2 = (a2ks) -1 , respectively, where H denotes the flow depth, ks denotes the characteristic bed roughness height, and a1 and a2 are coefficients [12]. Suu was calculated for a range of wavenumbers, kw, using the considered power law ranges with coefficients selected to ensure that the variance from Suu was equal to the squared value of the desired streamwise turbulent intensity, σu 2 . The Suu values were then converted to the frequency, f, domain using f = 2πkw/Uc, where Uc denotes the eddy convective velocity, which was determined as a function of kw, the bulk flow velocity, Ubulk, and the local time-averaged velocity, u [14]. An example Suu spectrum is shown in Figure 1a, where f1 and f2 correspond to kw1 and kw2.
Employing the generated Suu spectrum, the inverse fast Fourier transform (IFFT) was used to quickly simulate a time series realization of the fluctuating component of the instantaneous streamwise velocity, u′ [15]. Finally, the instantaneous velocity time series, u, was determined from u′ and u by using the typical decomposition of instantaneous velocity into fluctuating and time-averaged components, i.e. u = u + u′. A snapshot of a simulated u 2 time series is presented in Figure 1b. An arbitrary u 2 cr limit for impulse events is also superimposed in Figure 1b and the characteristics of an example impulse event are illustrated. Table 1 summarizes the input parameters required for the present model and describes the purpose(s) of each parameter.

Model Validation
Validation of the proposed impulse model was carried out using published experimental data from [11], wherein impulse distributions were studied and directly connected with the probability of particle entrainment. The most relevant model output statistics have been used for performance assessment. These are the mean and standard deviation (STDEV) of the impulse magnitude (Ii) and the impulse sub-components, which are the average force magnitude (<u 2 >i) and impulse duration (Ti). In addition, the number of predicted particle entrainments (nent) in a given time period was considered. The values for all model input parameters (i.e., Table 1) were available in [11] except for the coefficients a1 and a2. Model optimization was therefore conducted by adjusting a1 and a2 in order obtain good agreement between experimental and simulation statistics for one condition presented in [11]. Values of a1 = 2.0 and a2 = 0.18 were fitted, are within the range of expectation based on [16] and were maintained constant for all other conditions reported in [11] to assess prediction ability of the present model. The full comparison of experimental and simulation statistics is provided in Figure 2, where the agreement is very good for all metrics (< 20% error) except for very low magnitude values (i.e., STDEV of <u 2 >i and low values of nent). Figure 3 presents histograms of u 2 , <u 2 >i, Ti and Ii that were generated from simulations at conditions corresponding exactly to test E1 reported in [11]. An almost identical figure was published in [11] 3. Histograms of u 2 , impulse sub-components (<u 2 >i and Ti) and impulses (Ii) for condition E1 at τ * = 0.011 in [11]. Figure is directly comparable to Figure 5 in the manuscript of [11].

Model Predictions
In order to aid Lagrangian tracking of an individual particle through intermittent movement and resting, the primary outputs of concern from the present model are the resting time between particle entrainment events, tR, and the magnitude of hydrodynamic momentum transferred to the particle during entrainment, Ient. Therefore, multiple simulations were conducted for all of the flow conditions reported in [11] in order to identify the times of particle entrainments and the impulse magnitudes that caused them. In total, 1000 samples of both tR and Ient were collected. Figure 4 presents the mean and STDEV of tR and Ient as a function of the dimensionless bed shear stress, τ * = τ0/[(ρs-ρ)gd50]. Figure 4a shows that model simulations estimate both the mean and STDEV of tR to decrease rapidly with increasing τ * , as shown through the approximately "-8" power law. Furthermore, the mean and STDEV of tR have comparable magnitude, showing a wide range of variability in tR for a given τ * condition. Figure 4b shows that Ient statistics increase gradually with increasing τ * and that the STDEV of Ient increases more quickly than the mean.

Fig. 4.
Mean and STDEV of (a) the particle resting time, tR, and (b) the impulse at particle entrainment, Ient, for the conditions reported in [11].

Conclusions
To fill the critical gaps remaining before Lagrangian particle tracking models can predict the characteristics of intermittent particle movement, a model for predicting the timing and hydrodynamic characteristics of particle entrainment events has been developed herein. The model is based upon the conceptual framework of the impulse entrainment method, which suggests that both the magnitude and duration of strong turbulent events are important for predicting particle entrainment. To predict impulse statistics, we simulate realizations of turbulent time series with reasonable temporal characteristics over the range of turbulent scales by employing a generic, regime-based streamwise velocity spectrum.
Model validation was carried out by directly comparing simulation results with the experimental results from [11]. The differences between simulation and experimental impulse statistics are less than 20% except for very low magnitude values. The model outputs were analyzed to predict the statistics of particle resting time, tR, and hydrodynamic momentum transfer at entrainment, Ient. It was found that an increase in shear stress was correlated with a rapid decrease in tR and an increase in Ient.