Rheology of mobile sediment beds sheared by viscous, pressure-driven flows

We present a detailed comparison of the rheological behaviour of sheared sediment beds in a pressure-driven, straight channel configuration based on data that was generated by means of fully coupled, grain-resolved direct numerical simulations and experimental measurements reviously published by Aussillous {\it et al.} (J. Fluid Mech., vol. 736, 2013, pp. 594-615). The highly-resolved simulation data allows to compute the stress balance of the suspension in the streamwise and vertical directions and the stress exchange between the fluid and particle phase, which is information needed to infer the rheology, but has so far been unreachable in experiments. Applying this knowledge to the experimental and numerical data, we obtain the statistically-stationary, depth-resolved profiles of the relevant rheological quantities. The scaling behavior of rheological quantities such as the shear and normal viscosities and the effective friction coefficient are examined and compared to data coming from rheometry experiments and from widely-used rheological correlations. We show that rheological properties that have previously been inferred for annular Couette-type shear flows with neutrally buoyant particles still hold for our setup of sediment transport in a Poiseuille flow and in the dense regime we found good agreement with empirical relationships derived therefrom. Subdividing the total stress into parts from particle contact and hydrodynamics suggests a critical particle volume fraction of 0.3 to separate the dense from the dilute regime. In the dilute regime, i.e., the sediment transport layer, long-range hydrodynamic interactions are screened by the porous media and the effective viscosity obeys the Einstein relation.


Introduction
Understanding and predicting the behaviour of mobile, granular sediment beds exposed to shearing flows is essential for a number of natural phenomena (e. g. sediment transport in rivers and oceans) but also for numerous engineering processes (e. g. slurry transport in the mining and petroleum industries). While many of these flows are turbulent, this is not always the case, and laminar flows are also present in many situations, such as the † Email address for correspondence: b.vowinckel@tu-braunschweig.de arXiv:2103.08935v1 [physics.flu-dyn] 16 Mar 2021 debris flow of highly concentrated suspensions or the creeping motion of soils down a hill slope (Jerolmack & Daniels 2019). The underlying physical mechanisms leading to the morphology of sediment beds, i.e. ripples and dunes, observed in the turbulent case seem also to bear many similarities to those seen in the laminar case (Lajeunesse et al. 2010). Studies for laminar cases can thus be important to create analogues of phenomena that occur in turbulent flows at larger field scales but are also of interest by themselves. The present study focuses on this laminar regime.
Modelling sediment transport requires the understanding of the rheological behaviour of the granular material (Vowinckel 2021, and references therein). Gravity plays an important role as it controls the level of stress experienced by the grains. The motion of the grains is caused by the shearing forces exerted by the fluid at the surface of the sedimented bed but the grain packing is controlled by gravity and is free to dilate as the shearing forces are increased. This rheological situation termed pressure-imposed has been the subject of significant recent advances. It has been shown that a description in terms of a frictional rheology can be applied to both dry granular flows and viscous suspensions, despite the fact that the interactions between particles may be different. Inter-particle collisions and friction between contacting particles dominate in granular flow while hydrodynamic interactions are important in viscous suspension although the role of contacts becomes increasingly predominant with increasing concentration (see e. g. Guazzelli & Pouliquen 2018).
In the inertial case of a dry granular material sheared at a shear rateγ under an imposed granular pressure p p , the rheology is determined by the particle volume fraction, φ, and the macroscopic friction, µ = τ /p p , where τ is the shear stress, which both are functions of a single dimensionless inertial number I = d pγ ρ p /p p , where d p is the particle diameter and ρ p the particle density (GDR Midi 2004;Forterre & Pouliquen 2008). A similar formalism can be applied to viscous suspensions of non-Brownian spheres but with a viscous number J = η fγ /p p in place of the inertial number I (Boyer et al. 2011), where η f is the dynamic viscosity. This frictional formulation is equivalent to the more classical presentation using viscosities depending solely on φ, where the relative shear and normal viscosities can be obtained as η s = τ /η fγ = µ/J and η n = p p /η fγ = 1/J. The transition from the viscous to the inertial regime is far from being completely understood but is supposed to occur when the Stokes number St = I 2 /J = ρ pγ d 2 p /η f is ∼ O(1 − 10) (Bagnold 1954;Ness & Sun 2015).
Traditionally, the rheology of dense suspensions has been assessed in rheometry experiments of neutrally buoyant spheres, by imposing a constant volume fraction, φ, to obtain the shear viscosity η s as a function of φ (see e. g. Stickel & Powell 2005;Guazzelli & Pouliquen 2018). The pressure-imposed rheometry described in the preceding paragraphs is a recent addition and has been found to be particularly useful in the range of large φ which is less amenable to conventional rheometry (Boyer et al. 2011;Dagois-Bohy et al. 2015;Tapia et al. 2019;Guazzelli & Pouliquen 2018). This latter approach where the suspension is free to expand (or to contract) under shear is particularly well suited to study the rheology close to the jamming transition and to measure the maximum volume fraction, φ c , where the viscosities diverge. Another valuable aspect of this pressurecontrolled rheometry is the direct measurement of the particle pressure, p p , which is usually not accessible to conventional rheometry.
The rheological measurements can be described by empirical correlations relating the shear viscosity, η s , to the volume fraction, φ (see e. g. Stickel & Powell 2005). There are also phenomenological relations for the normal viscosity, η n , versus φ, in particular that proposed by Morris & Boulay (1999) to match experimental results on shear-induced migration. Of particular interest to sediment transport wherein two different phases need to be addressed are the relations proposed by Boyer et al. (2011) as the expressions for µ(J) and φ(J) and equivalently for η s (φ) and η n (φ) contain two terms, one coming from the hydrodynamic interactions and one coming from direct contacts. The viscous term is constructed to yield the Einstein viscosity at low φ while the contact term is similar to that found for dry granular media and produces the observed power-law divergence at the jamming transition. It is important to note that these phenomenological relations are inherently empirical, because they involve adjustable parameters that have been determined by best fit to experimental data of sheared dense suspensions with particle volume fractions φ/φ c > 0.5. Hence, using these correlations to describe sediment transport may be problematic, because the bed load transport layer can easily reach values lower than the data range of the rheometry experiments. Consequently, there is a need to investigate sediment transport by means of highly resolved data to test the validity of the empirical correlations as constitutive laws.
Modelling sediment flows on a continuum scale indeed requires applying a two-phase approach, and thus using the appropriate constitutive relationships for the stresses of the fluid and particle phases is essential. This search started with the pioneering studies of Bagnold (1956), who applied the results of his rheological experiments (Bagnold 1954) to the non-uniform case of grains flowing over a gravity bed, i.e. a sediment bed stabilized by gravity. Recognising the necessity of a frictional view of the problem has been found to be instrumental. In particular, Ouriemi et al. (2009) used a frictional rheology similar to that proposed for dry granular media to describe the stress of the particle phase as they considered that the grains were mainly interacting through contact forces inside the bed. They also took a Newtonian rheology for the fluid phase with an Einstein dilute viscosity, as for grains in contact higher-order hydrodynamic interactions are shielded and the viscosity just reduced to its dilute value. Two-phase modelling using a µ(J) frictional rheology has been tested against bedload experiments in channel flows by Aussillous et al. (2013) and was found to be successful in predicting the flow inside the mobile bed. However, the rheological coefficients were adjusted to match the experimental velocity and concentration profiles and overall differ from those found by Boyer et al. (2011). Another method followed by Houssais et al. (2016) was to consider heavy particles sheared in an annular channel and to infer the rheological properties of the settled suspension. The rheology of Boyer et al. (2011) was recovered, but with the addition of a pressure at the interface of the free fluid flow and the sediment bed, which corresponded to some fraction of the weight of an individual particle. These experiments which aim at studying the rheological behaviour of mobile sediment beds are particularly difficult as they require great accuracy in the measurements of the packing fraction and of both the particle and fluid motion. There are additional difficulties coming from the choice of the shearing flows. In annular Couette flows (Mouilleron et al. 2009;Houssais et al. 2016), the thickness of the flowing suspension is merely two particle diameters, and using a continuum description may not be fully justified. The mobile layer is much thicker in Poiseuille-type flows and is better suited to a comparison with a continuum modelling. This setup is also a more realistic representation of a natural channel because unlike in annular flows the sidewalls have no curvature (Lobkovsky et al. 2008;Aussillous et al. 2013;Allen & Kudrolli 2017). However, using a straight channel configuration may also prove to be problematic as there is a continuous erosion in the channel, and consequently the generated data always contain a transient component where erosion and deposition may not be in full equilibrium.
Consequently, uncertainties remain as to how rheological models derived for Couette flows are applicable to the situation of natural sediment transport, where the flow is typically driven by a volume force, such as a pressure gradient or a downhill force. For these types of flows, the key differences are (i) total shear increasing with flow depth, which yields a wide range of viscous numbers J, (ii) settled particles where the particle pressure p p and viscous number J vary vertically, and (iii) non-homogeneous particlevolume fractions throughout the bed-load transport layer at the interface between the free-fluid flow and the sediment bed. In particular, this latter region close to the interface remains poorly accessible by means of experimental measuring techniques. A promising alternative route for investigating the rheological behaviour of sediment beds is to use direct numerical simulations (DNS) at the particle scale. There are only a few contributions in the context of beload transport. In particular, the study of Kidanemariam & Uhlmann (2014), which uses the immersed boundary method (IBM) for the fluidsolid coupling and a soft-sphere approach for solid-solid contact (DEM), investigated the flow-induced motion of a thick bed of spherical particles and found excellent agreement with the mean flow properties reported by Aussillous et al. (2013) even though the Reynolds number in the simulations was two orders of magnitude larger than in the experiment. Furthermore, Kidanemariam (2017) presented a first attempt to directly infer the rheology of a sheared sediment bed. The present paper is following this route to verify if previous rheological considerations from pressure-imposed rheometry remain applicable for the setup of sediment transport in a pressure-driven flow.
The objective of the present contribution is, hence, to quantify the stress exchange of the fluid and particle phase in the mixture to access the highly-complex fluid particle interactions in the bedload transport layer. This analysis is crucial to compute directly the rheological behavior of sediment beds exposed to a pressure-driven flow and allows for a comparison to previous rheological studies of annular Couette flows with neutrally buoyant particles. Towards this goal, we employ particle-resolved DNS using the IBM (Uhlmann 2005;Kempe & Fröhlich 2012a) and the DEM validated by Biegert et al. (2017a). The simulation framework is applied to reproduce the experimental configuration of Aussillous et al. (2013), albeit for higher flow rates but by still remaining in the viscous regime of flow with St < 10, see § 2 and 3. We then apply in § 4 the strategy of Biegert et al. (2018) and Vowinckel et al. (2019) to capture quantities unreachable in experiments such as the stress balances for the fluid and particle phases and the whole fluid-particle mixture in the streamwise and vertical directions. Rheological quantities are also inferred from these highly-resolved data and are compared to the reconstructed rheological data from the experimental data of Aussillous et al. (2013)

