Turbulent and viscous sediment transport – a numerical study

Sediment transport is studied as a function of the grain to fluid density ratio using two phase numerical simulations based on a discrete element method (DEM) for particles coupled to a continuum Reynolds averaged description of hydrodynamics. At a density ratio close to unity (typically under water), sediment transport occurs in a thin layer at the surface of the static bed, and is called bed load. Steady, or ‘saturated’ transport is reached when the fluid borne shear stress at the interface between the mobile grains and the static grains is reduced to its threshold value. The number of grains transported per unit surface therefore scales as the excess shear stress. However, the fluid velocity in the transport layer remains almost undisturbed so that the mean grain velocity scales with the shear velocity u∗. At large density ratio (typically in air), the vertical velocities are large enough to make the transport layer wide and dilute. Sediment transport is then called saltation. In this case, particles are able to eject others when they collide with the granular bed. The number of grains transported per unit surface is selected by the balance between erosion and deposition and saturation is reached when one grain is statistically replaced by exactly one grain after a collision, which has the consequence that the mean grain velocity remains independent of u∗. The influence of the density ratio is systematically studied to reveal the transition between these two transport regimes. Finally, for the subaqueous case, the grain Reynolds number is lowered to investigate the change from turbulent and viscous transport.


Introduction
Despite a wide literature, some fundamental aspects of sediment transport in turbulent flows are still only partly understood.In particular, derivations of transport laws, relating the sediment flux to the flow velocity, have a strong empirical or semi-empirical basis (see e.g.among many others Meyer-Peter and Müller (1948), Ribberink (1998), Camemen and Larson (2005), Greeley et al. (1996), Iversen and Rasmussen (1999), Kok and Renno (2009) and references therein), thus lacking more physics-related inputs.Also, the dynamical mechanisms limiting sediment transport, in particular the role of the bed disorder (Charru, 2006) and turbulent fluctuations (Marchioli et al., 2006;Baas, 2008;Le Louvetel-Poilly et al., 2009), remain matter of discussion.
Here we investigate the properties of steady homogeneous sediment transport using a novel numerical description of particle-laden flows, using two-phase numerical simulations based on a discrete element method for particles coupled to a continuum Reynolds averaged description of hydrodynamics.In particular, we examine the transition from bed-load to saltation by studying the influence of the grain to fluid density ratio ρ p /ρ f .A similar approach has recently been used to study the onset of aeolian saltation (Carneiro et al., 2011).This paper mostly summarises the talk addressing these questions, given at the workshop SALADYN, Institut de Physique du Globe de Paris, 5-7 November 2012.More details on this work, as well as a more developed bibliography on the subject, can be found in Durán et al. (2012).

The model
The idea is to use a continuum description of hydrodynamics, averaged at a scale larger than the grain size.This means that the feedback of the particles on the flow is treated in the mean field manner.This method allows us to perform very long numerical simulations (typically 1000 √ d/g), using a (quasi) 2D large spatial domain (typically 15000 spherical grains in a xyz box of respective dimensions 1000 d × 1 d × 1000 d), while keeping the complexity of the granular phase.Periodic boundary conditions are used in the x (flow) direction.We will now detail the different ingredients of the model -see table 1 for notations.

Forces on particles
The grains have a spherical shape and are described by their position vector r, velocity u and angular velocity ω.A given grain labelled p inside a fluid obeys the equations of motion, where g is the gravity acceleration, I = md 2 /10 is the moment of inertia of a sphere, f p,q is the contact force with grain q, n p,q is the contact direction, and f p fluid encodes forces of hydrodynamical origin.
We model the contact forces following a standard approach for the modeling of contact forces in MD codes (see e.g.DEM book (2011) and references therein), where normal and tangential components are described by spring dashpot elements.A microscopic friction coefficient is also introduced but cohesion is neglected.For simplicity we assume that the net hydrodynamical force (f p fluid ) acting on a grain p due to the presence of the fluid is dominated by the drag and Archimedes forces, f p drag and f p Arch , respectively.The lift force, lubrication forces and the corrections to the drag force (Basset, added-mass, Magnus, etc.) are neglected, as they do not have any qualitative influence on the results presented here.
Drag force -We hypothesize here that the drag force exerted by a homogeneous fluid on a moving grain only depends on the difference between the grain velocity u p (x, z) and the fluid velocity u(z) at grain's height z.Introducing the Reynolds number R u based on this fluid-particle velocity difference R u = |u − u p |d/ν, the drag force can be written under the form where C d (R u ) is the drag coefficient.We use the following convenient phenomenological approximation Ferguson and Church, 2004).C ∞ d 0.5, is the constant drag coefficient of the grain in the turbulent limit (R u → ∞).R c u 24 is the transitional Reynolds number below which the drag force scales linearly with the velocity difference.
Archimedes force -This force results from the stress which would have been exerted on the grain, if the grain had been a fluid.Thus, where π 6 d 3 is the grain volume and σ f ij = −p f δ ij +τ f ij is the undisturbed fluid stress tensor (written in terms of the pressure p f and the shear stress tensor τ f ij ).In first approximation, the stress is evaluated at the center of the grain.

