Single sediment dynamics in turbulent flow over a porous bed – insights from interface-resolved simulations

We use interface-resolved direct numerical simulations to study the dynamics of a single sediment particle in a turbulent open channel flow over a fixed porous bed. The relative strength of the gravitational acceleration, quantified by the Galileo number, is varied so as to reproduce the different modes of sediment transport – resuspension, saltation and rolling. The results show that the sediment dynamics at lower Galileo numbers (i.e. resuspension and saltation) are mainly governed by the mean flow. Here, the regime of motion can be predicted by the ratio between the gravity and the shear-induced boundary force. In these cases, the sediment particle rapidly takes off when exposed to the flow, and proceeds with an oscillatory motion. Increasing the Galileo number, the frequency of these oscillations increases and their amplitude decreases, until the transport mode switches from resuspension to saltation. In this case, the sediment travels by short successive collisions with the bed. Further increasing the Galileo number, the particle rolls without detaching from the bed. Differently from the previous modes, the motion is triggered by extreme turbulent events, and the particle response depends on the specific initial conditions, at fixed Reynolds number. The results reveal that close to the onset of sediment motion, only turbulent sweeps can effectively trigger the particle motion by increasing the stagnation pressure upstream. We show that for the parameters in this study, a criterion based on the streamwise flow-induced force can successfully predict the incipient movement.


A24-2
A. Yousefi, P. Costa and L. Brandt and many engineering and industrial processes. A fundamental understanding of these transport mechanisms is therefore relevant to physicists and engineers, with a key role often played by the identification of the threshold conditions for the onset of sediment motion, also known as incipient motion. Despite several attempts to address this problem, it has remained elusive owing to the large number of parameters determining the particle dynamics. In particular, the complexity rises from the fact that the particle entrainment is highly sensitive to the bed surface morphology, particle size and shape, exposure and packing conditions, and, not least, to the intermittent behaviour of turbulence in near-bed flows (Dey & Ali 2018).
In general, sediments moving through a fluid flow experience two types of transport: bed load and suspended load. In bed load transport, the sediments slide, roll and hop (or saltate) along the bed; in suspended load transport, sediment motion away from the bed is sustained by the bulk flow (Bagnold 1966;Chanson 2004). A first approach to the fundamental understanding of this problem focuses on the simplified problem of the flow-induced dynamics of a single sediment particle, initially resting on a sediment bed.
Among the deterministic approaches proposed, the most widely used parameter to predict the onset of incipient motion originates from the seminal work of Shields (1936). There, it is argued that the incipient motion occurs once a threshold nondimensional boundary shear stress is exceeded. This non-dimensional number is the socalled Shields parameter Sh ≡ τ 0 /((ρ p − ρ)gD), where τ 0 is the boundary shear stress, ρ p and ρ are the mass density of the sediments and of the carrying fluid, g is the gravitational acceleration and D the sediment particle diameter. The Shields parameter is based on the time-averaged characteristics of the fluid-solid interaction forces, and is therefore not capable of accounting for the fluctuating forces typically encountered in natural and engineering turbulent flows . The incipient motion of individual sediment particles may be observed when the deterministic bedload predictor is below the threshold of motion (Kleinhans & van Rijn 2002), often identified empirically. Although many have tried to improve the initial idea behind the Shields stress criterion (see e.g. Ling (1995), Zanke (2003), Vollmer & Kleinhans (2007) and Ali & Dey (2016)), a distinct threshold has not been established due to the complex interactions between the fluid and the particle bed, a function also of the several parameters characterising the problem.
To overcome these issues, stochastic approaches have been suggested as a viable alternative to estimate the probability of events that exceed the threshold of incipient motion. Cheng & Chiew (1998) presented a theoretical derivation of the pickup probability of sediment entrainment on a hydraulically rough flow by considering the instantaneous streamwise velocity to follow a Gaussian distribution. Wu & Chou (2003) extended this model by proposing a theoretical formulation for the rolling and lifting probabilities for sediment entrainment in hydraulically smooth and transitional flows. This model accounts for the relation between lift coefficient and particle Reynolds number, and criteria for the modes of particle transport in contact and detachment are obtained by applying the equilibrium of forces and moments. For a given relative particle size, the probability of transport in contact mode increases with increasing Shields stress, until the transition from contact to detachment mode occurs, where it starts to decrease. On the other hand, the probability of transport in detachment mode increases monotonically with an increase of the Shields stress (Dey & Ali 2017).
Another modelling approach uses the concept of force impulse to deal with the bed particle entrainment (see Diplas et al. 2008). In this approach, the conditions for Single sediment dynamics in turbulent flow over a porous bed 893 A24-3 incipient motion are determined not only by the magnitude of the maximum force acting on individual grains, but also by the time interval a force exceeding a certain threshold is applied, all for a given bed material and configuration (Celik, Diplas & Dancey 2014). Celik et al. (2010) emphasise that the turbulence-induced impulse must exceed the threshold impulse for the motion to start. In addition to the concept of an impulse-based criterion, energy-based approaches have been suggested to tackle the problem. In this case, the geometrical pocket formed by the bed particles is viewed as a potential well with the gravitational and frictional mechanisms imposing an energy barrier for the sediment particle to fully escape the pocket (see e.g. Lee, Ha & Balachandar (2012) and Valyrakis, Diplas & Dancey (2013)).
Other studies have focused on a more phenomenological picture. The idea is that coherent structures, moving from the turbulent boundary layer towards the sediment bed, disrupt the viscous sublayer and impinge directly onto the first layer of grains. The interaction between these structures and the individual sediments results in their entrainment. This hypothesis was first introduced by Sutherland (1967), who attributed the (re)suspension mechanism to turbulent bursts. Cleaver & Yates (1973) suggested a two-stage mechanism of resuspension due to turbulent bursting: first a relatively slow movement away from the surface in the viscous sublayer, followed, by a rapid movement into the bulk due to a second burst. The main events considered responsible for incipient motion of sediment particles are related to large negative values of the Reynolds shear stresses: the upward motion of slowly moving-fluid parcels (ejections) and downward inrush of rapidly moving fluid parcels (sweeps). Séchet & Le Guennec (1999) argue that the dominant mechanism for particle transport is by ejections, as the time and space-scale characteristics of the particle motion in the period between two jumps are commensurable with the corresponding ejections. Moreover, Ninto & Garcia (1996) measured the instantaneous particle velocities during ejections in both smooth and transitionally rough flows, and showed that the streamwise velocity component tends to be much lower than the local mean flow velocity, while their vertical component tends to be rather intense. These fluctuations are much larger than the local standard deviation of the vertical flow velocity, which would indicate that such particles are responding to rather extreme ejection events. On the other hand, Dey, Sarkar & Solari (2011) suggested that sweep events with downward advection of the streamwise Reynolds normal stress are prevalent during sediment entrainment. The experiments of Dwivedi, Melville & Shamseldin (2010) showed the presence of high vertical and horizontal pressure gradients on the sediment particle during the occurrence of a sweep event. The fact that the contribution from sweep events to the Reynolds stresses near porous media increases (Suga, Mori & Kaneda 2011), strengthens this viewpoint.
The continuous growth in computing power, together with the development of efficient numerical methods, made it possible to study fluid-sediment interactions in turbulent flows. Studies where the ratio of particle size to smallest turbulent structures is small highlight the role of ejections on detachment, or even suspension of these particles (see e.g. Soldati &Marchioli (2009) andVinkovic et al. (2011)). For larger particle Reynolds number, interface-resolved approaches are required. Vowinckel et al. (2016) used interface-resolved simulations to study sediment transport through a loosely packed fixed bed; the authors highlight the erosion mechanism, triggered by collision with already eroded particles (see also Frey & Church (2011)). The study of the collective behaviour of sediment particles and pattern formation in a sediment bed has also become accessible through fully resolved numerical simulations (see e.g. Ji et al. (2014), Kidanemariam and Uhlmann (2014) (2019)).
In this work we investigate the dynamics of a finite-size sediment particle (with diameter D ≈ 20ν/u τ , where u τ is the bed friction velocity and ν the fluid viscosity) on a porous bed. We fix the size ratio of particle to smallest turbulent structures, and the ratio of particle inertia to fluid inertia, and vary the relative strength of fluid force to buoyant force. Parameters are chosen to include three distinct modes of sediment transport: suspension, saltation and rolling. We aim to investigate the dynamics of the sediment particle for each case by performing fully resolved numerical simulations, which allow us to explore the flow physics with the three-dimensional and time-resolved insight that is still difficult to achieve in a real experiment.
Our results show that for relatively high values of the Shields parameter, Sh, the dynamics are governed by the balance between fluid forces induced by the mean flow and buoyancy. Conversely, for lower values of the Sh -well below the corresponding critical value for the onset of sediment motion -we observe that the particle motion is solely triggered by extreme, strong sweep events. The associated pressure distribution close to the particle shows a clear and consistent fingerprint: an upstream region of higher-than-mean stagnation pressure (a consequence of an increased streamwise velocity), together with a smaller high-pressure spot close to the contact point of the sediment particle with the bed.
This manuscript is organised as follows. In § 2, we present the governing equations and numerical methodology, together with the flow configuration. The role of the Galileo number and the underlying mechanisms of detachment are discussed in § 3. Finally, a summary of the main findings and conclusions are drawn in § 4.