Experimental data
We use the experimental data of Aussillous et al. (2013) who examined the mobile layer of a granular bed in laminar flows in a rectangular-channel flow. This database yields depth-resolved profiles of particle velocity for a range of fluid heights h f as the height of the clear-liquid layer above the granular bed and flow rates in the laminar regime.
We summarise below the main features of the experimental apparatus used to obtain this database. Further details can be found in Aussillous et al. (2013). Two batches of particles, consisting respectively of borosilicate spheres having a mean diameter d p = 1.1±0.1mm and a specific density (ρ p −ρ f )/ρ f = 1.1 and of PMMA spheres having a mean diameter d p = 2.04 ± 0.03mm and a specific density (ρ p − ρ f )/ρ f = 0.1, were selected, where ρ f is the fluid density. The particles were immersed in a viscous fluid (mainly composed of Triton X-100 and water) having the same refractive index as the respective particles. A dye (Rhodamine 6G) that fluoresces when illuminated by the laser in the wavelength range greater than 555 nm was added to the fluid. The flow setup consisted of a horizontal rectangular channel of length L x = 100 cm, height L y = 6.5 cm, and width L z = 3.5 cm, where x, y, and z denote the streamwise, vertical and spanwise coordinate, respectively. To create a sediment bed, the channel was filled up with monodisperse particles and fluid, then turned upside down and tilted to consolidate the sediment bed. Afterwards, the channel is set horizontally and flipped back to its original position and a pressure gradient is applied to generate a small flow rate. This procedure was applied to fill solely the channel entrance with sedimented particles, leaving an empty buffer space near the outlet of the channel. A constant fluid flow rate, Q f , was then applied to erode the sedimented bed of particles into the empty buffer space. In this way, the fluid shear stress at the top of the bed decreased as the fluid-particle interface (upstream from the buffer region) receded (i. e. h f increased) with time. The fluid-sediment interface is detected by computing the maximum change of slope of the averaged grey level profile (green line in figures 1a and c). Several runs were conducted for each particle type, by varying the imposed flow rate. Data for the velocity profile (particle velocity, u p , and fluid velocity, u f , in solely the pure fluid zone for the PMMA particles) were collected by averaging 10 images over 0.5 s every 5 s. In addition, bulk quantities were deduced such as the time evolution of the particle bed height, h p , and the flow rate per unit width, q f,exp , by neglecting the role of the moving granular layer. The experiments were conducted in the viscous regime with a Reynolds number defined as Re = ρ f Q f /η f L z ranging between 0.2 < Re < 1.2. The Stokes number based on the shear rateγ was St = ρ p d 2 pγ /η f ≈ 0.01 on average. The most difficult quantity to extract from the experiments is the particle volume fraction, φ. In Aussillous et al. (2013), the volume-fraction profile for the borosilicate particles was evaluated by using the averaged grey-level profile, P I (y), of the same images as those utilised to infer the velocity profile and by scaling it by the grey-level profile of the immobile initial bed, which is assumed to be at the constant maximum particle-volume fraction φ c = 0.585. However, the laser intensity was likely to have changed during a run and the maximum packing in the bed measured with this method was found to vary with the fluid height. We thus decide to revisit these data and to consider another method (also based on the averaged grey-level profile) to evaluate the volume fraction. As in Aussillous et al. (2013), we suppose that the grey level intensity gradient is mainly due to a linear broadening of the laser line due to the use of a laser line generator with a finite fan angle (Dijksman et al. 2012). First, we adjust the intensity decay in the pure fluid zone by a linear function, P I = P 0l (1 − y/y 0 ), where P 0l and y 0 are constants which depend on the laser and the fluid properties, see figure 1(a,c) corresponding respectively to run 1 and run 12 of Aussillous et al. (2013). Second, for the liquid-particle mixture we assume that P I (y) = [(1 − φ)P 0l + φP 0s ](1 − y/y 0 ) where P 0s is related to the particle properties. We then introduce A = P I (y)/(1 − y/y 0 ) which is averaged over the same box size as that used for the velocity profile. If we suppose that φ(h c ) = φ c and A(h c ) = A c at the bottom of the mobile granular layer h c , i.e. the height for which particle motion ceases, the volume fraction is given by φ = φ c [A(y)−1]/(A c −1). In figure 1(b,d), we have plotted the normalised volume fraction, φ/φ c , versus the vertical position made dimensionless using the particle diameter, y/d p , for (b) borosilicate particles (Q f = 2.710 −6 m 3 /s and h f = 6.3 mm) and (d) PMMA particles (Q f = 2.710 −6 m 3 /s and h f = 15.3 mm). We observe that the volume fraction increases rapidly downward from the fluid-particle interface and reaches quickly a constant value after approximately two sphere diameters. To examine the rheological behaviour of the sheared sediment bed, we need to infer the depth-resolved shear rate,γ, the mixture shear stress, τ , and the particle pressure, p p , from the particle velocity profile, u p , and the volume fraction profile, φ. First, we interpolate linearly the velocity and volume fraction profiles to obtain their value at the fluid-bed interface, u p,in and φ in , see the red x and blue + on the dashed line in figure 1(b,d). The shear rate profile,γ(y) = du p /dy, is simply deduced numerically from the particle velocity profile using a second-order finite difference approximation. To evaluate the total shear stress and the particle pressure, we use the two-phase modelling developed in Aussillous et al. (2013). In this approach, a flat particle bed of thickness h p is subjected to a Poiseuille flow driven by a pressure gradient, ∂p f /∂x, in a horizontal channel. The flow is assumed to be two-dimensional, stationary, uniform in streamwise direction, parallel, and laminar. From the fluid phase equation (i. e. the Brinkman equation), the volume-averaged velocity in the horizontal x-direction is found to be U = φ u p + (1 − φ) u f ≈ u p ≈ u f in the bed due to the small permeability. In the laminar regime, the momentum equations for the mixture (particles plus fluid) then write where g is the gravitational acceleration and ∆ρ = ρ p − ρ f corresponds to the density difference between the two phases. To compute the particle pressure profile inside the bed, equation (2.2) is integrated numerically along the vertical y-direction. It shows that the pressure of the particle phase is proportional to the apparent weight of the solid phase and increases when penetrating inside the bed. The total stress profile inside the bed is deduced from equation (2.1), provided the stresses applied by the fluid at the fluid/bed interface, τ f (h p ), and the imposed fluid pressure gradient, ∂p f /∂x, are given. Due to the lack of measurements of the fluid velocity in the pure fluid zone for the borosilicate particles, we need to reconstruct the fluid velocity profile for this batch of spheres based on the following assumptions: (i) The fluid velocity, u f , is zero at the top wall (y = L y ) and matches the particle velocity, u p,in , at the bed height (y = h p ).
(ii) The fluid velocity profile is parabolic for y > h p which yields the analytical solution for a mixed Couette-Poiseuille-flow: where h f = L y − h p is the height of the clear fluid layer .
(iii) The fluid flux, q f,exp , is written as considering that the fluid velocity matches the particle velocity inside the bed. Using the experimental flow rate per unit width, q f,exp , and the experimental vertical profiles of φ and u p , the pure fluid flow rate q f is deduced from equation (2.4) by numerical integration and the pressure gradient for the parabolic velocity profile is then obtained Inserting this value in equation (2.3), we can calculate the fluid velocity profile in the pure fluid zone, see the black full curves in figures 1(b-d). The agreement with the experimental fluid velocity measurement for the PMMA particles is found to be good, see figure 1(d). This validates the present reconstruction method and its use for the borosilicate particles. Note that this method does not enforce a continuity of stress from the clear fluid phase to the particle phase; the obtained fluid and particle velocity profiles exhibit a change in slope at the fluid/bed interface. The fluid shear stress at the bed interface is then given by in /h f and the vertical profile of the shear rate is simply given by equation (2.1).
The extraction and reconstruction methods described above provide all the quantities needed to investigate the rheology of the mobile sediment bed: µ(J) and φ(J) or equivalently η s (φ) and η n (φ).