Hydrodynamics and coupling
In the presence of particles occupying a volume fraction φ, the hydrodynamics is described by the two-phase flow Reynolds averaged Navier-Stokes equations: where D t u i ≡ ∂ t u i +u j ∂ j u i denote the fluid inertia.τ f ij is the total shear stress tensor resulting both from viscous diffusion of momentum (viscous stress) and transport of momentum by turbulent fluctuations (Reynolds stress).F is the body force exerted by the grains on the fluid.In the steady and homogeneous case investigated here, These RANS equations simplify into where we note τ f = τ f xz the fluid shear stress, and later on u = u x for the fluid horizontal velocity.
and turbulent shear stress The coupling term F can then be obtained by averaging the hydrodynamical force f p fluid acting on all the grains moving around altitude z, in a horizontal layer of area A and thickness dz: We take for A the total horizontal extent of the domain (i.e.1000d × 1d).The symbols .denote ensemble averaging.
Here, we retain its x-component only, which simplifies into where the grain's volume fraction φ is defined as Eq. 6 integrates as τ f (z) = ρ f u 2 * − τ p (z), where we have introduced the shear velocity u * , defined by the undisturbed (grain free) wall shear stress, and the grain borne shear stress τ p , computed from the integration of ( 8) over sufficient vertical extension to count all moving grains.
In order to relate the fluid borne shear stress to the average fluid velocity field, we adopt a Prandtl-like turbulent closure.Introducing the turbulent mixing length , we write ν is the viscosity (a constant independent of the volume fraction).As for the mixing length , we know it should vanish below some critical Reynolds number R c and should be equal to the distance to the surface z, far above the transport layer.
To avoid the need of a somewhat arbitrary definition of an interface between he static and mobile zones of the bed, we propose the differential equation where κ 0.4 is von Karman's constant.We have checked that, in the case of a turbulent flow over a smooth and flat surface (no grains), we recover the prediction computed with the phenomenological expression for the mixing length suggested by van Driest (Pope, 2000), which reproduces well classical experimental results.Comparison to measurements determines the dimensionless parameter R c 7. Other empirical expressions for the exponential term in Eq. 11, e.g.
, with other values of γ give qualitatively similar results.
Starting integration deep enough in the static bed to be in the asymptotic limit that can be analytically derived, we obtain the different hydrodynamical fields.They are displayed in Fig. 1, in the case of subaqueous transport (ρ p /ρ f = 2).

Sediment flux
Steady and homogeneous sediment transport is quantified by the volumetric saturated flux q sat , i.e. the volume of the particles (at the bed density) crossing a vertical surface of unit transverse size per unit time.It has the dimension of a squared length per unit time.In the simulations, we compute it as A key issue is the dependence of q sat on the shear velocity or, equivalently, its dimensionless counterpart the Shields number , defined by simulations presented here have been performed for R = 10, except in section 6 where viscous regime is investigated with R = 0.1.We show in Fig. 2 the saturated flux in both cases (water and air).In agreement with experimental observations (e.g.Meyer-Peter and Müller (1948), Ribberink (1998), Lajeunesse et al. (2010), Rasmussen et al. (1996), Creyssels et al. (2009)), we find that q sat scales asymptotically as (or u 2 * ) for saltation, while q sat scales as 3/2 (or u 3 * ) underwater.This figure also reveals the existence of a threshold shear velocity below which the flux vanishes.More precisely, we define the dynamical threshold Shield number d from the extrapolation of the saturated flux curve to 0, which gives in our case d 0.12 for water (ρ p /ρ f = 2) and d 0.004 for air (ρ p /ρ f = 2000), respectively.These values are consistent with experimental ones within a factor of 2 (see the data gathered in the review of Durán et al. (2011)).