Governing equations
The numerical algorithm solves the Navier-Stokes equations for an incompressible Newtonian fluid with density ρ and kinematic viscosity ν, where u is the fluid velocity vector, p the modified fluid pressure (i.e. relative to the local hydrostatic load) and ∇p e a constant pressure gradient that drives the flow. The last term on the right-hand side of (2.2) acts close to the surface of the particles to model their presence and is due to the direct forcing immersed boundary method (IBM) used for the fluid-solid coupling. The dynamics of the rigid, spherical particles are governed by the Newton-Euler equations for conservation of linear and angular momentum: with u p and ω p the particle linear and angular velocity. Here, m p , I p , ρ p and V p denote the particle mass, moment-of-inertia, mass density and volume, r is the position vector with respect to the particle centre and n the outward-pointing normal to the particle Single sediment dynamics in turbulent flow over a porous bed 893 A24-5 surface ∂Ω p . The fluid stress tensor is given by τ = −pI + νρ(∇u + ∇u T ). Finally, g denotes the gravitational acceleration, and F col and T col the force and torque resulting from short-range particle-particle and particle-wall interactions, such as lubrication and collisions. The sets of equations governing each phase are coupled through a no-slip and nopenetration condition at the particle surface, i.e. u| ∂Ω p = u p + ω p × r.
(2.5) This is achieved by extending a standard Navier-Stokes solver with a direct-forcing IBM, as described below.
2.2. Numerical method We use the direct forcing IBM developed by Uhlmann (2005) and modified by Breugem (2012), briefly described hereafter. The Navier-Stokes equations governing the fluid phase are solved on a uniform ( x = y = z), staggered, Cartesian grid. The spherical particles are discretised by a set of Lagrangian points, uniformly distributed along their surface. The IBM forcing scheme requires three steps: (i) the fluid prediction velocity is interpolated from the Eulerian to the Lagrangian grid, (ii) the IBM force required for matching the local fluid velocity and the local particle velocity is computed on each Lagrangian grid point and (iii) the resulting IBM force is spread from the Lagrangian to the Eulerian grid. The interpolation and spreading operations are done through the regularised Dirac delta function of Roma, Peskin & Berger (1999), which extends over three grid cells in all coordinate directions. Regularisation of the particle-fluid interface can result in a loss of accuracy. Breugem (2012) showed that by slightly retracting the Lagrangian grid, second-order accuracy for particle-averaged quantities (e.g. drag) can be attained. Due to the kernel support, the forces obtained from different neighbouring Lagrangian points can be spread to the same Eulerian point, which reduces the accuracy of the forcing. This issue is circumvented by iterating the forcing scheme as proposed by Luo et al. (2007). This multi-direct forcing scheme is particularly relevant in the present set-up, because the kernel's spreading range for forcing points pertaining to two nearly touching spheres may also partially overlap.
When the distance between two particles or a particle and a wall is smaller than x, the IBM fails to resolve the short-range hydrodynamic interactions. The particle dynamics are therefore solved with a subgrid-scale model for short-range interactions, based on canonical lubrication interactions in the Stokes regime, and a soft-sphere model for solid-solid contacts; see Costa et al. (2015). This model is described below for the sake of clarity.
Unresolved normal lubrication interactions and roughness effects are accounted for by a two-parameter model. For a normalised gap distance ε x (with ε ≡ δ ij,n /R p , with δ ij,n being the gap distance between particles i and j), a lubrication correction force is computed. Here ε x corresponds to the threshold below which the IBM fails to resolve canonical lubrication interactions in the Stokes regime. The closure uses a correction force F lub = −6πµR p u ij,n (λ(ε) − λ(ε x )), added to the integration of the particle Newton equation (2.3), where λ is the Stokes amplification factor. Roughness effects are accounted for by making this correction independent of the normalised gap width ε when this quantity is smaller than a threshold ε σ , i.e. F lub ∝ u ij,n (λ(ε σ ) − λ(ε x )), for 0 < ε < ε σ . This correction is applied until the surfaces undergo actual solid-solid contact; then the soft-sphere collision model takes over.