Simulation data
In the present work, we use the framework described in detail in Biegert et al. (2017a,b) to execute several simulations in an attempt to compare to the different experimental results of Aussillous et al. (2013) at different flow rates and fluid heights. In order to keep the paper self-contained, we provide a brief summary of the computational approach.
The particle-laden flows of interest require us to solve the Navier-Stokes equation , where p f represents the fluid pressure with the hydrostatic component subtracted out and E the identity matrix. The righthand side includes the volume forces The former is a source term used to create the pressure gradient driving the flow and the latter an immersed boundary force used to enforce the no-slip condition on the particle surface. We discretize the equations of motion for the fluid on a cubic finite difference mesh (λ = ∆x = ∆y = ∆z).
The numerical treatment is based on the IBM for fluid-particle coupling (Uhlmann 2005;Kempe & Fröhlich 2012a) and the scheme of Biegert et al. (2017a) for particleparticle interaction.
We solve for the translational velocity, u p , and the angular velocity, ω p , of spherical particles, where m p is the particle mass, I p the particle moment of inertia, V p the particle volume, and g = (0, −g, 0) T the gravitational acceleration vector. The fluid acts on the particles through the hydrodynamic stress tensor τ f , where r p represents the vector from the particle center to a point on the surface Γ p and n is the unit normal vector pointing outwards from that point. The net force and torque acting on the particle center of mass due to particle interactions are given by F i,p and T i,p , respectively. We evaluate the particle interaction forces and torques according to Biegert et al. (2017a). The resulting collision model involves normal contact forces, F n , tangential (frictional) contact forces, F t , and short-range hydrodynamic lubrication forces, F l , to provide the total collision force where N tot is the total number of particles and the subscript pq indicates interactions of particle p with particle q. According to Cox & Brenner (1967), we model the unresolved hydrodynamic component of the lubrication forces in our simulations as where ζ n is the gap size, R eff = R p R q /(R p + R q ) the effective radius, and R p and R q are the radii of the two interacting particles. Furthermore, ζ min = 3 · 10 −3 R p is a limiter as calibrated by Biegert et al. (2017a) that can be interpreted as the roughness of the particle surface and g n is the normal component of the relative velocity of the two colliding particles. The repulsive normal component is represented by a nonlinear spring-dashpot model for the normal direction where k n and d n represent stiffness and damping coefficients that are adaptively calibrated for every collision/contact to prescribe a restitution coefficient of Here, u o and u i indicate the normal components of the relative particle speed immediately after and right before the particle contact, i.e. ζ n = 0. The forces in the tangential direction are modeled by a linear spring-dashpot system capped by the Coulomb friction law as where k t and d t are stiffness and damping computed according to Thornton et al. (2011) and Thornton et al. (2013), µ f = 0.15 represents the friction coefficient between the two surfaces, ζ t is the tangential displacement integrated over the time interval for which the two particles are in contact (Thornton et al. 2013) and t is a unit vector pointing into the tangential direction. The empirical parameters e dry = 0.97 and µ f = 0.15 have been taken from experiments involving glass spheres (Gondret et al. 2002;Joseph & Hunt 2004). Using these values, the contact-model has been validated in detail by Biegert et al. (2017a) against seminal experimental benchmark data involving the same material (Foerster et al. 1994;Gondret et al. 2002;Ten Cate et al. 2004;Aussillous et al. 2013). Note that these parameters were not measured for the particles used in the experiments of Aussillous et al. (2013) and may, hence, be different.
In the present work, we consider a computational setup very similar to the one presented in Aussillous et al. (2013) We discretize the domain with a regular grid that is equidistant in all three directions and has a resolution of 25.6 grid cells per particle diameter, which is a fine enough resolution to obtain results independent of the grid cell size as shown by Biegert et al. (2018). We generate the bed by allowing 4,339 monodisperse particles to settle under gravity, without the influence of the surrounding fluid, onto a layer of 200 fixed particles whose centers randomly vary in height above the bottom wall within a range of d p , providing an irregular roughness (Jain et al. 2017). The resulting bed fills the domain to about a height of h p ≈ 20d p from the bottom wall, where h p is the particle bed height, leaving a gap of about 10d p between the top wall and Lines indicate the average streamwise (x-direction) fluid velocity (dark gray, red online) and particle velocity (light gray, yellow online), where the consecutive averaging in the horizontal plane and time is defined by (3.9).

