Numerical Analysis of Motion Characteristics of Sliding or Rolling and Saltation of Sediment Particles under Turbulent Flow

: The processes of sediment particle movement were studied through numerical simulation using a coupled method with focus on discussing the characteristics of sliding or rolling and saltation sediment particles, respectively. Turbulent ﬂow was simulated using large eddy simulation (LES). The sediment particle was simulated using the combined ﬁnite-discrete element method (FDEM). The interaction forces of turbulent ﬂow and sediment particle were calculated using the immersed boundary method (IBM). It indicated that the collisions of saltating particle with low concentration increase the saltation length and ﬂight time. In response, sediment particle velocity also increases. The particle angular velocity is largest at the takeoff moment, and decreases gradually in the saltation progress. The drag and lift forces near the bed are large, and away from the bed decrease and trend to be a stable value, gradually. From the relative magnitudes of the drag and lift forces, the lift force plays a more important role than the drag force in the sediment saltation. The relative magnitudes of drag and lift forces inﬂuence the incident and takeoff angles. The sediment transport rate calculated based on the characteristics of saltation sediment particles is overestimated, ignoring the effect of sliding or rolling sediment particles and inter-particle collisions.


Introduction
Bed load movement is one of the most important sediment transport phenomena and can generally be classified into the three different modes, namely, rolling, sliding, and saltation. The different transport types depend on the sediment size, density, and flow condition among others. Above a certain threshold shear stress, these transport types may occur simultaneously. The insight into the mechanism of bed load movement is crucial to river, estuarine, and coastal processes, such as sediment transport, morphological change, structure erosion, and deposition.
Bed load transport has been investigated extensively over the decades via different approaches. Einstein [1] found bed load movement within a thin layer of two particle diameters above the bed and the saltation length to be about 100 times the particle size. Bagnold [2] described bed load as the unsuspended transport of particles above the sediment bed and found the sediment saltation to be dominated by the gravity force and not affected by the effect of turbulence and hydrodynamic lift force. The statistical characteristics of saltation sediment particle over a fixed bed such as saltation length, height, streamwise particle velocity, incidence, and takeoff angles of sediment particle collision with the bed were analyzed by experimental studies [3][4][5][6][7][8][9]. Based on the experimental data, numerical model, statistical features of flow on a rough bed, processes of sediment particle movement, and the distribution of statistical characteristics of sediment particle movement. "Results and Discussion of Sediment Particle Movement" presents the comparative analysis of the model and experimental results. Finally, the section titled "Conclusions" summarizes the conclusions.

Numerical Method
The simulation of bed load movement on the fixed sediment bed involves the simulation of turbulent flow, the particle movement, the deformation and collision with bed, and the interactions between fluid and particle. In this study, the turbulent flow is simulated using the named CgLes computational fluid dynamics (CFD) code [44], which is a three-dimensional grid-based code that solves the Navier-Stokes equations using a finite volume method and has a second-order numerical accuracy in both time and space. This code has been applied to simulate the turbulent flow over a rough bed [45] and has a good robustness and parallelizing efficiency. To simulate the movement, deformation, and collision of particles, the combined finite-discrete element method (FDEM) proposed by Munjiza et al. [46] is used. The interaction forces between fluid and particle are calculated by the immersed boundary method (IBM) [47]. The hydrodynamic forces exerted on the particle cause the particle movement; in return, the particle forces exerted on the fluid are then added into the momentum equation of fluid motion in order to achieve the fluid-particle interaction.

Governing Equations of Flow
In the Cartesian coordinate system, the continuity equation and the incompressible fluid Navier-Stokes equations using the IBM with the second-order Adams-Bashforth temporal-discretized are expressed as ∇ · u n+1 = 0 (1) where u is the fluid velocity, h is the convective and diffusive terms, υ is the eddy viscosity coefficient, δt is the time step length, ∇ is the gradient operator, t is the matrix transposition, n−1, n, n + 1/2, n + 1 is the time step, and f is the fluid-particle interaction forces. The pressure Poisson equation is solved using a multigrid preconditioner and conjugate gradient method.

Mathematical Model of Particle Motion
In the present studies, the sediment particle is regarded as a discrete element, and the deformability of discrete particles and the interactions between discrete particles are represented by a set of contacting finite elements. When two discrete particles are in contact and penetrating each other, the magnitude of the contact forces is dependent on the overlapping area of the finite elements. In the present study, the combined finitediscrete element method [40] for contact detection, deformability, and interaction between particles is adopted; the algorithmic procedures for fracture and fragmentation algorithms of individual particles in the FDEM code are switched off considering that deformation of sediment particles is in the rigid grain limit. A distributed contact force algorithm introduced by Munjiza and Andrews [48] is applied to calculate realistic contact forces over a finite contact element.
Generally, the contact forces between two discrete particles or between discrete particles and the wall can be separated into the normal force F n and the tangential force F s . The normal contact force is calculated by the penalty function method proposed by Munjiza and Andrews [48], the expression is given as where β ci and β tj represent the i th and j th finite elements of the contactor and target discrete elements, respectively; n and m represent the total number of discretized finite elements of the contactor and target discrete elements, respectively; n Γ β ci ∩β tj represents the outward unit normal to the boundary Γ of the overlapping volume β ci ∩ β tj ; and ϕ ci and ϕ tj represent potential functions for the contactor and target discrete elements, respectively. In the present study, the simplest finite element of linear four-node tetrahedron has been adopted. The potential function at an arbitrary point p inside a tetrahedral finite element is calculated by where k p is the penalty parameter and is approximately equal to 10E; E is the Young's modulus. V is the volume of the tetrahedral finite element, and V i (1,2,3,4) is the volume of the corresponding sub-tetrahedrons at the point p. If the point p is at the tetrahedron center, the potential function ϕ equals to 1; if the point p is at the tetrahedron surface, the potential function ϕ equals 0. The tangential contact force in the FDEM by taking into account the sliding friction force is calculated using the formula proposed by Xiang et al. [41], which is expressed as where k s = Ed/2 is the tangential spring stiffness constant, and Roux and Combe [49] pointed out that the stiffness constants can be chosen with the correct order of magnitudes, comparing them to typical estimated values in contacts under some given stress level, which is a mere computational trick to forbid grain interpenetration. To save computational cost, a much smaller value (10 MPa) for Young's modulus is used in this simulation. Schmeeckle [29] simulated the sediment transport of medium sand and investigated the effect of sediment movement on the turbulent flow using the value (0.5 MPa) for Young's modulus. Simeonov and Calantoni [50] pointed out that the dynamics of granular flows are not very sensitive to the value of E in the range of 0.1-12 MPa and lead to a difference in the stiffness coefficient greater than 4%. Furthermore, η d is the viscous damping coefficient, η d = 2h e Eρ s ; h e is the smallest size of the finite elements; and δ s and V s are the tangential relative displacement and velocity between particles, respectively. When the two particles are in contact with each other, if the tangential elastic contact force is bigger than the friction force, the tangential force is then calculated using the total normal elastic contact force F n based on Coulomb friction law, which is defined as where the µ is the sliding friction coefficient, which is adopted as 0.4.