A24-6
A. Yousefi, P. Costa and L. Brandt The soft-sphere collision model computes the normal component of the collision force from the following linear spring-dashpot model: where δ ij and u ij,n are the gap width and relative particle velocity projected in the line of centres n ij ≡ (x j − x i )/( x j − x i ) where x i/j corresponds to the centroid position of particles i/j. The spring and dashpot coefficients are given by k n = m e (π 2 + ln 2 e n,d ) where e n,d is the dry coefficient of restitution and m e = (m −1 i + m −1 j ) −1 the reduced mass of the particles. Here, N t is the collision time, set as a multiple N of the time step of the Navier-Stokes solver, t. This allows the flow solution to gradually adapt to the sudden changes in particle velocity (Costa et al. 2015;Biegert, Vowinckel & Meiburg 2017). The model for the tangential component of the collision force uses a similar force formulation, but with a Coulomb friction model for sliding motion Similarly to the formulation of the normal force, the spring and dashpot coefficients are given by with m e,t = (2/7)m e , where e t,d denotes the dry tangential coefficient of restitution and µ c is the Coulomb coefficient of sliding friction. We refer to Costa et al. (2015) for more details.

Computational set-up
We simulate the dynamics of a single mobile particle in a channel flow, placed on top of a porous bed of fixed spheres, of the same size as the mobile particle (figure 1). Direct numerical simulations (DNS) of fully developed turbulent open channel flow are performed in a domain with dimensions of L x × L y × L z = 40D × 13D × 14D, with x, y and z being the streamwise, bed-normal and spanwise directions and D the particle diameter. The fixed bed consists of 2560 spheres in four layers, arranged in a hexagonal close packing configuration. The particles are resolved with D/ x = 35 grid points per particle diameter, which corresponds to an inner-scaled grid spacing of x + = u τ x/ν = 0.6, which resolves all flow scales; here, u τ = √ τ tot /ρ is the friction velocity with τ tot the total stress, i.e. the sum of the viscous and the Reynolds stress evaluated at the fluid-porous interface (Breugem, Boersma & Uittenbogaard 2006). The flow is periodic in the streamwise and spanwise directions, with no-slip/nopenetration and free-slip/no-penetration boundary conditions imposed at the bottom and top boundaries. A uniform pressure gradient (∇p e ) forces the flow such that the Single sediment dynamics in turbulent flow over a porous bed 893 A24-7 Shows the whole domain, the moving particle is depicted in orange; D is the particle diameter, h is the height of fluid section and h b is bed height. (b) Shows a close up of the moving particle and bed configuration. For clarity, only the top layer of the bed is shown. The reference frame adopted here has its origin on the plane tangent to the top layer of particles, and at the streamwise and spanwise position of the centre of the sediment particle.
bulk fluid velocity U b remains constant. Here the bulk Reynolds number is defined as The moving (non-Brownian) particle has a density ratio ρ p /ρ = 1.05 (e.g. polystyrene in water (van Hout 2013)), and an inner-scaled diameter D + = u τ D/ν = 21.5. The initial position of the sediment particle centre with respect to the computational box is [x c , y c , z c ] = [20.5D, (4 The reference frame adopted here is depicted in figure 1. In the present study we vary the Galileo number Ga ≡ (ρ p /ρ − 1)gD 3 /ν, which quantifies the relative importance of gravitational to viscous forces. This is achieved by varying the gravitational acceleration (pointing downward in y). An initial set of simulations was performed to seek the Galileo numbers corresponding to the different modes of sediment transport considered. The particle collisional parameters are chosen to be close to what is often measured in collision measurements of rigid spheres: e n,d = 0.97, e t,d = 0.1 and µ s = 0.15 (see Foerster et al. 1994;Yang & Hunt 2006). The collision time is set to N = 8 time steps of the Navier-Stokes solver, corresponding to a collision time of around 0.017h/U b . Other physical and computational parameters are summarised in table 1.
To obtain the initial turbulent field, the simulations start from laminar Poiseuille flow with a high amplitude disturbance consisting of streamwise counter-rotating vortices to trigger the transition to turbulence efficiently (Henningson & Kim 1991). The position and rotation of the moving particle are fixed at this stage to allow the turbulent field to develop around the particle. When fully developed turbulence is obtained, the mobile particle is released, i.e. the particle dynamics are governed by (2.3) and (2.4). In the text, time has been normalised by the bulk fluid velocity U b and the particle diameter D, i.e. tU b /D denotes the non-dimensional time.

Results
We start by characterising the turbulent flow over the porous bed of spheres. The profiles of averaged Eulerian fluid statistics reported here correspond to mean denotes the Shields parameter, which is proportional to the ratio of fluid force on the particle to the weight of the particle. The critical Shields parameter for this particle Reynolds number is equal to Sh cr = 0.0344 (see Cao, Pender & Meng 2006). Here, N ini denotes the number of initial conditions tested for each case and N inc the number of the realisations in which incipient movement is detected.
intrinsic averages. The intrinsic average of a quantity ξ is computed as where Ψ ijk,t is the fluid volume fraction at the grid cell ijk and instant t, and y j the wall-normal location of the averaging bin, which extends over the entire domain in the two homogeneous directions. The fluid velocity/pressure fluctuations are defined consistently by subtracting the corresponding fluid intrinsic average from the fluid phase velocity. The mean fluid velocity and the bed-normal velocity fluctuations profiles, normalised by u τ , are depicted in figure 2. The bed is weakly permeable in this case, with porosity = 0.28, defined as the volume fraction of voids in the bed, so it behaves similarly to a solid wall as the mean velocity profile is reminiscent of a typical channel-flow profile. Nonetheless, close to the bed, the mean profile is altered by the particle form-induced drag and reaches U = 0 only below the nominal surface level. This is in agreement with the results of Breugem et al. (2006) and Rosti, Brandt & Pinelli (2018), considering the porosity of the bed in our simulations ( = 0.28).
Comparison with the DNS data of Kim et al. (1987) for smooth wall, shows slightly larger peaks in the profiles of the velocity fluctuations in the cross-stream direction, while the peak of the streamwise root mean square (r.m.s.) velocity profile near the bed is hardly changing. Breugem et al. (2006) argue that this behaviour can be attributed to the weakening of the wall-blocking and wall-induced viscous effects; increasing the wall porosity, the peak in the streamwise r.m.s. profile decreases and elongated streaky structures weaken, replaced by almost two-dimensional structures. Note that the difference in the normal component between our results (blue line) and those of Kim et al. (1987) (blue circles) at a normal location close to y/h = 1 is due to the difference in the boundary condition of our simulations for an open (half-)channel flow compared to that of Kim et al. (1987) for a channel flow with a top wall.
As the size of the moving particle is small compared to the size of the domain, its presence and its eventual motion does not affect the mean flow statistics and the The colour indicates the non-dimensional spanwise angular velocity. The red line is fixed on the two-dimensional projection of the particle at beginning of the simulations and changes according to the solid-body rotation of the sediment particle; it therefore depicts the change in the angular position of the particle.
turbulent flow can be considered to be the same for the different cases, i.e. for the different Galileo numbers investigated. The effect of the Galileo number on the dynamical response of the moving particle on top of the bed is visualised first in order to introduce the classification used in the next subsections. We display in figure 3 the trajectory of the moving particle in the bed-normal x-y plane for the different Galileo numbers under investigation.

A24-10
A. Yousefi, P. Costa and L. Brandt All the cases in the figure use the same turbulent flow field, which resulted in incipient movement in all cases, as initial condition. The figure clearly shows that the streamwise displacement of the moving particle varies significantly with the Galileo number. Note also that the maximum off-plane displacement (i.e. in the spanwise direction and computed as average between the different instances for each case; not shown in the figure) is approximately 5D for Ga25, around 1D for Ga150 and less than 0.2D for Ga250 and Ga310. As concerns the bed-normal motion, the particle detaches instantly from the bed as exposed to the turbulent field in the case Ga25. The particle spends most of the time detached from the bed and bounces on the bed two times during our simulation, at x/D ≈ 23 and 28; at this time, it gains significant angular velocity due to the tangential collision force with the bed. The trajectory pertaining to the case Ga25 is very similar to the two-stage lift-off pattern observed in van Hout (2013). The particle with Ga150, see figure 3, experiences a sequence of relatively small excursions and quick returns to the wall. Due to these contacts, the particle has, on average, a higher angular velocity than in the case Ga25, as shown by the colour scale. For cases Ga250 and Ga310, the particle is in contact with the bed for the whole duration of the simulation. Note that the particle angular velocities reveal that in both cases the particles display an intermittent roll/no-roll movement. The duration of the periods of zero velocity is expectedly larger for the case Ga310, as the particle travels a smaller distance, 10.5D, whereas the particle with Ga250 moves 32D during the same time interval.
For a more detailed analysis of each case, we categorise the cases depending on the contact with the bed; in other words, cases Ga25 and Ga150, having no contact and limited contact with the bed, and cases Ga250 and Ga310, when the particles are long in contact with the bed, will be examined together. For each case, the criterion based on the Shields parameter for predicting the initiation of motion is tested (see table 1). One should note that this parameter emphasises the time-averaged characteristics of the applied fluid forces  and the criterion based on this parameter is currently in use for determining the onset of motion when a significant number of sediment particles are subjected to erosion. Despite this, we use this parameter in our study of a single-particle dynamics to put our results in perspective against this classical criterion. Also, we use the terminologies resuspension, saltation, incipient motion and rolling to address the single-particle dynamics when it resembles that of the transported particles in each of the mentioned modes of sediment transport. Indeed, in a natural flow when the Shields parameter is high enough, a significant bed-load motion and suspension would take place (see e.g Revil-Baudard et al. 2015;Berzi & Fraccarollo 2016). However, in our study, as we focus on the single-sediment dynamics, the cases labelled as resuspension and saltation would mimic the movement of a sediment particle, which, in a real scenario, would be accompanied by a significant number of other sediments.

Lower Galileo numbers: resuspension and saltation
We start with cases Ga25 and Ga150, when the particle does not remain in contact with the bed. Each of these cases was run with four different turbulent initial conditions. Figure 4 shows the time history of the bed-normal position of the particle centre for all the realisations of case Ga25. For all cases, we observe that the sediment particle dislodges immediately after its release. This suggests that the onset of sediment motion and the subsequent dynamics can be predicted by the mean flow characteristics. From these realisations, the one denoted 'initial condition 1' is selected for further investigations of case Ga25. . Time history of the bed-normal position of the centre of the sediment particle for Ga25, coloured by the bed-normal fluid force exerted on the particle F y , scaled with (ρu 2 τ A p ); A p = πD 2 /4 is the projected area of the particle onto the x-z plane. The inset shows a detailed view of the initial stage of particle motion. Figure 5 reports the time history of the bed-normal position of the centre of the free sediment for Ga25. The colour shows the y-component of the total fluid surface force exerted on the particle F y = ( ∂Ω p τ · n dA) y , normalised by the boundary shearinduced force scale ρu 2 τ A p , where A p = πD 2 /4 is the area of the particle projected onto the x-z plane. Alternatively, one can scale the fluid force with the particle submerged weight (ρ p − ρ)gV p , where V p = πD 3 /6 is the volume of the sediment particle. Note that the ratio between these scales is proportional to the classical Shields parameter, which here succeeds in predicting the onset of sediment transport. The figure shows that the fluid surface force is of the same order of magnitude as the boundary shearinduced force, in other words, the vertical force can be related to the mean turbulent profile. The ratio between these two scaling factors, the boundary shear-induced force and the particle submerged weight, is 0.96 for the case Ga25, and if we scale the fluid force with the latter, the result will appear as those in figure 5. The particle vertical acceleration is therefore relatively small, and the particle hovers in the near wall region (y/D < 3). Figure 6 depicts the time history of the bed-normal position of the centre of the sediment in case Ga25 together with its streamwise (apparent) slip velocity, defined as the difference between the mean fluid velocity and the particle velocity, evaluated at the particle centre in panel (a). Panel (b) of the same figure depicts the bed-normal acceleration of the particle (black line) together with the mean shear du/dy experienced by the particle (red line), i.e. the bed-normal derivative of the streamwise fluid velocity, averaged in a spherical shell with the diameter 2D Time history of the bed-normal particle acceleration (black) together with the mean shear du/dy experienced by the particle (red) for case Ga25; the mean flow velocity is evaluated at the particle centre, while the shear rate is averaged in a spherical shell with the diameter 2D around the particle centre.
around the particle centre. The observed cyclic dynamics are qualitatively explained by a balance between shear-induced forces and gravity. First, high shear and slip velocity cause a large shear-induced lift force (see e.g. Saffman (1965), Bagchi & Balachandar (2002), Thomas et al. (2015)) that triggers the vertical motion. Second, as the particle moves upward, the local shear rate and apparent slip velocity starts to decrease; then the particle submerged weight starts to dominate and the particle is eventually driven downward. Finally, during its downward movement, the high momentum acquired during the previous ascent is dissipated, which increases the slip velocity. This, together with the increase in local shear rate, increases again the relative importance of the shear-induced lift force. Note that in panel (b), the particle vertical acceleration (i.e. the net force acting on the particle), does not correlate perfectly with the streamwise slip velocity and mean shear rate experienced by the sediment particle. This highlights the effect of turbulent fluctuations and contact forces on the actual dynamics of the sediment particle (the particle collides with the bed twice at tU b /D ≈ 38 and 45, which results in high amplitude oscillations in the particle vertical acceleration in panel b). The initial lift force due to the large streamwise slip velocity and high local mean shear rate, and the subsequent reduction of the bed-normal velocity as the particle takes off agree with the experimental observations of Barros, Hiltbrand & Longmire (2018), for a spherical particle in a turbulent boundary layer. The Shields parameter for this case is Sh = 0.695 (see table 1), approximately 20 times larger than the corresponding critical value, and a criterion based on the Shields parameter can successfully predict the fully suspended trajectory of the particle. Let us now investigate the sediment transport in saltation mode. The time history of the normal position of the particle centre for the four realisations of case Ga150, all with incipient motion, is depicted in figure 7. One of these fields, 'initial condition 1', is selected for further investigation. Figure 8  FIGURE 8. Time history of the bed-normal position of the centre of the sediment particle for Ga150, coloured by the bed-normal fluid force exerted on the particle F y , scaled with (ρu 2 τ A p ); A p = πD 2 /4 is the projected area of the particle on x-z plane. The inset shows a detailed view of the initial stage of particle motion. particle for Ga150. The colour represents the fluid surface force, normalised by the boundary shear-induced force. This case shows a cyclic succession of bounces and wall collisions, with an average amplitude of 0.6D and duration between two subsequent contacts with the bed of approximately tU b /D = 5; note that the duration of resuspensions for Ga25 is approximately tU b /D = 80, i.e. 16 times larger. During most of the motion, the integral of the pressure and shear stress force is one order of magnitude larger than the estimate of the shear-induced lift force, which characterised the particle motion for Ga25. The particle lift force is of the order of the gravitational force (as discussed below), and therefore larger than the shear-induced force at this higher value of the Galileo number. The difference with the case Ga25 is attributed to the increased role of the history (or Basset) force and added mass, because of the larger accelerations. A relatively large downwards vertical force is exerted on the particle, which is due to the reaction of the fluid to the particle's sudden acceleration just after the collision with the bed, hence a large drag force; these large negative values are indicated in blue colour in figure 8. The estimated shear-induced boundary force is at most 30 times smaller than the bed-normal component of the integral hydrodynamic force, and does not seem to significantly affect the saltating motion of the particle. This force is, however, responsible for the initial particle dislodgement, see the insets of figures 5 and 8, which show that the magnitude of the dimensionless lift during the initial stages of motion for runs Ga25 and Ga150 is comparable.
The Shields parameter for this case is close to the corresponding critical value (see table 1), i.e. it is expected for the particle to be close to the onset of motion according . Time history of bed-normal position of the centre of the free particle (black) together with (a) bed-normal velocity of the particle, v p , (b) bed-normal particle acceleration a y,p , normalised with the gravitational acceleration for case Ga150. The dashed line depicts a y,p /g = −1. (c) Streamwise component of the collision force between the sediment particle and the bed.
to this classical criterion. This is consistent with the observed prompt resuspension, which does not depend on the initial turbulent field. Figure 9 illustrates the time history of the bed-normal position of the sediment, together with its bed-normal velocity (panel a), acceleration (panel b) and the streamwise bed-collision force exerted on the sediment particle (panel c). For each cycle, the dynamics qualitatively resemble those of the bouncing motion of a spherical particle colliding onto a planar wall in a quiescent liquid (Gondret, Lance & Petit 2002). First, at each impact with the bed, the sediment particle gains a high acceleration due to the collision. This sudden acceleration is quickly damped by the drag force, as shown by the fact that the negative acceleration is slightly larger in magnitude than the gravitational acceleration, g. Second, the gravitational force becomes dominant and the velocity decreases almost linearly. At this stage, the difference between the gravity and the shear-induced lift force dictates the particle acceleration, which is approximately 0.6g, see figure 9(b). As the particle gets closer to the bed and shear-induced lift force increases, the particle acceleration decreases slightly, see figure 9(b). The bouncing amplitude always decreases due to viscous dissipation, unless it happens that the sediment particle undergoes an oblique collision with high (negative) streamwise component, due to the bed geometry, see figure 9(c). In this case, the sediment slip velocity increases, and saltation is triggered by a mechanism similar to that discussed in the previous section for case Ga25.
These two cases, resuspension and saltation, quantitatively showed similar responses to different initial turbulent fields given the same geometry and mean flow; for each Galileo number we tried four different initial conditions as reported in the table 1. The turbulent fluctuations appear to cause just local distortions of the particle dynamics, leaving the main features unaffected. In the next section, we will investigate the particle dynamics at higher Galileo numbers. In this case, instead, turbulent fluctuations cause the initiation of particle motion.

Higher Galileo numbers: incipient motion and rolling
We explore the dynamics of heavier sediment particles by further increasing the Galileo number, still keeping constant the Reynolds number, particle size and the density ratio, i.e. the ratio between the particle submerged weight and the viscous forces is increased, while the relative size of the particle to the flow structures and the particle response to the flow (e.g. added mass effects) are the same. Two additional cases are considered, with Ga = 250 and 310, while fixing the remaining governing parameters. For the sediment at Ga = 250, we consider four different initial turbulent fields and let the simulation run for a time interval of approximately tU b /D = 250. Incipient motion is detected in three of these simulations, while no movement was observed in the fourth case. Case Ga310 is assessed in the same way, but using 12 different turbulent initial conditions. Particle movement was observed in 50 % of the runs, which is an indication that this case may be used as a reference for the onset of sediment motion for the parameters considered in this work. Indeed, test simulations with slightly larger values of the Galileo number showed no particle movement. The dependency of the response of the sediment particle on the turbulent initial field implies that the particle movement can only originate from extreme turbulent events, unlike the cases discussed in the previous section. Figure 10 shows the time history of the streamwise position of the particle centre for the cases when particles are dislodged for Ga = 250 in panel (a) and Ga = 310 in panel (b). Among these realisations, one for each value of the Galileo number is chosen for further investigation. Figure 11 shows the time history of the streamwise position of the sediment for Ga250 (panel a) and Ga310 (panel b), coloured by the streamwise flow-induced force on the sediment, while figure 12 shows the time history of the normal position of the sediment for Ga250 (panel a) and Ga310 (panel b), coloured by the normal flow-induced force on the sediment. The most obvious difference between these cases and those of the previous section (Ga25 and Ga150) is that now the particle is in Ga310; the averaged velocities are defined as the mean streamwise travelled distance divided by the total duration of the simulations. Consistently, the streamwise fluid surface force which moves the particle is at most 30 times larger than the boundary shear force scale (ρu 2 τ A p ) in case Ga250, and at most 50 times larger for case Ga310. On the other hand, the normal force when the sediment particle starts rolling, e.g. at tU b /D ≈ 210 for Ga250 and tU b /D ≈ 50 for Ga310, is approximately 20 times the estimated shear-induced boundary force ρu 2 τ A p (see figure 12). We performed two additional simulations, to investigate the effect of the turbulent activity on the dynamics of the sediment particle with Galileo numbers of 250 and 310. Figure 13 shows the flow configuration of these cases. In these simulations, the bottom porous bed is replaced by a geometrically rough wall, while the bulk Reynolds number based on the flow depth was kept constant. The roughness elements consist of hemispheres with the height of 13 in viscous units. Due to the effect of the vicinity of the wall to the hemisphere top elevation, the streamwise and normal fluctuating components of the fluid velocity have increased by 3 % and the friction velocity u τ by around 9 % compared to cases with the porous bed (Chan-Braun, García-Villalba & Uhlmann 2011), i.e. while the Galileo number for the cases with the porous bed and the roughness elements are the same, the Shields parameter of the cases with roughness elements has increased by 19 %. In both cases, the regimes of the sediment transport resemble those of the porous bed, with approximately 50 % larger streamwise displacement of the sediment particle. These results show the relative importance of the streamwise component of the turbulent fluctuations compared to its normal counterpart. Although both velocity components have increased by 3 %, the effect on the dynamics of the sediment particle is the increase of the streamwise displacement, while the particle still sustains the contact with the bed for the whole duration of the movement. The results of these simulations motivate us to look for those events which impose strong streamwise velocity fluctuations, and thus have the most impact on the dynamics of a sediment particle for the parameters considered here. Indeed, the recent experimental study by Cameron, Nikora & Marusic (2019) shows that, for a sediment particle with high protrusion, the drag force variance is determined by the correlation between the streamwise velocity component and the drag force, quantified via a coherence function.
To understand the events causing incipient motion, let us first inspect the turbulent structures near the sediment particle just before the dislodgement. We will focus on the six instances detected for case Ga310. Figure 14(a) shows the Q-criterion (i.e. isosurfaces of positive second invariant of the velocity gradient tensor, Hunt, Wray & Moin (1988)), coloured by the bed-normal velocity fluctuations, while the sediment particle is depicted in green. From the plot, one can clearly recognise a hairpin-like structure (see Zhou et al. 1999, for more details on the role of these structure in wall-bounded turbulence) upstream of the sediment particle as a common feature for all instances. As the outer part of the hairpin leg, inducing the so-called sweeps, passes over the particle, dislodgement is caused. This effect is quantified in panel (b) of the same figure. This panel shows the scatter plot of the instantaneous streamwise and bed-normal velocity fluctuations in a subdomain of size 5D × 3D × 3D, centred at [x, y, z] = [−3D, 2D, 0], right upstream the sediment particle. Indeed, the flow sampled inside this box corresponds to that of a strong sweep; i.e. a higher than mean streamwise velocity, and negative bed-normal velocity. This is in line with previous experimental observations that sweeps initiate the sediment motion (see e.g. Nelson et al. (1995), Hofland (2005) and Dwivedi et al. (2010)). The resulting pressure distribution on the surface of the sediment particle for the same instant is shown in panel (c) of the same figure. We clearly observe two high-pressure spots near the surface of the particle: one upstream, and a second one right below, close to the contact point of the sediment particle with the bed. These hotspots are caused by the increase in stagnation pressure due to the increased streamwise velocity. The forces created by this pressure imbalance result in a net torque which effectively displaces the particle. Figure 15 gives a clearer overview of the evolution of the pressure distribution close to the surface of the particle. Panel (a) shows the time history of the spanwise angular velocity of the sediment particle from one of the detected incipient-motion events for case Ga310. The blue vertical lines depict the three moments for which the normal stress contours are illustrated in panel (b): (i) tU b /D = 5 time units before the incipient motion, (ii) at the moment of the incipient motion and (iii) tU b /D = 5 time units after the particle started to move. The red lines in panel (b) depict the trajectory of the sediment particle over the bed after the dislodgement. As the subpanel (i) shows, the pressure distribution before the incipient motion consists of a stagnation high-pressure zone upstream and a negative pressure zone downstream of the sediment particle due to its wake. At the moment of the incipient movement (subpanel ii), the stagnation pressure increases significantly by the incoming sweep and a second hotspot appears close to the contact point of the sediment particle with the bed. Shortly after the incipient motion, as the hairpin like vortex shown in figure 14(a) is passing the sediment particle, the magnitude of the pressure in the stagnation zone decreases (subpanel iii).
To better quantify the extreme nature of the pressure distribution at the moment of the dislodgement, we present in figure 16(a) the time-averaged normal stress contours, plotted close to the surface of the sediment particle. To obtain this contour plot, the position and orientation of the sediment particle are fixed at the resting pocket depicted in figure 1(a) and 110 samples are collected over a time period of tU b /D = 2939. Comparing the time-averaged contours to those of the instances of incipient motion (figure 14c), shows that due to the strong upstream sweep, the magnitude of the high-pressure zone is significantly larger, spanning a larger surface area. Figure 16(b) shows the probability density function (p.d.f.) of the pressure distribution in a spherical shell around the sediment particle (with a slight outward Instances of incipient motion detected in the case Ga310. (a) Isosurfaces of Q-criterion, coloured by v (flow direction is from left to right). The sediment particle is depicted in green; (b) the scatter plot of the streamwise and bed-normal velocity fluctuations in a subdomain of 5D × 3D × 3D upstream of the sediment particle, with the red mark depicting the averaged value; (c) contours of the normal stress plotted on the surface of a sphere with same centre as the sediment particle and radius D/2 + 3 x.
(point of view of an observer aligned with the positive streamwise direction). offset of 3 x to ensure that only pressure values outside the particle are sampled). The samples are collected from the same data used for the results in panel (a), with the points inside the particles forming the bed excluded from the sampling. The instances of negative pressure in the wake of the particle cause the negative tail of the p.d.f., while the positive tail of the curve illustrates the high-pressure realisations around the stagnation point due to the upstream flow. Indeed, events of high stagnation pressure, with p > 0.3ρU 2 b , responsible for the incipient motion, are extremely rare.
Single sediment dynamics in turbulent flow over a porous bed 893 A24-21 FIGURE 17. Sketch of the forces acting on the sediment particle and their lever arms, used in (3.2). Here, O 1 and O 2 are the contact points between the sediment particle and the neighbouring bed particles, F x is the streamwise fluid-induced surface force, F b the particle submerged weight, X and Y the lever arms and θ the angle between the streamwise direction and the normal to O 1−2 . (a) Top view, (b) side view normal to O 1−2 .

Criteria for incipient motion
Having described in detail the cause for turbulence-induced incipient motion, we now focus on assessing different major criteria for the sediment incipient motion in this rolling mode. One of the main difficulties in predicting the onset of incipient motion is determining the critical force above which the sediment starts to roll (Ghodke & Apte 2018). As a first criterion, we use a critical value for the streamwise force based on the moment balance about the contact points of the sediment particle, similar to the method used in Celik et al. (2010). On should note that this criterion is based on the actual instantaneous hydrodynamic forces, and is not a general criterion that is simply based on global quantities, like the Shields parameter. Figure 17 shows a schematic of the sediment particle and the bed configuration together with the considered forces and their respective lever arms about the contact points. Cross-sectional components of the fluid force are neglected for now, to test the criterion based only on the streamwise component. This is in line with the work of Schmeeckle, Nelson & Shreve (2007), stating that for particles highly exposed to the flow (as in the present work), the streamwise force governs the onset of motion. By considering a moment balance about the contact points (O 1 and O 2 in figure 17), equation (3.2) yields the critical streamwise force, F x,cr , where f v = [1 + 0.5(ρ/ρ p − ρ)] is the hydrodynamic mass coefficient (Papanicolaou et al. 2002), F b is the particle submerged weight, X and Y are lever arms for horizontal and vertical forces and θ = π/6 is the angle between the streamwise direction and the normal to O 1−2 .
A more recently refined criterion suggests that the combined effect of the magnitude and duration of the hydrodynamic forces, quantified by the force impulse (i.e. the time integral of the flow-induced force), determines whether the particle will roll (see e.g. Diplas et al. (2008), Celik et al. (2010) and Celik et al. (2014)); in other words, although the force magnitude be larger than the critical value for some time, the duration may not be sufficient to induce the rolling motion.

A24-22
A. Yousefi, P. Costa and L. Brandt FIGURE 18. Time series of the streamwise fluid surface force exerted on the particle for four of the different instances of the incipient movement in the case Ga310 (solid blue line). The dashed red line depicts the critical force value based on (3.2), the cyan line shows the total moment of the fluid-induced force and the particle submerged weight around O 1−2 , and the dashed black line shows the spanwise angular velocity of the sediment particle.
To compare these two criteria, figure 18 shows the history of the streamwise fluid force, the total moment of all forces acting on the sediment particle and the spanwise angular velocity for four instances of incipient motion at Ga310. As one can infer from the angular velocity, the time history is shown until the sediment particle starts to roll. The cyan solid line depicts the moment of the total fluid force (all components) together with that of the particle submerged weight at the contact points (around O 1−2 ). The data clearly show that the incipient motion occurs as soon as the torque of the fluid-induced surface force overcomes that of the particle submerged weight. The streamwise force history is compared to the corresponding critical value F x,cr obtained from (3.2). In all cases, once the critical force value is reached (i.e. the solid blue line crosses the dashed red line), the particle rolls, which is consistent with the force-based criterion. This confirms that considering only the streamwise component of the fluid-induced force suffices for predicting the onset of sediment motion. The fact that the streamwise component of the fluid-induced force is predominant for the dislodgement of highly exposed sediment particles can also be inferred from figure 14(b), because at the start of rolling the positive streamwise velocity fluctuation contributes much more to the Reynolds shear stress than the negative bed-normal component. From the present dataset, we could not observe any instance where the force-based criterion fails to predict dislodgement, while the impulse-based criterion succeeds. In other words, we did not observe a force larger than the critical value that could not dislodge the particle due to a too short duration. Yet, we should note that the sampling time in the present work is much smaller than that of the experiments of Celik et al. (2010), due to the high costs of interface resolved simulations. Moreover, the ratio between the particle diameter and the smallest flow structures is approximately 20 times larger in the high-Reynolds-number experimental study (D + ≈ 20 in the present work, compared to ≈400 in the experimental study).
To confirm the hypothesis that a force-based criterion for sediment dislodgement may fail when a larger particle is considered, we have performed an auxiliary DNS with a particle with diameter D twice as large, while keeping the bulk Reynolds number and Shields parameter equal to those of case Ga310, i.e. a sediment particle with size of D + ≈ 40 and Galileo number of 620. Figure 19(a) shows a segment of the time history of the streamwise flow-induced force on the particle after its release, while the surface contours of the pressure for an instance with a streamwise fluid force greater than that of the critical value, are illustrated in the panel (b) of the same figure. We observe here and in several other instances that the force exceeds the critical value, but the particle remains stationary in its resting pocket. The shaded area in the figure represents the impulse, which the criterion in the study of Celik et al. (2010) is based on. Indeed, when considering a larger particle we have found instances where, though the streamwise force is sufficiently high to promote dislodgement, the time scale at which it acts is not sufficiently long. Chan-Braun et al. (2011) and Mazzuoli & Uhlmann (2017) among others showed that the interaction of small vortices with a relatively large sphere (i.e. for large values of D + ) results in fluctuations of the force acting on the particle due to the filtering of vortices smaller than the particle size. Therefore, the dynamics of the particles is mostly driven by vortices of length scale comparable with that of the particles, so larger particles and FIGURE 20. Spectra of the streamwise fluid surface force exerted on the particle for different particle size: D + = 21.5 (blue) and D + = 43 (red), normalised by the critical force value for each size. The vertical lines denote the angular frequency at each particle size.
larger vortex sizes are associated with stronger forces and longer time scales. This is quantified in figure 20, where the spectra of the streamwise fluid surface force exerted on the particle normalised by the critical force value for the two different particle sizes considered (D + = 21.5 in blue and D + = 43 in red) are compared. The vertical lines in the figure denote the angular frequency at each particle size. The data in the figure show that, for the bigger particle, the effect of the flow structures with a time scale comparable or larger than that of the bigger particle (angular frequencies smaller than the red vertical line) is larger than the effect of structures of similar time scale on the smaller particle, which confirms that the filter at the particle scale dissipates more energy in the case of smaller particles.

Conclusions
We have used interface-resolved simulations to study the dynamics of a sediment particle, initially at rest on a fixed porous bed, in a turbulent open channel flow. A relatively small particle (diameter D ≈ 20ν/u τ ) is subjected to a turbulent flow, for different relative strengths of buoyancy to shear-induced boundary force, quantified by the Shields parameter. In practice, we have explicitly varied the Galileo number, since it can be computed a priori from the flow governing parameters.
We have shown that, for the cases with lower Galileo numbers (below 150), the sediment dislodgement is driven by the mean flow: the particle promptly lifts-off when exposed to the turbulent field, irrespective of the specific turbulent initial condition. Two classical modes of sediment transport, resuspension and saltation, have been reproduced by the simulations and analysed.
In the resuspension mode, the particle travels suspended in the flow with an oscillatory path. These oscillations are explained by the balance between the shear-induced lift force, which decreases when increasing the distance from the wall, and gravity. This regime is characterised by a flow-induced lift force of the same order of magnitude as the particle submerged weight ((ρ p − ρ)gV p ) and the bed shear-induced force (ρu 2 τ A p ), that the particle weight is mainly balanced by the lift force induced by the mean flow. Moreover, we show that although the bulk dynamics Single sediment dynamics in turbulent flow over a porous bed 893 A24-25 of the sediment particle is correlated with the apparent slip velocity and the local shear rate experienced by the particle, the effect of turbulent fluctuations and contact forces on the particle acceleration is not negligible.
The saltation mode, case Ga150, is characterised by a high-frequency bouncing trajectory, with its vertical component resembling the bouncing motion of a sphere colliding onto a planar wall (Gondret et al. 2002). Here the acceleration is almost constant throughout most of the motion (approximately 0.5g), except very close to the bed. Then, as the particle approaches the wall, it experiences short-range hydrodynamic interactions, and ultimately collisions which induce, expectedly, a very strong vertical acceleration. This positive acceleration is quenched by the gravitational force and by an increased drag until the particle starts falling again. In this case, the fluid surface force has the same order of magnitude as the particle submerged weight, but is almost 30 times larger than the shear-induced boundary force. Despite this, whenever the sediment particle has a significant apparent slip velocity (with respect to the mean flow), the shear-induced lift force plays an important role on the sediment particle dynamics. This is the case when the streamwise velocity of the particle is reduced considerably due to a collision (see figure 9c), or in the beginning of the simulation, since the particle is at rest.
We have then considered cases where the flow governing parameters are such that the particle motion is dictated by extreme turbulent events. Two different values of the Galileo number, 250 and 310, have been examined here. The main difference between the trajectories for the two cases is the resting times, and the probability of onset of motion. The Shields parameter is approximately 0.005 for these cases -much smaller than the corresponding critical value. The flow events triggering the sediment motion are captured by sampling the flow field upstream of the particle, just prior to its movement. Consistently with previous experimental studies (see e.g. Nelson et al. (1995), Hofland (2005), Dwivedi et al. (2010)), our results show that intense streamwise velocity fluctuations, associated with the sweeps, are the main cause of sediment motion. Visualisations of the flow structure at the moment of incipient motion showed that these events are associated with the tail of a hairpin-like vortex, travelling near the particle. Taking further advantage of the three-dimensional nature of the DNS data, we have investigated the pressure distribution close to the surface of the sediment particle at the moment of incipient motion. For all the instances detected, this distribution shows a consistent footprint of the passage of the hairpin vortex: an increased stagnation pressure zone upstream due to a higher than mean streamwise velocity fluctuation, and a high-pressure region close to the contact point of the particle with the bed.
Finally, previously proposed criteria for the onset of sediment transport, based on the concept of a critical impulse (Celik et al. 2010), have been tested against the numerical results for these cases at higher Galileo numbers. Our results clearly demonstrate that, for the present set-up, a critical threshold streamwise force effectively predicts the incipient motion. Although the principle behind the impulse-based criteria for predicting the onset of sediment transport seems more rigorous, in the present set-up a force exceeding the corresponding critical value was always able to move the sediment particle, regardless of its duration. This might be a consequence of the relatively low Reynolds number of the DNS, where extreme turbulent events with a typical time scale of D/u τ always act on a long enough time to move the sediment particle. Indeed, a test simulation with a larger particle size at the same Reynolds number showed instances of the force exceeding the critical value, while the sediment particle remained stationary. Thus, simulations at higher Reynolds number and varying particle sizes could provide further insights on the necessary conditions for turbulence-induced impulses capable of dislodging relatively larger sediment particles.