Number of particles Ntot 4339
Galileo number Ga = ρ f η f ρp ρ f − 1 g d 3 p 0.85 ρp/ρ f 2.1 Timestep CFL = 0.5 Domain size (Lx/dp × Ly/dp × Lz/dp) 20 × 30 × 10 Domain grid size (Lx/λ × Ly/λ × Lz/λ) 512 × 768 × 256 Domain boundary conditions periodic × no-slip × periodic Initial h f /dp 10.0 Particle resolution, dp/λ 25.6 We enforce different volumetric flow rates, governed by the volume force f b,x . As it can take a long time for a simulation to reach a statistically steady state when initialized from rest, we found it to be more efficient to obtain a steady state by initializing the   ). After this initialization phase, the imposed pressure gradient is reduced to produce simulations Re17 and Re33. Re8 is carried out by continuing Re17 with an even lower imposed pressure gradient. This procedure allows us to quickly reach a steady state for runs Re17 and Re8, as can be seen in figure 3a, whereas Re33 is still in the process of dilation.
To compare the simulation results to rheological models, we have to analyze the discrete information of the particle-phase in our simulations, such as particle velocities and forces, from a continuum viewpoint. To this end, we employ the coarse-graining method (CGM) based on the works of Goldhirsch (2010) and Weinhart et al. (2012). This coarse-graining method conserves quantities of interest, but additionally smoothes out the resulting continuum field. For a given discrete particle quantity θ p , we define its coarse-grained continuous counterpart, θ cg , as where x p (t) is the position of the center of particle p, and W(r) is the conservative coarse-graining function that smears a given local quantity in a spherical volume of radius |r|.
The main property is that R 3 W(r) dr = 1. We implemented a coarse-graining function based on the Dirac delta function of Roma et al. (1999), which is the same delta-function used in the IBM (Uhlmann 2005). The filter |r| has to be as small as possible to fully exploit the highly-resolved simulation data for our rheological analysis, but large enough to smooth out the sub-particle scale. Here, we choose |r| = 1.5d p , which was deemed to be a good compromise to satisfy these requirements.
The steady-state configuration for the moving bed is steady only in a time-averaged sense, because particle collisions and positions continuously fluctuate. We therefore define the time average of a given continuous field θ to be where the overbar and angular brackets represent averaging in time and space, respectively (Vowinckel et al. 2014. Note that this averaging operator applies for both, continuous fluid quantities and coarse-grained particle quantities. We present the values for t avg,1 and t avg,2 in table 2. These time-averaging windows were chosen to capture the steady-state results, or as large a time span as possible for as similar a particle flux as possible (figure 3). We will show, however, that our results for the rheology are independent of transient behavior for all Reynolds numbers investigated.
We can use (3.8) to obtain a continuous field of φ and (3.9) to generate vertical profiles of φ and u f by computing double-averaged values (in space and time, cf. figure 4). For the analysis of the rheology, we exclude values for y/y ref < 0.68 to eliminate boundary effects from the artificial bed roughness at the bottom all. We observe the same evolution for the volume fraction as in the experiments with a rapid increase downward from the fluid-particle interface to reach a constant value after a few sphere diameters. Here, we follow the definition of Biegert et al. (2017a) to determine the fluid-sediment interface at φ = 0.05. The oscillations of φ in this region reflect the sub-particle scale as the particle layering within the sediment bed (figure 4a). We, therefore apply the coarsegraining radius of |r| = 1.5d p to the grid resolved data shown in this figure to smooth out the layered structure in these vertical profiles. Note that figure 4b yieldsγ = ∂ u f /∂y as another rheological quantity. Looking at the double-averaged fluid velocity profiles in figure 4b, the three cases represent three distinctively different regimes (Jenkins & Larcher 2017). The sediment bed of Re8 reaches a quasi-static regime of the granular suspension, whereas the sediment motion in Re17 and Re33 can be considered "layered" and "collisional", respectively. There is a clear qualitative difference between run Re8, whose velocity profile is concave and goes to zero within the bed at y/y ref ≈ 0.5, and run Re33, whose velocity profile is convex and goes to zero only at the fixed particles at the lower wall. Note that the fluid velocity is to a good approximation equal to the particle velocity, which is consistent with the observation of Aussillous et al. (2013). The Stokes number ranges from 0.19 St 0.77, which is 20 to 80 times larger than the values obtained from the experimental data §2, but still expected to be in the viscous regime (Ness & Sun 2015).
While we can readily extract the quantities for φ andγ from figure 4, an additional investigation of the total stress balance of the fluid-particle mixture in x-and y-direction is needed to compute the total shear stress τ and the granular pressure p p as will be detailed in the next section.

Stress balance of the simulation data
This section presents the analysis of the stress balance for the fluid-particle mixture to compute wall-normal profiles of the rheological quantities τ and p p . To this end, we follow the argument of the previous numerical work Biegert 2018;Vowinckel et al. 2019), where the full derivation of the stress balances in both shearing (x) and wall-normal (y) direction are given from first principles. In addition, a detailed analysis of the present datasets for all components entering the stress balances in shear and wall-normal direction, respectively, can be found in Biegert (2018).