Mechanisms at work within the transport layer
Bed load and saltation mainly differ by the vertical characteristics of the transport layer.At small density ratios the motion of grains is confined within a thin layer of few grain diameters.By contrast, for large density ratios, grains experience much higher trajectories: the transport layer is much wider and the flux density decreases exponentially with height with a characteristic size of the order of 50d, roughly independent of the shear velocity.The transport layer thickness is effectively determined by the hop length for ρ p /ρ f 10.Below this cross-over value, this thickness is given by the grain diameter d, as trajectories are almost horizontal.The transition from bed load to saltation therefore takes place when the vertical velocities of the particles are sufficiently large for these particles to escape the traps formed by the grains on the static bed.
Another difference between bed load and saltation is how the grain's feedback on the flow is distributed within the steady state transport layer.Fig. 3 presents the vertical profiles of the fluid shear stress, rescaled by the dynamical threshold τ d = d (ρ p /ρ f −1)gd (as defined by the saturated flux), for different shear velocities.For bed load (Fig. 3a), the different profiles of the fluid shear stress seems to converge to the threshold value very close to the surface (z = 0).In this transport layer, the fluid momentum decays over few grain sizes, in agreement with the vertical extension of the transport layer.By contrast, the fluid shear stress is below the threshold in the bed (z < 0) but some (weak) transport still occurs there, which is sustained not by the fluid itself but by the momentum transferred to the surface by grain col-This general picture is still valid for saltation (Fig. 3b), however now the dynamical threshold is reached much farther from the surface (at z 10 d) which implies that the kinetic energy of impacting grains is large enough as to sustain the transport below this height.Above it, the transport is driven by the fluid and most of its momentum is dissipated in a much larger layer (comprising tens of grain diameters) again in agreement with the size of the saltation layer.Notice that although this surface sublayer below 10 d contains most of the grains, it still represents a small fraction of the overall transport layer.
An important consequence of this distinction in the vertical structure of the grain's feedback is that although for bed load transport is equilibrated when the fluid shear stress reaches its dynamical threshold below the transport layer, this condition is not enough for saltation to equilibrate.For saltation there is a sub-layer where transport is not directly driven by the fluid and thus its equilibration is not dictated by the threshold.There, the properties of grain's collisions become relevant and the equilibrium is described by the conservation of the number of saltating grains i.e. when the number of grains entering the flow exactly balance those grains trapped by the bed.

Scaling laws
The saturated flux can then be decomposed as the product of the number n of transported grains per unit area by the mean grain horizontal velocity ūp : q sat = 1 These quantities are plotted as functions of the Schields number in Fig. 4. A scaling law n ∝ − d is well verified over two decades, independently of ρ p /ρ f .By contrast, the density ratio has a strong effect of ūp .The mean grain velocity is independent of for large ρ p /ρ f (aeolian case), whereas it varies linearly with the fluid shear velocity at low density ratio (sub-aqueous case).Interestingly, ūp remains finite at the threshold, at a value independent of ρ p /ρ f .These behaviours are in agreement with experimental findings in the case of bedload (Lajeunesse et al., 2010).We can derive these scaling laws from simple models.Following Bagnold's original ideas for the case of bedload (Bagnold, 1956), we write the grain born shear stress τ p as proportional to the moving grain density n and to the drag force f d acting on a moving grain.As these grains are in steady motion, f d balances a resistive force due granular friction, collisions with the bed, etc.These different dissipative mechanisms can be modeled as an overall effective friction force characterized by a friction coefficient µ d , leading to f d = π 6 µ d (ρ p − ρ f )gd 3 .Saturation is reached when the fluid shear stress equals the transport threshold at the surface of the static bed, i.e. when τ p = ρ f u 2 * − τ d , with As consequence, the number of transported particles per unit area is solely determined by the excess shear stress: n = (ρ f u 2 * − τ d )/f d .Assuming that the transported grains do not disturb the flow, the flow velocity around grains u must be proportional to the shear velocity, so that u/u d = √ / d .One can then deduce: , where µ s is a friction coefficient characterising the drag force necessary to set into motion a static grain.This predicts that the grain velocity does not vanish at the threshold, if friction is lowered during motion (µ d < µ s ).The velocity at threshold can be interpreted as the velocity needed by a grain to be extracted from the bed and entrained by the flow.
We can proceed in a similar manner for the aeolian saltation regime, following ideas initially proposed by Owen (1964) and Ungar and Haff (1987).The momentum balance τ p = ρ f u 2 * − τ d still holds, so that n has the same form as in the bed-load case, but with a different effective drag force f d , not related to friction anymore but to grain velocities.For saltation, steady transport also implies that the number of grains expelled from the bed into the flow exactly balances those trapped by the bed, i.e. a replacement capacity equal to one.Due to the grain feedback on the flow, in contrast with bed load, grains in the transport layer feel a flow independent of the wind strength (see Fig. 3b).Thus, new moving grains come only from high energy bed collisions.Since the number of ejected grains is a function of the impact energy (or equivalently, of the impact velocity), the mean grain velocity ūp must be constant, independent of the shear velocity, scaling with u d .In fact, all particle surface velocities also scale with u d , so that f d is a constant too, leading again to n ∝ − d .
These scaling laws for n and ūp as functions of the Shields number explain the different behaviour of q sat ( ) in the  sub-aqueous bedload and aeolian saltation cases, as shown in Fig. 2.