The Interaction of Flow and Particle
The sediment particles are discretized into finite elements, and the interactions between flow and discrete particles are modeled by the immersed boundary method. The two-way coupling forces between the finite element points on the surface of the sediment particle and the Cartesian grid points are achieved by interpolation/distribution functions (I/D). The interpolation function I projects the physical field from the Cartesian grids to the immersed boundary points. On the other hand, the distribution function D maps the physical field from the immersed boundary points back to the Cartesian grids. The added mass force f exerted on the Cartesian grids is expressed as where D and I represent the distribution and interpolation functions, respectively; V presents the velocity of the immersed boundary points on the surface of the sediment particle. The variable Φ(X i ) (such as U, F, etc) on the immersed boundary points is calculated as where g h represents the set of Cartesian grids both inside and outside of individual sediment particle; and δ h presents the discrete delta function proposed by Peskin [51].
The variable φ(x) on the Cartesian grids (such as u, p, f, etc.) both inside and outside the immersed boundary is calculated as where N ibp represents the total number of immersed boundary points; and ∆V i represents the discrete volume around the immersed boundary points X i . Uhlmann [52] suggested that these volumes should form a thin shell of thickness equal to one mesh width around each immersed boundary point. For a more detailed calculation procedure of the interaction between solid and flow using the immersed boundary method, readers should refer to the paper of Ji et al. [47].

Setup of Numerical Model
The numerical simulation of sediment particle movement on plane sediment bed in the present study was carried out in a three-dimensional (3D) domain, with the Cartesian coordinates x, y, and z aligned in the streamwise, vertical (upward positive), and spanwise directions, respectively ( Figure 1). In order to investigate the different modes of sediment particle movement on the rough bed, a layer of sediment particle was fixed and arranged in a hexagonal lattice in the bottom of the computational domain. The no-slip boundary condition was applied to the bottom of the computational domain; a free-slip rigid lid boundary condition was applied to the top of the fluid domain; and periodic boundary conditions were imposed in the upstream, downstream, and cross-stream sides of the computational domain. Considering the large amount of computation costs, and the periodic boundary conditions of the streamwise and lateral side in the model, the streamwise, vertical, and spanwise length of the computational domain were set as 42 √ 3 d, 15 d, and 18 d, respectively. The setup of the numerical model was based on the setup of other similar numerical models, which is that the width of the computational domain is slightly larger than the height [53]. computational domain. Considering the large amount of computation costs, and the peri-odic boundary conditions of the streamwise and lateral side in the model, the streamwise, vertical, and spanwise length of the computational domain were set as 42 3 d, 15 d, and 18 d, respectively. The setup of the numerical model was based on the setup of other similar numerical models, which is that the width of the computational domain is slightly larger than the height [53].

Sediment Transport and Hydraulic Parameters of Numerical Model
The motion characteristics of the sediment particle transport depend on a variety of parameters, such as the density and shape of the moving particles, solid volume fraction, Shields number, particle Reynolds number, Galileo number, and bed structure [14,53,54]. In this study, the sediment particles were assumed to be the "smooth" sphere with a diameter of d = 1 mm; the density was 2.650 g/cm 3 . The net water depth f h was 14.25 d. There were 100 particles in the computational domain for statistical convenience, and the global solid volume fraction was set to

Sediment Transport and Hydraulic Parameters of Numerical Model
The motion characteristics of the sediment particle transport depend on a variety of parameters, such as the density and shape of the moving particles, solid volume fraction, Shields number, particle Reynolds number, Galileo number, and bed structure [14,53,54]. In this study, the sediment particles were assumed to be the "smooth" sphere with a diameter of d = 1 mm; the density was 2.650 g/cm 3 . The net water depth h f was 14.25 d. There were 100 particles in the computational domain for statistical convenience, and the global solid volume fraction was set to φ s = 8 × 10 −5 . The important parameters were calculated as follows: the friction velocity was calculated as u * = τ/ρ f , where τ is the shear stress estimated using a linear extrapolation of the total shear stress profile to the rough channel bed [45]. The Shields number was defined as θ = τ/((ρ s − ρ f )gd); the bulk velocity u b = h y b u f dy/h f ; the Reynolds number based on bulk flow velocity Re b = u b h f /ν; and the Reynolds number based on friction velocity Re τ = u * h f /ν, where ν is the fluid kinematic viscosity. The particle Reynolds number was calculated as Re p = u * d/ν. The Galileo number was defined as Ga = ρ s /ρ f − 1 gd 3 /ν 2 . The computational domain was discretized by a uniform isotropic grid with grid spacing in terms of wall units ∆x + = ∆y + = ∆z + = ∆xρu * /ν. Table 1 summarizes the important physical and numerical parameter values adopted in this study.

Computational Procedure of Numerical Model
The computational procedure of the simulated cases involved two stages. In the first stage, the sediment particles were fixed at their positions, and the streamwise gravity component g x was imposed on the fluid to generate fully developed turbulence in the computational domain. In the second stage, the sediment particles were released to move under the hydrodynamic forces. The motion trajectories and the hydrodynamic forces exerted on particles were tracked at the different flow intensities.