Stress balance in the x-direction to obtain the total stress
To obtain the total shear stress, we consider the momentum balance of the fluid phase, i.e. the Navier-Stokes equations (3.1), over the control volume Ω CV in the x-direction spanning from the top wall Γ w to some arbitrary height y with the lower boundary being Γ y and involving multiple particles (figure 5). We distinguish the free fluid from the particle interior as exemplified in figure 5 using Ω + CV and Γ + y for parts of the domain that are occupied by fluid and the lower boundary, respectively.
We can simplify the momentum balance of the fluid phase for the following reasons. Owing to the periodic boundary conditions and the fact that we impose the Poiseuille flow via the volume force f b , the pressure term (∂p f /∂x in (3.1)) does not contribute to the x-momentum, since our control volume covers the entire streamwise extent of the domain. At the top wall Γ w , the vertical velocity, v f , is zero, so that only η f ∂u f /∂y contributes to the fluid stress. At the lower boundary, Γ y , the pressure again does not play a role and the convective term vanishes due to the laminar flow conditions, so that the long-range hydrodynamic stress originates from the viscous term. Finally, we include f IBM,x as the sum of all hydrodynamic stresses arising from pressure and viscous forces that act on the particle surfaces L CV enclosed in Ω CV . These assumptions yield (4.1) The last term Particle force provides the linkage to the momentum balance of the particle phase, which can be obtained by appling the coarse-graining method (3.8) to the particle equation of motion (3.2). This yields where a cg , F cg h , F cg g F cg i are the coarse grained forces due to the particle acceleration, hydrodynamic stress due to the IBM, gravity and particle interaction, respectively. Similar to the fluid momentum balance, we can analyze the coarse-grained particle forces within a control volume spanning the entire domain in the streamwise and spanwise directions and extending from the top wall to an arbitrary height y. Integrating (4.2) over this volume, we obtain If particles are in a steady state, either naturally or through double-averaging, then the acceleration term vanishes. Furthermore, our setup of a horizontal channel yields F g,x = 0. We can, therefore, use the fact that f IBM,x = F cg h,x = −F cg i,x and apply (3.4) to further distinguish between normal and tangential forces due to particle contact and friction (F cg n,x and F cg t,x ) and short-range hydrodynamic lubrication forces (F cg l,x ). Using the averaging operator (3.9) and dividing by the horizontal area of the domain, we can rewrite (4.1) as where γ is an indicator function for the fluid (γ = 1 outside the particle and γ = 0 inside the particle), resulting in double-averaged equations for laminar flows akin to Nikora et al. (2013) and Vowinckel et al. (2017a,b). We have also used the fact that η f , ρ f , and f b,x are constant throughout the domain. It is important to note that we are explicitly separating the stresses arising from hydrodynamic interactions and particle contacts, respectively (Gallier et al. 2014(Gallier et al. , 2018, which will be analyzed in more detail in §5. The right-hand side of (4.4) comprises the long-range fluid stress due to viscous effects evaluated at height y at the boundaries of the control volumes Γ + CV occupied by fluid as well as the short-range lubrication effects within the control volume Γ + CV . The contact stress comprises both, normal contact forces and tangential frictional forces. Note that it was shown in Biegert et al. (2018) and Biegert (2018) that there is also a viscous and convective stress inside the particles Γ − CV as a by-product of the IBM. For the present study, this effect does not contribute a significant part to the total stress balance. The external stress on the left-hand side of (4.4) consists of the viscous stress at the top wall and the stress from the body force acting throughout the control volume. It is also equal to all other stresses arising from the movement of the fluid and particles and is hence equivalent to the total stress τ f needed to compute the rheological quantities µ(J) and η s .
For the analysis of the stress balance of the fluid phase, we focus on runs Re8 and Re33 to get a sense of the results for different flow conditions. Figure 6 shows the momentum balance of the fluid phase, given by (4.4), for runs Re8 (figure 6a) and Re33 (figure 6b), in which we expect the external stress to match the sum of the hydrodynamic and contact stresses. As expected, the external stress at the top wall is close to σ x /σ ref = −1. In the upper part of the flow (y/y ref > 2.15 and y/y ref > 2.3 for Re8 and Re33, respectively), there are no particles, and the hydrodynamic stress matches all of the external stress. Within the particle bed (y/y ref < 2.3), however, the majority of the external stress is taken up by particle contact. The total stress comes out to be a linear profile. Hence, the present results support the conceptual model (2.1) proposed by Aussillous et al. (2013) (figure 4b in this reference) and we can use these data to compute η s and µ in the following.
Naturally, the hydrodynamic stress of the laminar flow is entirely made up of the viscous term alone. Run Re8 differs from Re33 in that the hydrodynamic stress reaches a higher positive value above the particle bed and quickly drops to zero within the particle bed. The hydrodynamic stress for Re33, on the other hand, increases with increasing depth within the particle bed. Since the entire sediment bed is set in motion for this case, the lubrication component makes a significant contribution to the hydrodynamic stress. These results are consistent with the velocity profiles in figure 4b, where the concavity of the profile for Re8 results in a high shear stress at the fluid/particle bed interface and low stresses within the bed, while the convexity of the profile for Re33 leads to a large shear stress at the lower wall. The total stress, on the other hand, is completely dominated by particle contact within the sediment bed. This result is consistent with the locations of the sharp gradients in the stress balance, so that the hydrodynamic and contact stresses together close the x-momentum balance. It also shows that effects from the local acceleration term are negligible and that the unsteady flow conditions (figure 3) are expected to have no impact on the reported results even for cases of dilating and consolidating sediment beds (e.g. Re33).

4.2.
Stress balance in the y-direction to obtain particle pressure The rheological quantities η n and µ(J) require information about the particle pressure p p , which we can obtain by further analyzing the particle phase. Another way to interpret the bed particle pressure p p is to think of it as the total submerged weight of the particles (Vowinckel et al. 2019). Indeed, this has been done by Stickel & Powell (2005) and Ouriemi et al. (2009) for continuum modeling. We can demonstrate the validity of this reasoning by analyzing the momentum balance for the particle phase in the y-direction. To this end, we utilize the coarse-grained momentum balance of the particle phase (4.2) and apply the averaging operator (3.9) to recast (4.3) as a line integral in the wall-normal direction (4.5) Since we are interested in the granular pressure, i.e. the bed weight, we consider only the vertical y-component. We again neglect the acceleration term, but keep the gravity term that is acting in the same vertical direction. Furthermore, we subdivide F cg i into the short-range hydrodynamic lubrication and particle contact component, (4.6) Since F g,p = V p (ρ p − ρ f )g, coarse-graining the particle weight yields the left-hand side as the submerged bed weight, which is equivalent to the particle pressure: We also note that (4.7) is the integral of (2.2), so that the present derivation provides a direct linkage to the two-phase modeling of Aussillous et al. (2013). Owing to the fact that we obtain full information of the solid volume fraction (cf. figure 4b), we do not need to introduce an artificial pressure at the fluid sediment interface. This was suggested by Houssais et al. (2016) (called P 0 in that study), but using such an artificial pressure would neither be in line with the two-phase modeling, nor would it obey the definition of our control volume sketched in figure 5. Hence, omitting P 0 naturally yields J → ∞ for φ → 0, because p p → 0. This behavior is consistent with the approach of Boyer et al. (2011). Figure 7 shows the coarse-grained particle phase stresses, given by the time average of (4.6) for runs Re8 (left column) and Re33 (right column). In figures 7a and 7b, the bed weight increases almost linearly from the top of the particle bed down to the lower wall, balanced by the sum of the hydrodynamic and collision stresses. Again, this observation confirms the conceptual model of (2.2) of Aussillous et al. (2013). The fact that we have found a linear profile for this physical quantity simplifies the situation from a modeling perspective, i.e. one needs the sediment height and the total submerged weight of the sediment bed to reconstruct depth-resolved profiles of p p . The results for the stress balance in y-direction show clear differences between runs Re8 and Re33. Owing to the normalisation, the dimensionless y-momentum collision stress is three times larger for Re8 than for Re33. Moreover, the hydrodynamic stress is positive for Re8 and negative for Re33, so that the collision stresses are less than and greater than the bed weight, respectively, for the stress balance to be in equilibrium. These deviations from zero stem from the fact that Re8 and Re33 are in the process of consolidation and dilation, respectively.

Rheology
We now turn to the examination of the rheological quantities extracted from the experimental measurements in § 2 and from the numerical simulations in § 3 and 4. We display these data within the pressure-imposed view by plotting the effective friction coefficient, µ = τ /p p , and the bulk volume fraction, φ, as a function of the viscous number, J = η fγ /p p , in figure 8. The macroscopic friction coefficient, µ, is found to be an increasing function of the viscous number, J, whereas the volume fraction, φ, is a decreasing function of the viscous number, J, for both, the experimental data coming from the two batches of spheres made of borosilicate (+) and PMMA (×), where no appreciable differences can be seen between the runs of these two particle materials, and the numerical data coming from the three different runs at different Reynolds numbers.
It is important to note, however, that figure 8 is not meant as a validation of the numerical results by comparing it to the post-processed experimental data of Aussillous et al. (2013) for the following reasons. First, the two data sets represent different conditions. The particle properties e dry and µ f chosen for the simulations may not correspond to the (unknown) values of the experiments and the flow conditions expressed by the Stokes number differ by more than one order of magnitude. Second, owing to the difficulty to capture experimentally the sediment transport layer close to the bed interface of only a few particle diameters thickness, the experimental data present significant scatter. Finally, recall that our data processing in §2 assumes u f = u p,in to reconstruct the experimental velocity profile, which slightly underestimates the fluid velocity in the sediment transport layer. On a similar note, the averaging operator using coarse-grained quantities (3.8) for the simulations yields an averaging window of three diameters in height to smooth out the sub-particle scale. For all these reasons, we can observe a rather large discrepancy between the experimental and numerical data for large J, i.e. in the sediment transport layer, but the qualitative trend is correctly reproduced.
Despite the discrepancy at large J, both data sets show asymptotic behavior for small J in figure 8b. The results can therefore be used to deduce the critical parameters µ c and φ c and provide proper scaling for our rheological analysis. Among the numerical results, only the data for the smallest Reynolds number (run Re8) reach the limit of small J, i. e. J < 10 −2 and, hence, smallγ. The data coming from run Re8 and from the experiments show that both φ and µ tend to finite measurable values, φ c and µ c respectively, at vanishing J, i. e. at the jamming point where the suspension ceases to flow. The critical (or maximum flowable) volume fraction, φ c , and friction coefficient, µ c , can be measured by fitting the data using an asymptotic linear regression in √ J for small J (for J < 10 −2 ) as done by Tapia et al. (2019). The critical volume fraction is found to be φ c ≈ 0.59 for both the experiments and simulations and is similar to the value found for suspensions made of spheres having only small surface roughnesses or equivalently having moderate inter-particle sliding friction coefficients (Boyer et al. 2011;Dagois-Bohy et al. 2015;Tapia et al. 2019). The critical friction coefficient is found to be µ c ≈ 0.18 for the experiments and µ c ≈ 0.27 for the simulations. These values are smaller than those found previously (µ c ≈ 0.30 − 0.37) in pressure-imposed rheometry (Boyer et al. 2011;Dagois-Bohy et al. 2015;Tapia et al. 2019) but are of the same order of magnitude or even slightly larger than that given (µ c ≈ 0.17) by the correlation of Morris & Boulay (1999) which was meant to match experimental results on shear-induced migration.
The critical values φ c and µ c are not universal parameters but depend on particle properties, i. e. the values of φ c and µ c depend on the particle size distribution but also on their surface properties. It is thus convenient to plot the data of figure 8 by scaling φ by φ c and µ by µ c for a comparison with pressure-imposed rheological measurements (Boyer et al. 2011;Dagois-Bohy et al. 2015;Tapia et al. 2019) and correlations (Boyer et al. 2011;Morris & Boulay 1999). These scaled data shown in figure 9 collapse reasonably well onto the same constitutive curves for the small J-range (J 10 −2 ) but present discrepancies at larger J, i. e. low φ. However, the discrepancy at low φ is not surprising, because the empirical correlations of Boyer et al. (2011) and Morris & Boulay (1999) have been derived by fitting to experimental data of volume fractions φ/φ c > 0.5 for neutrally buoyant particles that are homogeneously distributed in a flow cell. In fact, it is worth noting that there even exist some disparities within the pressure-imposed measurements, as while the data of Dagois-Bohy et al. (2015) and Tapia et al. (2019) are rather similar, they differ from the data of Boyer et al. (2011) at large J. This is reflected in the correlation of Boyer et al. (2011) which is designed to agree with the experimental measurements of Boyer et al. (2011). This correlation seems to provide a lower bound for the whole collection of data while that of Morris & Boulay (1999) is more like an upper bound for the data of µ(J) reported in figure 9a. The plotted data, therefore illustrates the range of uncertainty related to the extrapolation of empirical relations. It is important to stress that these two correlations have different critical values: φ c = 0.585 and µ c = 0.32 for the correlation of Boyer et al. (2011) and φ c = 0.68 and µ c = 0.17 for the correlation of Morris & Boulay (1999). Note that value of φ c = 0.68 is rather unrealistic as it exceeds the randomly closed packing fraction. We also provide the classical volume-imposed view by plotting the shear, η s = τ /η fγ = µ/J, and normal, η n = p p /η fγ = 1/J, viscosities as a function of the scaled volume fraction, φ/φ c , in figure 10. The collection of data shows that both viscosities increase with increasing φ and diverge at φ c . The data coming from the erosion experiments follow this general trend despite the large scatters. Again, the correlations provide bounds for the whole collection of data, but this time the correlation of Boyer et al. (2011) acts as an upper bound while that of Morris & Boulay (1999) as a lower bound. The simulation data are notably lower than the pressure-imposed rheological measurements of Boyer et al. (2011), Dagois-Bohy et al. (2015), and Tapia et al. (2019) in particular for the smaller φ-range.
The log-log plots shown in the insets of figure 10 are made to analyze the asymptotic behaviors close to the jamming transition. The data for η n coming from run Re8 in the inset of figure 10 (b) present a divergence in (1 − φ/φ c ) −2 in agreement with previous work (see e. g. Boyer et al. 2011;Tapia et al. 2019); the other runs are too far from the jamming point to allow for any conclusions. The log-log plot shown in the inset of figure 10 (a) reveals that the same divergence dominates close to jamming for η s . This last point is also demonstrated by the finite value of µ = η s /η n at vanishing J for the data coming from run Re8 and from the erosion experiments in figure 8 (a). Apart from the difficult task to accurately measure within the sediment transport layer of only a couple particle diameters thickness, the issues described above to recover the rheological quantities at large J and small φ can be attributed to the fact that sediment transport can be sub-divided into two different regimens (Revil-Baudard et al. 2015). On the one hand, the dynamics of dense suspensions with large φ are dominated by frictional contact. In fact, this is the regime that has been investigated in the experiments of Boyer et al. (2011), Dagois-Bohy et al. (2015), and Tapia et al. (2019) and for which the empirical correlations shown in figures 9 and 10 have been derived. On the other hand, low volume fractions represent the dilute regime of mostly binary collisions between particles. It was therefore concluded by Maurin et al. (2016) that those empirical correlations may need adjustments for these flow conditions.
Hence, our observations are in line with the study of Houssais et al. (2016), who conducted experiments of sediment transported in an annular flume and measured depthresolved profiles very similar to the data presented here. In this study, the authors had to introduce a confinement pressure to correct µ for large values of J. This pressure was added to the granular pressure and can be interpreted as the artificial weight of the top plate in a pressure-imposed rheometer. Houssais et al. (2016) found good agreement with the correlations of Boyer et al. (2011) for µ(J) using this measure, but unfortunately, the agreement for the correlation φ(J) was not as good in that study. As already mentioned in §4.2, we omit using this strategy of an additional confinement pressure in the present analysis. The high resolution of the numerical data did not require such a treatment, because all relevant quantities could be recovered by our statistical analysis in a straightforward manner. This way, we ensure to stay consistent with the two-phase flow model (2.1) and (2.2) from Aussillous et al. (2013). The same applies for the revisited data of Aussillous et al. (2013), but for these data we use the additional assumption of u f = u p,in to reconstruct the fluid velocity profile from a mixed Couette-Poiseuille-flow, which yields a slight underestimation of the fluid velocity in the range of large J and low φ.
To further investigate the governing effects in the dilute and dense regime as a function of the particle volume fraction, we make use of the highly-resolved data and the stress balance (4.4) derived in §4 and separate the stresses due contact and hydrodynamic interaction. These are plotted in figure 11 (without the normalization by φ c ) together with a comparison of the data of Gallier et al. (2014), who carried out three-dimensional simulations of neutrally buoyant, non-Brownian, frictional spheres in a Couette cell of constant size, i.e. volume-imposed rheometry. For our simulation results, the data of the three simulation runs collapse onto a single master curve for both the contact and the hydrodynamic stress components. For φ > 0.3, there are no significant longrange hydrodynamic effects possible and the particle contact is the main contribution to the total stress, whereas hydrodynamic effects prevail below this threshold value of φ. The analysis therefore reveals the reason why the considerations for pressure-imposed rheometry also hold for sediment transport in the dense regime, i.e. φ > 0.3, where the rheology is dominated by particle contact. In the dilute regime, i.e. φ < 0.3, the hydrodynamic component scales very well with the Einstein relation η s = 1 + 5 2 φ (Einstein 1956) and the contribution from particle contact becomes negligible. Our analysis furthermore shows that the deviation of the hydrodynamic component from its clear fluid value (η s = 1) is induced by the lubrication forces (not shown here). This observation differs from the simulation results of Gallier et al. (2014) obtained using volume-imposed rheometry for 0.1 φ 0.45, where particles are distributed homogeneously in the simulation domain. In the present simulations, the dilute regime represents the thin layer of active sediment transport, where particles are still touching and aligned in the shearing direction so that the particle volume fraction decreases rapidly from densely packed to zero within a thin layer that is only a couple of diameters thick. Hence, long-range hydrodynamic effects are screened by the highly anisotropic distribution and the steep gradient of φ in this layer. As a result, the hydrodynamic component resembles a porous media flow behavior as anticipated in the two-phase modelling of Ouriemi et al. (2009) and Aussillous et al. (2013).
To the knowledge of the authors, the only numerical study of particle-resolved DNS that was able to perform the analysis of the rheological behavior of a sheared sediment bed has been presented by Kidanemariam (2017). Interestingly, this study underestimates the empirical correlation for φ(J)/φ c for large J, whereas our simulation data overestimates the volume fraction in this regime (cf. figure 9 (b)). The low values of φ for large J in the data of Kidanemariam (2017) can in parts be attributed to the difference in handling particle collisions and contact. Kidanemariam (2017) did not resolve lubrication forces which imposes the constraint that one has to keep a minimal distance of two grid cells between particles. Consequently, the granular packing has to be less dense so that unusually low values of 0.43 φ c 0.53 were reported for the maximum packing fraction. Hence, the data and the analysis presented here demonstrate a substantial improvement over previous efforts to capture the rheology of sediment beds over a wide range of J.
The rheological data of figures 8, 9, 10 and 11 are given as supplementary materials.

Concluding remarks
Serving the need to use constitutive laws for macroscopic sediment transport models, the present paper introduced a new means to generate highly resolved data for the rheological behavior of granular sediment beds sheared by a viscous, pressure-driven flow. The rheology is assessed by a statistical analysis to extract the physical quantities needed to compute the rheological parameters. The highly-resolved simulations provide the data needed to assess the stress exchange of the fluid-particle mixture by deriving the stress balance from first principles. To better understand naturally occurring sedimentladen flows in pipes and rivers, we focus on pressure-driven flows, where the total stress increases with flow depth. The analysis verified the conceptual two-phase model of Ouriemi et al. (2009) and Aussillous et al. (2013) that has been derived from the Brinkman equations, and deduces the total fluid shear and the granular pressure as the relevant rheological quantities. We compare our numerical results to the rheology of corresponding experiments by revisiting the data set of pressure-driven flows investigated by Aussillous et al. (2013) and found reasonably good agreement between the numerical and experimental data. The results presented here clearly show that sediment transport by pressure-driven flows yields results that are similar to those of previous studies of annular Couette-type flows, even though these studies were conducted with dense suspensions of neutrally buoyant particles in pressure-imposed rheometers. The analyzed data agree well with the empirical correlations of Boyer et al. (2011) and Morris & Boulay (1999) derived therefrom. In the more dilute regime, the simulation data fall in between the bounds set by those two extrapolated correlations. The good agreement of the data obtained from pressure driven flows with those from Couette flows potentially justifies the use of these empirical correlations as constitutive equations for two-phase flow solvers for sediment transport applications in the dense regime. Separating the total shear stress of our simulation results into the parts from hydrodynamic interaction and particle contact revealed a critical particle volume fraction of φ ≈ 0.3, for which particle contact becomes the dominant component, so that classical empirical relations from rheometry become applicable. In the dilute regime, we obtain a scaling in agreement with the Einstein relation, which illustrates a flow behavior that differs from classical rheometry studies, because the long-range interactions are screened by the porous media. This effect corresponds to the thin layer of sediment transport over which hydrodynamic and contact stresses but also φ vary rapidly, which is very difficult to grasp, but it constitutes the main difference between classical rheometry of neutrally buoyant suspensions and the rheology of sediment transport. More work of this kind will be needed in the future to address different particle properties and flow situations that are more complex than the two fundamental flow types Couette and Poiseuille flow.