Viscous transport
We have also investigated the case of sediment transport in the viscous regime, setting the grain Reynolds number at a small value R = 0.1 (in comparison to R = 10 for all other simulations), and for a water-like density ratio ρ p /ρ f = 2.The data displayed in Fig. 5a show that the sediment flux increases more rapidly with the Shields number than in the turbulent case for which q sat ∼ 3/2 above the threshold.This is mainly due to the behaviour of the mean grain velocity whose scaling with the shear velocity

√
is not linear any more but rather quadratic (Fig. 5c).By contrast, the linear scaling of the moving grain density n with − d is still fairly verified, except close to the threshold (Fig. 5b).
These dependences of ūp and n on are what we can expect in the viscous case.In this regime, the flow velocity at the distance z = d close to the bed is u u 2 * d/ν, so that u/u d = / d .From the expression of the drag force in the limit R u 1, one can then deduce: in agreement with the data.As for the moving grain density, the momentum balance τ p = nf d = ρ f u 2 * − τ d , leading to n ∝ − d should hold.The reason for the discrepancy between the data and this scaling law close to the threshold could be due to a non-frictional behavior of the transport layer, but clearly necessitates further investigation.Finally, taking the product of n and ūp to get q sat , we obtain a flux that scales asymptotically like 2 , which is in agreement with the experiments of Charru et al. (2004) performed with a bed of particles sheared by a viscous flow.
Extending these experiments, these authors have further investigated the flow and velocity profiles inside and close to the bed of particles (Mouilleron et al., 2009).We can directly compare their data to the simulations.As shown in Fig. 6, the overall agreement is good.One can, however, note that the model predict slightly larger velocities very close to the bed than in the experiment, both for the fluid and the particles.Again, this requires further studies in this regime.

Conclusions
The aim of this paper was to present a novel numerical approach for sediment transport based on a discrete element method for particles coupled to a continuum Reynolds averaged description of hydrodynamics.We have studied the effect of the grain to fluid density ratio ρ p /ρ f and showed that we can reproduce both (sub-aqueous) bed load at ρ p /ρ f close to unity, where transport occurs in a thin layer at the surface of the static bed, and (aeolian) saltation at large ρ p /ρ f , where the transport layer is wider and more dilute.Scaling laws for the density of moving grains, and for the average velocity of these grains, as functions of the Schields number are found in agreement with experiments, and support simple mechanisms at work in steady and homogeneous transport.
Further work will be focused on transient situations, in order to study the time and length scales encoding the relaxation properties of out-of-equilibrium transport.In particular, it would be crucial to analyse how both n and ūp individually relax, in order to understand what are the physical mechanisms responsible for the response of the flux to a given flow.Also, it would be interesting to investigate the case of bimodal or more polydisperse grains (Houssais and Lajeunesse, 2012).

Fig. 3 .
Fig. 3. Vertical profiles of the fluid borne shear stress τ f for different values of the shear velocity ratio √ / d (see legend), in water (a) and air (b).
n ūp In the numerical simulations, we compute n and ūp as n =

Fig. 4 .Fig. 5 .
Fig. 4. (a) Number of transported grains per area n and (b) mean velocity of these grains ūp as functions of the Shields number , for different values of the density ratio ρ p /ρ f (see legend).

Table 1 .
Units used in the model, expressed in terms of the grain density (ρ p ), the fluid density (ρ f ), the gravity (g) and the mean grain diameter(d)