Validation of Numerical Models
The mean streamwise velocity profile normalized with the bed shear velocity u * is presented in Figure 2. Y + is the normalized effective height (y − y b )u * /ν, y b is the effective location of the bed, which is normally estimated by fitting the mean velocity curve to the logarithmic law-of-the-wall region. Grass et al. [55], Defina [56], and Singh et al. [45] determined equivalent sand roughness k s = 79.2, 288, and 102 by fitting mean velocity profiles to the logarithmic law-of-the-wall. Using a similar approach, equivalent sand roughness k s was equal to 110 in the present study. Through a comparison with the experimental results of Defina [56] and Grass et al. [57], it is indicated that the flow profile followed the standard pattern of velocity defect increasing with wall roughness.

Computational Procedure of Numerical Model
The computational procedure of the simulated cases involved two stages. In the firs stage, the sediment particles were fixed at their positions, and the streamwise gravity component gx was imposed on the fluid to generate fully developed turbulence in the computational domain. In the second stage, the sediment particles were released to mov under the hydrodynamic forces. The motion trajectories and the hydrodynamic forces ex erted on particles were tracked at the different flow intensities.

Validation of Numerical Models
The mean streamwise velocity profile normalized with the bed shear velocity * i presented in Figure 2. Y + is the normalized effective height y is the ef fective location of the bed, which is normally estimated by fitting the mean velocity curv to the logarithmic law-of-the-wall region. Grass et al. [55], Defina [56], and Singh et al. [45 determined equivalent sand roughness ks = 79.2, 288, and 102 by fitting mean velocity pro files to the logarithmic law-of-the-wall. Using a similar approach, equivalent sand rough ness was equal to 110 in the present study. Through a comparison with the experi mental results of Defina [56] and Grass et al. [57], it is indicated that the flow profile fol lowed the standard pattern of velocity defect increasing with wall roughness. The Reynolds shear stress normalized by the bed shear stress is shown in Figure 3 together with the total shear stress on the smooth bed. The profile of the Reynolds shea stress presented a straight line away from the bed and no significant deviations compared with the total shear stress profile for the smooth. It is indicated that the turbulent flow wa indeed fully developed on the rough particle bed. The Reynolds shear stress normalized by the bed shear stress is shown in Figure 3, together with the total shear stress on the smooth bed. The profile of the Reynolds shear stress presented a straight line away from the bed and no significant deviations compared with the total shear stress profile for the smooth. It is indicated that the turbulent flow was indeed fully developed on the rough particle bed.
Water 2022, 14, x FOR PEER REVIEW The turbulence intensities normalized by the bed shear velocity are shown i 4, together with the experimental results from Grass [57] and Nezu [58], and the n simulation results from Singh et al. [45]. The simulated turbulence intensities in  The turbulence intensities normalized by the bed shear velocity are shown in Figure 4, together with the experimental results from Grass [57] and Nezu [58], and the numerical simulation results from Singh et al. [45]. The simulated turbulence intensities in three directions are also in good agreement with both experimental and DNS data. The turbulence intensities normalized by the bed shear vel 4, together with the experimental results from Grass [57] and Nez simulation results from Singh et al. [45]. The simulated turbulen rections are also in good agreement with both experimental and

Processes of Sediment Particle Movement
The sediment particle movement near the bed under a turbu complex problem. In this study, the high-precision numerical vestigate the dynamic behaviors and the corresponding hydrod derived from the average value of the surface points on the se ing or rolling and saltation sediment particles. Figure 5 show instantaneous contours of the streamwise velocity, the Q-criterio herent structures, and the particles distribution. It is indicated th close to the bed are induced to move under the velocity fluctuat and coherent structures.

Processes of Sediment Particle Movement
The sediment particle movement near the bed under a turbulent flow is an inherently complex problem. In this study, the high-precision numerical models were used to investigate the dynamic behaviors and the corresponding hydrodynamic forces (they are derived from the average value of the surface points on the sediment particle) of sliding or rolling and saltation sediment particles. Figure 5 shows the model results of the instantaneous contours of the streamwise velocity, the Q-criterion eddy iso-surface of coherent structures, and the particles distribution. It is indicated that the sediment particles close to the bed are induced to move under the velocity fluctuation of high and low strip and coherent structures. Figure 6 presents the typical continuous trajectories and dynamic behaviors of individual particle movement and the corresponding hydrodynamic forces exerted on sediment particles with a resolution of 0.00005 s. From the trajectory variations of sediment particle movement, it clearly indicates that a sediment particle moves forward in the mode of sliding or rolling and saltation; the different motion modes are intermittent with time. From the variation of the dynamic responses of sediment particle movement, it indicates that the streamwise and vertical velocity of the particle are disturbed under the instantaneous impulses of hydrodynamic forces during coherent structures. From the variations of the hydrodynamic forces exerted on sediment particles, it shows that the hydrodynamic forces of the moving particle present relatively obvious changes during the process of the moving particle's contact with the bed. At the moment of sediment particle saltation, the absolute value of the drag force increases drastically then decreases gradually. The lift force plays a more important role than the drag force in the saltation process. Particle rotation is mainly controlled by the particle's collisions with the bed. From the variation of the particle's angular velocity (the clockwise rotation along the z direction is positive), when the particle collides with the fixed particle, the particle's takeoff angular velocity is obviously larger than the particle's incidence angular velocity, and the particle's angular velocity increases drastically and reaches the maximum value, then decreases gradually owing to the viscous effects.  Figure 6 presents the typical continuous trajectories and dynamic behaviors of i vidual particle movement and the corresponding hydrodynamic forces exerted on s ment particles with a resolution of 0.00005 s. From the trajectory variations of sedim particle movement, it clearly indicates that a sediment particle moves forward in the m of sliding or rolling and saltation; the different motion modes are intermittent with t From the variation of the dynamic responses of sediment particle movement, it indic that the streamwise and vertical velocity of the particle are disturbed under the insta neous impulses of hydrodynamic forces during coherent structures. From the variat of the hydrodynamic forces exerted on sediment particles, it shows that the hydrodyna forces of the moving particle present relatively obvious changes during the process o moving particle's contact with the bed. At the moment of sediment particle saltation absolute value of the drag force increases drastically then decreases gradually. The force plays a more important role than the drag force in the saltation process. Par rotation is mainly controlled by the particle's collisions with the bed. From the varia of the particle's angular velocity (the clockwise rotation along the z direction is posit when the particle collides with the fixed particle, the particle's takeoff angular veloci obviously larger than the particle's incidence angular velocity, and the particle's ang velocity increases drastically and reaches the maximum value, then decreases gradu owing to the viscous effects. Bed loads move randomly and intermittently in sliding or rolling and saltation modes under the turbulence flow. Böhm et al. [59] and Frey [60] studied the bed load motion modes on the rough bed with three particle layers on average by experiment and distinguished the rolling and saltation of the sediment particle based on the threshold velocity (u 0 ) and distance (d n is the measured distance between the mass centers from the particle to the next neighbor particle based on averaged five consecutive frames). The parameter values of the threshold velocity (u 0 = 0.025 m/s) and distance (d n /d < ε, ε < 1.07) between the particles were determined by trial and error to minimize the differences between the different motion state. Auel et al. [61] investigated the isolate sediment transport modes on the fixed bed and differentiated the sediment particle rolling and saltation according to threshold height (if particle center exceeded a threshold distance (0.6 d) away from the bed surface, the particle was assumed to change its mode from rolling to saltation). Based on the studies [59][60][61], the threshold height and distance (Figure 7) are adopted to distinguish the sediment particle motion modes in this study. The two threshold parameters are determined compared with the trajectories and hydrodynamic forces of sediment particle movement ( Figure 6). If the vertical distance d v between moving particles and fixed particles is less than 1.1 d (corresponding threshold distance (0.6 d) of Auel et al., 2017 [61]), and the horizontal distance d l between the same particle's continuous contact with fixed particles is less than 1.74 d, then the sediment particle is assumed to be a rolling or sliding particle (SRP). If the vertical distance d v between moving particle and fixed particle is larger than 1.1 d, the sediment particle is assumed to be a saltation particle (SP). The distribution of statistical characteristics of different motion modes of sediment particles are investigated based on the simulated results, including the velocity and angle velocity, the saltation height and length, and the incidence angles and takeoff angles of sediment particles.
Water 2022, 14, x FOR PEER REVIEW 10 of 22 Figure 6. Typical continuous trajectories and dynamic behaviors of individual sediment particle movement (sliding, or rolling and saltation), and the corresponding hydrodynamic forces exerted on the particle. p x and p y = streamwise and vertical center coordinates of the particle, respectively; p u and p v = streamwise and vertical velocities of the particle, respectively; D F and L F = drag and lift forces of hydrodynamic forces exerted on the particle, respectively; x s = angular velocity of the particle; and G = particle submerged weight.
Bed loads move randomly and intermittently in sliding or rolling and saltation modes under the turbulence flow. Böhm et al. [59] and Frey [60] studied the bed load motion modes on the rough bed with three particle layers on average by experiment and distinguished the rolling and saltation of the sediment particle based on the threshold velocity (u0) and distance (dn is the measured distance between the mass centers from the particle to the next neighbor particle based on averaged five consecutive frames). The parameter values of the threshold velocity (u0= 0.025 m/s) and distance (dn/d < ε, ε < 1.07) between the particles were determined by trial and error to minimize the differences between the different motion state. Auel et al. [61] investigated the isolate sediment transport modes on the fixed bed and differentiated the sediment particle rolling and saltation according to threshold height (if particle center exceeded a threshold distance (0.6d) away from the bed surface, the particle was assumed to change its mode from rolling to saltation). Based on the studies [59][60][61], the threshold height and distance (Figure 7) are adopted to distinguish the sediment particle motion modes in this study. The two threshold parameters are determined compared with the trajectories and hydrodynamic forces of sediment particle movement ( Figure 6). If the vertical distance dv between moving particles and fixed particles is less than 1.1 d (corresponding threshold distance (0.6 d) of Auel Figure 6. Typical continuous trajectories and dynamic behaviors of individual sediment particle movement (sliding, or rolling and saltation), and the corresponding hydrodynamic forces exerted on the particle. x p and y p = streamwise and vertical center coordinates of the particle, respectively; u p and v p = streamwise and vertical velocities of the particle, respectively; F D and F L = drag and lift forces of hydrodynamic forces exerted on the particle, respectively; s x = angular velocity of the particle; and G = particle submerged weight. , and the horizontal distance dl between the same particle's continuous contact with fixed particles is less than 1.74 d, then the sediment particle is assumed to be a rolling or sliding particle (SRP). If the vertical distance dv between moving particle and fixed particle is larger than 1.1 d, the sediment particle is assumed to be a saltation particle (SP). The distribution of statistical characteristics of different motion modes of sediment particles are investigated based on the simulated results, including the velocity and angle velocity, the saltation height and length, and the incidence angles and takeoff angles of sediment particles.

Results of Saltation Sediment Particle Movement
The predicted results of the empirical formula for the height, length, and velocity of the saltation sediment proposed by some researchers based on the experimental data are quite scattered [14,62]. The reasons may be the sediment shape, size and density, the bed roughness, and the flow intensity. The other reason may have resulted from the different modes of sediment particle movement. The model results under the six different flow intensities are analyzed to investigate the statistical characteristics of sediment particle movement, including the saltation height and length, movement velocity, the incidence

Results of Saltation Sediment Particle Movement
The predicted results of the empirical formula for the height, length, and velocity of the saltation sediment proposed by some researchers based on the experimental data are quite scattered [14,62]. The reasons may be the sediment shape, size and density, the bed roughness, and the flow intensity. The other reason may have resulted from the different modes of sediment particle movement. The model results under the six different flow intensities are analyzed to investigate the statistical characteristics of sediment particle movement, including the saltation height and length, movement velocity, the incidence, and takeoff angles. In order to make a better comparison between the present simulated results and experimental results, some experimental data for the saltation characteristics of the sediment particle were obtained from the previous study. A brief description of the configuration and parameters of these experiments is presented in Table 2. In this study, the saltation height and length, and incidence and takeoff angle of the saltation sediment are studied and compared with the experimental results. Figure 8 shows the probability density distributions of dimensionless saltation length and height for SP. It is indicated that the dimensionless saltation length and height both follow the Γ function distribution, which is similar to the distribution pattern of the experimental results [15]. The distribution pattern is mainly related to the velocity fluctuations of the flow and the bed roughness. Lu and Willmarth [63] pointed out sediment particle saltation mainly occurs during sweeps (fourth quadrant) and outward interactions (first quadrant), while saltation is almost negligible during ejections (second quadrant) and inward interactions (third quadrant) in the quadrant events under the action of turbulent structures. Furthermore, the collisions between saltation particles increase the saltation length and height, as can be seen in Figure 9. Similarly, the collisions between saltation particles increase the probability of a larger flight time and cause the flight time following the logarithmic normal distribution, as can be seen in Figure 10.

REVIEW 12 of 22
and cause the flight time following the logarithmic normal distribution, as can be seen in Figure 10.      A comparison of the dimensionless saltation length between the present simulated results and experimental results is plotted in Figure 11. It indicates that dimensionless saltation length increases as the transport stage increases. At the lower transport stage, the simulated results of dimensionless saltation length are very close to the experimental results [9,17], which verify the correctness of the numerical model. At the medium transport stage, the simulated results are larger than the experimental results. Einstein and El-Sammi [64] pointed out that saltation length is affected by particle size and shape; their opinions are verified by the experiment of Lee et al. [18] at the low flow intensity, in     A comparison of the dimensionless saltation length between the pre results and experimental results is plotted in Figure 11. It indicates that saltation length increases as the transport stage increases. At the lower t the simulated results of dimensionless saltation length are very close to th results [9,17], which verify the correctness of the numerical model. A transport stage, the simulated results are larger than the experimental re and El-Sammi [64] pointed out that saltation length is affected by particle s their opinions are verified by the experiment of Lee et al. [18] at the low flo A comparison of the dimensionless saltation length between the present simulated results and experimental results is plotted in Figure 11. It indicates that dimensionless saltation length increases as the transport stage increases. At the lower transport stage, the simulated results of dimensionless saltation length are very close to the experimental results [9,17], which verify the correctness of the numerical model. At the medium transport stage, the simulated results are larger than the experimental results. Einstein and El-Sammi [64] pointed out that saltation length is affected by particle size and shape; their opinions are verified by the experiment of Lee et al. [18] at the low flow intensity, in which they pointed out that at the low flow intensity, the dimensionless saltation length and height increase with the particle diameter. Amir et al. [65] believed that the maximum saltation height increases with the increase of the flow depth and particle density. At the medium transport stage, the simulated results are larger than the experimental results. The reason may be that the moving sediment particles collide with each other under the turbulent flow structure, and the inter-particle collisions increase the saltation length, as can be seen in Figure 9. Lee et al. [66] investigated the inter-particle collision behaviors during the saltating process of multiple particles, and found that the maximum saltation height was about 2-2.5 times the values calculated by a single-particle saltation. It is also verified that the trajectories of saltation particles due to inter-particle collision has multiple peak positions, which are clearly different from the typical trajectory of single-particle saltation.
can be seen in Figure 9. Lee et al. [66] investigated the inter-particle collision be during the saltating process of multiple particles, and found that the maximum s height was about 2-2.5 times the values calculated by a single-particle saltation. verified that the trajectories of saltation particles due to inter-particle collision ha ple peak positions, which are clearly different from the typical trajectory of single saltation. Figure 11. Dimensionless saltation length with transport stage [8,9,12,18]. Figure 12 presents a comparison of the dimensionless saltation height betw present simulated results and the experimental results. It is found that the dimen saltation height increases with the increase of the transport stage. At the low tr stage, the dimensionless saltation height is in good agreement with the experim sults by Niño and García [9,17]. Similar to the saltation length, the difference of s height between the present simulated results and experimental results at the m transport stage may be induced by the inter-particle collision, particle size, a depth. Moreover, Auel et al. [61] indicated that for a given transport stage, the s height for a mild bed slope is larger than that for a steep slope, and the difference with the increase of the transport stage.  Figure 11. Dimensionless saltation length with transport stage [8,9,12,18]. Figure 12 presents a comparison of the dimensionless saltation height between the present simulated results and the experimental results. It is found that the dimensionless saltation height increases with the increase of the transport stage. At the low transport stage, the dimensionless saltation height is in good agreement with the experimental results by Niño and García [9,17]. Similar to the saltation length, the difference of saltation height between the present simulated results and experimental results at the medium transport stage may be induced by the inter-particle collision, particle size, and flow depth. Moreover, Auel et al. [61] indicated that for a given transport stage, the saltation height for a mild bed slope is larger than that for a steep slope, and the differences reduce with the increase of the transport stage. medium transport stage, the simulated results are larger than the experimental results. The reason may be that the moving sediment particles collide with each other under the turbulent flow structure, and the inter-particle collisions increase the saltation length, as can be seen in Figure 9. Lee et al. [66] investigated the inter-particle collision behaviors during the saltating process of multiple particles, and found that the maximum saltation height was about 2-2.5 times the values calculated by a single-particle saltation. It is also verified that the trajectories of saltation particles due to inter-particle collision has multiple peak positions, which are clearly different from the typical trajectory of single-particle saltation. Figure 11. Dimensionless saltation length with transport stage [8,9,12,18]. Figure 12 presents a comparison of the dimensionless saltation height between the present simulated results and the experimental results. It is found that the dimensionless saltation height increases with the increase of the transport stage. At the low transport stage, the dimensionless saltation height is in good agreement with the experimental results by Niño and García [9,17]. Similar to the saltation length, the difference of saltation height between the present simulated results and experimental results at the medium transport stage may be induced by the inter-particle collision, particle size, and flow depth. Moreover, Auel et al. [61] indicated that for a given transport stage, the saltation height for a mild bed slope is larger than that for a steep slope, and the differences reduce with the increase of the transport stage.  [8,9,12,18].
The sediment particles under the action of turbulent flow collide with the bed, and three important collision angles are used to characterize the saltation process, i.e., incidence angle θ in , takeoff angle θ out , and collision position angle θ b as shown in Figure 13. The angles of incidence θ in and takeoff θ out are defined with respect to a line parallel to the channel bottom, and they are calculated through the particle trajectory right before and after the collision [33]. The collision position angle θ b is the angle between the tangential collision surface and the bed bottom. Based on the magnitude of the collision position angle, the collision point is located on the upstream (θ b > 0) and downstream faces (θ b < 0) of the bed particle. In this study, the analysis on the incidence angle θ in and takeoff angle θ out are presented. The results on the collision position angle θ b involve the analysis of the "splash function", which will be studied in future research work.
angle between the tangential collision surface and the bed bottom. Based on the magnitude of the collision position angle, the collision point is located on the upstream ( b θ > 0) and downstream faces ( b θ < 0) of the bed particle. In this study, the analysis on the incidence angle in θ and takeoff angle out θ are presented. The results on the collision position angle b θ involve the analysis of the "splash function", which will be studied in future research work.  Figure 14 shows the frequency of distributions of the incidence angle in θ and takeoff angle out θ ; the distribution is also similar to the experimental results [18]. The distributions of the incidence and takeoff angles can be approximately described using the skew-normal distribution function. The distribution pattern is closely related to the random collision of sediment particle-bed and the flow turbulent fluctuations. Bhattacharyya et al. [67] pointed out that the shape of the normal distribution depends on the particle size and bed roughness. In general, the particle's takeoff angle at collision is larger than the particle's incidence angle. This observation is largely supported by the experimental data [12,18,61]. This phenomenon is related to the asymmetry of the forces exerted on the sediment particles. The saltation process of sediment particles is due to the combined interaction of lift and drag forces and the submerged weight. In the rising phase of the SP trajectory, the hydrodynamic lift L F and the submerged weight G of the sediment particle are directed downwards; in the recession phase, the hydrodynamic drag L F is directed upwards and opposes the submerged weight G.  Figure 14 shows the frequency of distributions of the incidence angle θ in and takeoff angle θ out ; the distribution is also similar to the experimental results [18]. The distributions of the incidence and takeoff angles can be approximately described using the skew-normal distribution function. The distribution pattern is closely related to the random collision of sediment particle-bed and the flow turbulent fluctuations. Bhattacharyya et al. [67] pointed out that the shape of the normal distribution depends on the particle size and bed roughness. In general, the particle's takeoff angle at collision is larger than the particle's incidence angle. This observation is largely supported by the experimental data [12,18,61]. This phenomenon is related to the asymmetry of the forces exerted on the sediment particles. The saltation process of sediment particles is due to the combined interaction of lift and drag forces and the submerged weight. In the rising phase of the SP trajectory, the hydrodynamic lift F L and the submerged weight G of the sediment particle are directed downwards; in the recession phase, the hydrodynamic drag F L is directed upwards and opposes the submerged weight G. The particle's incidence angles θin and takeoff angles θout change due to the combined interaction of lift and drag forces and the submerged weight under different flow intensities. Figure 15 presents the variation of incidence angles and takeoff angles with the transport stage. It indicates that the incidence and takeoff angles decrease with the increase of flow intensity. This is related to the relative magnitude of the drag and lift forces (Figure 16c). It also indicates that the present simulated results of incidence angles agree with the experimental results, but the takeoff angles are larger than the experimental results [9,18]. The reason may be related to the arrangement of bed and water depth. The particle's incidence angles θ in and takeoff angles θ out change due to the combined interaction of lift and drag forces and the submerged weight under different flow intensities. Figure 15 presents the variation of incidence angles and takeoff angles with the transport stage. It indicates that the incidence and takeoff angles decrease with the increase of flow intensity. This is related to the relative magnitude of the drag and lift forces (Figure 16c). It also indicates that the present simulated results of incidence angles agree with the experimental results, but the takeoff angles are larger than the experimental results [9,18]. The reason may be related to the arrangement of bed and water depth. transport stage. It indicates that the incidence and takeoff angles de crease of flow intensity. This is related to the relative magnitude of the (Figure 16c). It also indicates that the present simulated results of inci with the experimental results, but the takeoff angles are larger than th sults [9,18]. The reason may be related to the arrangement of bed and The drag and lift forces play significant roles in the sediment mo ence the motion modes of sliding or rolling, saltation, and suspension the variation and relative magnitude of the dimensionless drag and li tical direction above the bed. It indicates that the drag forces are large bed and increase with the increase of the flow intensity, while the lift fo significantly along the flow depth. From the relative magnitude of the the lift force plays a more important role than the drag force in the sedi  (Niño and Garcia,1998) θ out θ in (Lee et al., 2000) θ out θ in (Present simulation) Degree θ/θ c Figure 15. Incidence angles θ in and takeoff angles θ out with transport stage [9,18].

Discussion of Sliding or Rolling and Saltation Sediment Particle Movement
Previous studies mainly focused on the sediment saltation (Shim and Duan, 2017 [15]); therefore, the motion characteristics of different movement modes of bed load particles are investigated in this study. Figure 17 presents the frequency distributions of the dimensionless velocity of SRP and SP. Through curve fitting of the streamwise velocity distributions of sediment particles, skew-normal distributions are found for SRP and SP, respectively, which are similar to the distribution patterns of the experimental results [18] and simulated results [33]. The reason for this skew-normal distribution is related to the disturbance of flow velocity near the bed. The fluid resistance in the flow direction and bed friction resistance induced the asymmetric of particle velocity distribution. The difference between two frequency distributions of the dimensionless velocity is that the mean value of the skew-normal distribution of the SRP is smaller than that of the SP. This is related to the fact that the saltating particles achieved higher and longer trajectories, which induced in the SP obtained more speed and energy from the flow than the SRP. Furthermore, bed resistance has greater influence on the SRP than the SP. The drag and lift forces play significant roles in the sediment movement, and influence the motion modes of sliding or rolling, saltation, and suspension. Figure 16 presents the variation and relative magnitude of the dimensionless drag and lift forces in the vertical direction above the bed. It indicates that the drag forces are large near the sediment bed and increase with the increase of the flow intensity, while the lift forces do not change significantly along the flow depth. From the relative magnitude of the drag and lift forces, the lift force plays a more important role than the drag force in the sediment saltation [18].

Discussion of Sliding or Rolling and Saltation Sediment Particle Movement
Previous studies mainly focused on the sediment saltation (Shim and Duan, 2017 [15]); therefore, the motion characteristics of different movement modes of bed load particles are investigated in this study. Figure 17 presents the frequency distributions of the dimensionless velocity of SRP and SP. Through curve fitting of the streamwise velocity distributions of sediment particles, skew-normal distributions are found for SRP and SP, respectively, which are similar to the distribution patterns of the experimental results [18] and simulated results [33]. The reason for this skew-normal distribution is related to the disturbance of flow velocity near the bed. The fluid resistance in the flow direction and bed friction resistance induced the asymmetric of particle velocity distribution. The difference between two frequency distributions of the dimensionless velocity is that the mean value of the skew-normal distribution of the SRP is smaller than that of the SP. This is related to the fact that the saltating particles achieved higher and longer trajectories, which induced in the SP obtained more speed and energy from the flow than the SRP. Furthermore, bed resistance has greater influence on the SRP than the SP. difference between two frequency distributions of the dimen the mean value of the skew-normal distribution of the SRP is SP. This is related to the fact that the saltating particles achie trajectories, which induced in the SP obtained more speed an than the SRP. Furthermore, bed resistance has greater influen SP. The formulas of mean velocity of a saltating particle are pr investigators based on experimental data, as can be seen in Tabl comparison of the dimensionless streamwise saltation velocity b ulated results and the empirical formula [7,8,11,13]. It is found th of the empirical formula are quite scattered, except for the empir study, the present simulated results of the SRP and SP are prese the distinction of different movement modes. Figure 18 indicate lated results of the SRP are lower than the predicted results of the The formulas of mean velocity of a saltating particle are proposed by the different investigators based on experimental data, as can be seen in Table 3. Figure 18 presents a comparison of the dimensionless streamwise saltation velocity between the present simulated results and the empirical formula [7,8,11,13]. It is found that the predicted results of the empirical formula are quite scattered, except for the empirical formula [7,8]. In this study, the present simulated results of the SRP and SP are presented in order to analyze the distinction of different movement modes. Figure 18 indicates that the present simulated results of the SRP are lower than the predicted results of the empirical formula, and the present simulated results of the SP are between the curves of the empirical formula. For the SRP, the sediment particles move near the sediment bed and contact with the bed frequently, so the bed resistance decreases the streamwise velocity of the SRP. For the SP, the collisions between sediment particles under the strong turbulence structure increase the flight time, and cause saltation particles to absorb more energy from the flow and increase sediment particle velocity. The other possible reason for the difference between the present simulated results and the predicted results of the empirical formula is the shape and size of particles. The smooth particles move with the lower and shorter jumps more than the rough particles of the same mean diameter, and are less influenced by the flow velocity than the coarse particles; therefore, the velocities of smooth particles are less than those of the rough particles [68]. Table 3. Formulas of mean velocity u b of saltating particles given by different investigators.
It is well known that the sediment particles move in both translational and rotational motion, and the rotation of sediment particles has a significant effect on their trajectories [8,9,25]. Sediment particle rotation occurs under the combined effects of the bed roughness, particle-bed collision, and the velocity gradient of the carrier fluid. In the previous experiment, the particle's angular velocity was hard to describe exactly due to the limits of the measurement technique. The dimensionless angular velocity of particle Ω is defined as * / sd u Ω = (13) The frequency distributions of the instantaneous dimensionless angular velocity of the SRP and SP around z direction at the different vertical positions is presented in Figure 19. Through curve fitting of the frequency distributions of instantaneous dimensionless angular velocity, it shows that, for the SRP, the frequency distribution of dimensionless angular velocity follows skew-normal distribution, while for the SP, the frequency distribution of dimensionless angular velocity follows a logarithmic normal distribution. The different distribution patterns are mainly related to the motion modes of the sediment particle. Particle rotation is mainly derived from sediment  Figure 18. Comparison of dimensionless streamwise velocity of moving particle between the present model results and the empirical formula with Shields number [7,8,11,13].
It is well known that the sediment particles move in both translational and rotational motion, and the rotation of sediment particles has a significant effect on their trajectories [8,9,25]. Sediment particle rotation occurs under the combined effects of the bed roughness, particle-bed collision, and the velocity gradient of the carrier fluid. In the previous experiment, the particle's angular velocity was hard to describe exactly due to the limits of the measurement technique. The dimensionless angular velocity of particle Ω is defined as The frequency distributions of the instantaneous dimensionless angular velocity of the SRP and SP around z direction at the different vertical positions is presented in Figure 19. Through curve fitting of the frequency distributions of instantaneous dimensionless angular velocity, it shows that, for the SRP, the frequency distribution of dimensionless angular velocity follows skew-normal distribution, while for the SP, the frequency distribution of dimensionless angular velocity follows a logarithmic normal distribution. The different distribution patterns are mainly related to the motion modes of the sediment particle. Particle rotation is mainly derived from sediment particle-bed collisions. The SRP slide and roll close to the bed and frequently contact with the bed, while the SP are intermittently losing contact with the bed for a short time and are strongly influenced by the fluid velocity near the bed.  Figure 20 presents the variation of particle angular velocity with the transport stage. It indicates that the angular velocity of the SP decreases with the increase of the transport stage; this trend is similar to the experimental results of Niño and García [9]. The difference may result from the particle shape and bed roughness. The shape of the particle affects the angular velocity of the particle, because a particle with an elongated shape colliding with the bed results in a higher angular velocity than a spherical shape [68]. Niño and García [9] pointed out that the particle angular velocity is controlled mainly by the process of collision with the bed. This phenomenon can be seen in Figure 7. The incidence angle and takeoff angle at collision also influence the particle angular velocity [8]. It also indicates that the angular velocity of the SRP increases with the increase of the transport Figure 20 presents the variation of particle angular velocity with the transport stage. It indicates that the angular velocity of the SP decreases with the increase of the transport stage; this trend is similar to the experimental results of Niño and García [9]. The difference may result from the particle shape and bed roughness. The shape of the particle affects the angular velocity of the particle, because a particle with an elongated shape colliding with the bed results in a higher angular velocity than a spherical shape [68]. Niño and García [9] pointed out that the particle angular velocity is controlled mainly by the process of collision with the bed. This phenomenon can be seen in Figure 7. The incidence angle and takeoff angle at collision also influence the particle angular velocity [8]. It also indicates that the angular velocity of the SRP increases with the increase of the transport stage. The greater the transport stage, the faster the rolling sediment particles move, and accordingly, the greater the angular velocity of the SRP.  Figure 20 presents the variation of particle angular velocity with the transport stage. It indicates that the angular velocity of the SP decreases with the increase of the transport stage; this trend is similar to the experimental results of Niño and García [9]. The difference may result from the particle shape and bed roughness. The shape of the particle affects the angular velocity of the particle, because a particle with an elongated shape colliding with the bed results in a higher angular velocity than a spherical shape [68]. Niño and García [9] pointed out that the particle angular velocity is controlled mainly by the process of collision with the bed. This phenomenon can be seen in Figure 7. The incidence angle and takeoff angle at collision also influence the particle angular velocity [8]. It also indicates that the angular velocity of the SRP increases with the increase of the transport stage. The greater the transport stage, the faster the rolling sediment particles move, and accordingly, the greater the angular velocity of the SRP.

Conclusions
In this study, the motion trajectories, the dynamic responses, the corresponding hydrodynamic forces, and the statistical characteristics of sediment particle movement (sliding, rolling, and saltation) are studied through simulating the movement of sediment particles using a coupled CFD-FDEM model. The conclusions are summarized as follows: The bed load sediment particles have a different motion trajectory under the fluctuating velocity of turbulent flow. The inter-particle collision between moving sediment particles has an important effect on the characteristics of moving sediment particles. The inter-particle collision in the low concentration increases the saltation length, height, and

Conclusions
In this study, the motion trajectories, the dynamic responses, the corresponding hydrodynamic forces, and the statistical characteristics of sediment particle movement (sliding, rolling, and saltation) are studied through simulating the movement of sediment particles using a coupled CFD-FDEM model. The conclusions are summarized as follows: The bed load sediment particles have a different motion trajectory under the fluctuating velocity of turbulent flow. The inter-particle collision between moving sediment particles has an important effect on the characteristics of moving sediment particles. The inter-particle collision in the low concentration increases the saltation length, height, and flight time, correspondingly increasing the velocity of the saltating sediment particle. The frequency distribution of dimensionless angular velocity for the sliding and rolling sediment particle follows skew-normal distribution, while for the saltation sediment particle it follows a logarithmic normal distribution. Therefore, the inter-particle collision effect between particles should be considered in the calculation of sediment transport rate.
The collision between moving sediment particles and the bed has an important effect on the motion pattern and characteristics of moving sediment particles. The hydrodynamic forces exerted on a sediment particle away from the bed are smooth, while near the bed they present a relatively obvious disturbance owing to the contact between the sediment particle and the bed. The drag force near the bed is large, and away from the bed it decreases gradually during the whole movement process. The relative magnitudes of drag and lift forces influence the incidence and takeoff angles of saltation sediment particles. The particle's angular velocity is the largest at the takeoff moment, then decreases gradually owing to the viscous effects. The research results are helpful to further understand the distinction of motion characteristics of sliding, rolling, and saltating sediment particles, and to better predict the bed load transport rates.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the first author.

Acknowledgments:
The author would like to thank Chunning Ji, Tianjin University, and Baosheng Wu, Tsinghua University for their assistances in the numerical simulation. All the computational tasks were performed on the Explore 100 cluster system of Tsinghua National Laboratory for Information Science and Technology.

Conflicts of Interest:
The authors declare no conflict of interest. Young's modulus (kg m s −2 ) F fluid-particle interaction force (kg m s −2 ) F n normal force between particles (kg m s −2 ) F s tangential force between particles (kg m s −2 ) F D drag force of hydrodynamic forces exerted on particle (kg m s −2 ) F L lift force of hydrodynamic forces exerted on particle (kg m s −2 ) G gravity acceleration (m s −2 ) G particle submerged weight (kg m s −2 ) Ga Galileo number (-) H convective and diffusive terms (-) h e smallest size of the finite elements (m) h f net water depth (m) I interpolation function (-) k p penalty parameter (-) k s tangential spring stiffness (-) k s equivalent sand roughness (-) N ibp total number of immersed boundary points (-) Re b Reynolds number based on bulk flow velocity (-) Re t Reynolds number based on friction velocity (-) Re p particle Reynolds number (-) S angular velocity of the particle (Hz) SP saltation sediment particle (-) SRP rolling or sliding particle (-) T matrix transposition (-) u b bulk velocity (m s −1 ) u p streamwise velocity of the particle (m s −1 ) u * friction velocity (m s −1 ) U fluid velocity (m s −1 ) v p vertical velocity of the particle (m s −1 ) V volume of the tetrahedral finite element (m 3 ) V velocity of the immersed boundary points (m s −1 ) V i volume of the corresponding sub-tetrahedrons (m 3 ) V s tangential relative velocity between particles (m s −1 ) ∆V i discrete volume around the immersed boundary points (m 3 ) ∇ gradient operator (-) Γ overlapping volume between finite elements of the contactor and target discrete elements (m 3 ) Ω dimensionless angular velocity of particle (-)