Rheology of mobile sediment beds in laminar shear flow: effects of creep and polydispersity

Classical scaling relationships for rheological quantities such as the $\mu(J)$-rheology have become increasingly popular for closures of two-phase flow modeling. However, these frameworks have been derived for monodisperse particles. We aim to extend these considerations to sediment transport modeling by using a more realistic sediment composition. We investigate the rheological behavior of sheared sediment beds composed of polydisperse spherical particles in a laminar Couette-type shear flow. The sediment beds consist of particles with a diameter size ratio of up to ten, which corresponds to grains ranging from fine to coarse sand. The data was generated using fully coupled, grain resolved direct numerical simulations using a combined lattice Boltzmann - discrete element method. These highly-resolved data yield detailed depth-resolved profiles of the relevant physical quantities that determine the rheology, i.e., the local shear rate of the fluid, particle volume fraction, total shear, and granular pressure. A comparison against experimental data shows excellent agreement for the monodisperse case. We improve upon the parameterization of the $\mu(J)$-rheology by expressing its empirically derived parameters as a function of the maximum particle volume fraction. Furthermore, we extend these considerations by exploring the creeping regime for viscous numbers much lower than used by previous studies to calibrate these correlations. Considering the low viscous numbers of our data, we found that the friction coefficient governing the quasi-static state in the creeping regime tends to a finite value for vanishing shear, which decreases the critical friction coefficient by a factor of three for all cases investigated.


Introduction
The fluid mediated transport of granular sediment is a key process for the mass movement in a geophysical but also an engineering context (e.g. Frey & Church 2011). The transport typically occurs along a slope or by a fluid flow shearing the sediment (Jerolmack & Daniels 2019) and can lead to bedform evolution, such as ripples and dunes, even for laminar flow conditions (Lajeunesse et al. 2010). This consideration allows to characterize sediment transport in laminar flows in terms of the rheology to investigate the fluid-particle mixture's deformation behavior in shearing flows (Aussillous et al. 2013;Houssais et al. 2016;Kidanemariam 2016;Vowinckel et al. 2021). All these studies justified their approach by comparing the results to data previously obtained in rheometer studies with dense suspensions of neutrally buoyant particles (e.g. Morris & Boulay 1999;Boyer et al. 2011). For these classical rheological investigations, a shear rate is applied to a dense granular material suspended in a fluid with viscosity to investigate the total shear stress acting on the fluid-particle mixture in the shearing direction and the imposed particle pressure in the wall-normal direction. The total shear comprises hydrodynamic and frictional inter-particle stresses, with the latter becoming more important with increasing particle volume fraction (Gallier et al. 2014;Guazzelli & Pouliquen 2018;Vowinckel et al. 2021).
In this regard, two types of rheometer setups are possible. On the one hand, the volumeimposed rheometry confines the suspension by shearing walls with constant gap size (e.g. Morris & Boulay 1999). While Morris & Boulay (1999) were investigating shear induced migration to begin with, they were also able to measure the effective shear and normal viscosities, = / and = / , respectively, and to derive empirical correlations for these two quantities as functions of . On the other hand, a pressure-imposed rheometer, where a constant confining pressure is applied to a movable upper wall, allows for the dilation of the dense suspension under shear (e.g. Boyer et al. 2011;Dagois-Bohy et al. 2015). For laminar viscous flows, i.e. a Stokes number = 2 / smaller than 10 (Bagnold 1954; Ness & Sun 2016), where is the particle density and is the characteristic particle diameter, this measure allowed Boyer et al. (2011) to define a macroscopic friction coefficient = / that depends on the viscous number = / . Based on this, the authors were able to propose empirical correlations for ( ) and ( ) that distinguish between stress contributions from particle contact and hydrodynamic interactions. This framework has become known as the ( )-rheology. In this article, we will follow the nomenclature of Guazzelli & Pouliquen (2018) and use the symbol rather than for the viscous number to distinguish it more clearly from the inertial number defined for highly inertial granular flows.
The pressure-imposed rheometry also allows for the analogy to sediment transport, where the imposed particle pressure at some depth in the sediment bed is equal to the submerged weight of the overlying grains (Aussillous et al. 2013;Maurin et al. 2016;Vowinckel et al. 2019a). This analogy is important for two-phase fluid sediment transport modeling (Jenkins & Hanes 1998;Hsu et al. 2004), where the fluid-particle mixture is treated as two separated continua with interconnected conservation laws of mass and momentum (Ouriemi et al. 2009). The empirical correlations of the ( )-rheology can provide the constitutive equations needed to close this set of equations (Chauchat et al. 2017;Lee & Huang 2018;Lee 2021). Unfortunately, the empirical correlations ( ) and ( ) involve parameters that are not universal but were calibrated against the experimental data of Boyer et al. (2011) in the dense regime with non-vanishing shear (0.4 < < 0.58 and > 10 −6 ). It has been pointed out by Revil-Baudard et al. (2015) who investigated sheet-flow processes under turbulent flow conditions that these correlations need adjustments for more dilute systems, whereas Houssais et al. (2016) investigated viscous numbers as low as ≈ 10 −9 and found that the grains were still moving under creeping conditions even for these extremely low shear rates. It remained unclear, however, if this was a particle property or an effect originating from the curvature of the annual flume employed in this study. Hence, for cases, where the modeled flow conditions exceed the range of the calibration data, the ( )-rheology can even lead to ill-posed problems as reported by Barker et al. (2015), who then proposed an extension to tackle this problem (Barker & Gray 2017).
To increase the robustness of the ( )-rheology for two-phase fluid models, more work is needed to derive more universal constitutive equations (Denn & Morris 2014;Pähtz et al. 2019). A good starting point will be to address the coefficients that enter the models of the ( )-rheology and are known to depend on the particle properties. For the critical state of very low shear rates and dense systems, i.e. low and large , the frictional inter-particle forces may become large enough to inhibit grains sliding past one another. This quasi-static regime is determined by the particle properties critical friction coefficient 1 and maximum particle volume fraction . For example, Boyer et al. (2011) reported 1 = 0.32 and = 0.585 for the monodisperse case, but it has been shown by Tapia et al. (2019) for pressure-imposed rheometry that these two quantities decrease with increasing particle roughness. For obvious reasons, the critical volume fraction may also depend on the grain size distribution of the sediment as smaller particles can fill the void interstitial pore space provided in between larger grains (Guazzelli & Pouliquen 2018). This aspect has thus far been neglected in the framework of the ( )-rheology. In fact, most of the studies use sediment compositions of uniform grains, where the standard deviation of the grain size distribution is smaller than 10%. However, neither is this variance in grain size distribution large enough to see appreciable effects of polydispersity on the sediment transport (Biegert et al. 2017), nor does this variance reflect the grain size distribution of fluvial sediments.
In this regard, it is important to acknowledge that natural sediments are by no means monodisperse or bidisperse, but obey a certain continuous grain size distribution. For example, according to ISO 14688-1:2002, cohesionless sand grains can range from 0.063 to 2 millimeters in diameter. This calls for an extension of the ( )-rheology towards more realistic polydisperse sediment compositions.
As a first step, bidisperse suspensions were investigated in volume-imposed rheometers. For this scenario, the effective viscosities were reduced as compared to the monodisperse case (Chang & Powell 1994;Gondret & Petit 1997). In these studies, the non-uniformity of the bidisperse grains was up to ,max / ,min = 13.75, where ,max and ,min are the maximum and minimum diameter of the grains, respectively. The critical volume fraction that indicates the quasi-static regime was also increased from = 0.585 for the monodisperse case (Boyer et al. 2011) to = 0.64. Consequently, models for in bi-disperse volume-imposed rheometry were proposed by Dörr et al. (2013) and Mwasame et al. (2016) that can also be applied to polydisperse systems (Pednekar et al. 2018).
As a next step, 2D-DEM simulations with grains of continuous polydispersity have been carried out where the fluid drag was approximated by Stokes drag and lubrication (Trulsson et al. 2012;Ness & Sun 2016) and the variation of the grain size was kept constant at ,max / ,min = 3.0 and 1.4, respectively. A recent study by Amarsid et al. (2017) extended these considerations to a lattice Boltzmann -discrete element method for simulations in 2D for , / , = 1.67. Since, however, the focus of these studies was to investigate the transition from the viscous to the inertial regime, polydispersity was merely added to prevent artificial crystallization of the densely packed scenario and its role on the rheology was not discussed. To the knowledge of the authors, 3D-simulations with a systematic focus on the degree of polydispersity in pressureimposed rheometry or even sheared sediment beds have not been considered yet. The present study addresses this issue.
We employ the open-source simulation framework LB (Bauer et al. 2020a) to carry out fully-coupled particle-resolved direct numerical simulations of sediment beds sheared by a laminar Couette-type flow in the viscous regime, i.e. < 10. To this end, we utilize the combined lattice Boltzmann -discrete element method of Rettinger & Rüde (2017) and Rettinger & Rüde (2020). This extends our pore-resolved simulations of fluid flow through porous media (Fattahi et al. 2016;Gil et al. 2017;Rybak et al. 2020), and is in line with previous erosion studies using a similar methodology (Derksen 2011;Rettinger et al. 2017). We follow the approach by Vowinckel et al. (2021) to compute time-averaged, depth-resolved profiles to quantify the stress exchange between the fluid and the particle phase. This allows for a systematic simulation campaign of different sediment grain size compositions under exact control of the flow conditions and eradicates potentially unwanted effects from curved sidewalls, as present in existing laboratory experiments. The highly-resolved data yields all the relevant quantities, i.e. particle volume fraction, shear rate, total shear, and granular pressure, to infer the rheology of the polydisperse fluid-particle mixture down to viscous numbers of ≈ 10 −9 . The investigated sediment beds have a non-uniformity of ,max / ,min up to a factor of ten, which corresponds to a variety typically encountered in fluvial sediments of lowland rivers (e.g. Kuhnle 1993; Frings 2008).
The rather large disparity of the grain sizes is achieved using the efficient parallelization scheme of Eibl & Rüde (2018). These studies ultimately allow us to derive a robust parameterization strategy of the classical ( )-rheology to account for the sediment polydispersity by linking the non-uniformity to the critical volume fraction and propose a straightforward extension to creeping flow conditions that recovers the original ( )-rheology for higher shear rates.
The paper is structured as follows. We first provide a brief summary of the numerical framework in §2 and the simulation setup in §3. We then infer the pressure-imposed rheology and validate our simulation approach in §4 by comparing the monodisperse case to the experimental data of Boyer et al. (2011) andHoussais et al. (2016), including the classical empirical correlations of the ( )-rheology (Boyer et al. 2011). Finally, we utilize the data from our simulation campaign to present extensions of the ( )-rheology for polydispersity and creeping flow in §5 and §6, respectively.

Numerical Method
For the numerical studies presented here, we couple the lattice Boltzmann method for fluid flow with a discrete element method to account for particle interactions of polydisperse, spherical grains. This approach has proven to be accurate and efficient for geometrically fully-resolved particle flow simulations and has been thoroughly validated in Rettinger & Rüde (2020). Therein, a detailed presentation and discussion of the method is given. We briefly summarize the key aspects for completeness here. All parts of the employed numerical scheme are contained in the opensource high-performance framework LB (cf. Bauer et al. 2020a), and its implementation can be found in the official software repository †. A sketch of the numerical scheme is presented in figure 1.

Lattice Boltzmann method
The lattice Boltzmann method (LBM) is a relatively recent approach for the simulation of viscous fluid flow. It describes the evolution of particle distribution functions (PDFs) on a † https://walberla.net/ uniform computational grid and thereby fulfills the macroscopic Navier-Stokes equations. A detailed overview of the theory and various approaches can be found in Krüger et al. (2017). For the present studies, we employ the 3 19 two-relaxation-time model of Ginzburg et al. (2008). The relaxation times, connected via the parameter Λ = 3/16, determine the kinematic fluid viscosity and allow for accurate flow simulations. The local fluid pressure and velocity are obtained via zeroth-and first-order moments of the PDFs in a fluid cell. Commonly, all quantities are expressed in a normalized LBM unit system, the so-called lattice units, which results in the cell size Δ = 1, the time step size Δ = 1, and a reference fluid density of = 1. Those will be used in the remainder of this work.

Discrete element method
The motion of a spherical particle can be described by the Newton-Euler equations Here, , = , is the mass of the particle of density and volume , , and , = ( , 2 , )/10 is the moment of inertia for a sphere of diameter , . The temporal change of the particle's translational velocity is thus given by the acting forces , , with contributions from the collisions col , , the hydrodynamic interactions hyd , and external sources ext , . Similarly, the angular velocity changes according to the acting torque , , due to collisions and hydrodynamic interactions. These equations, together with the particle's position, are integrated in time via a Velocity Verlet scheme (Wachs 2019) with a constant time step size Δ = Δ /10. Consequently, ten particle simulation time steps are carried out within one fluid time step, which improves the overall accuracy of particle interactions and the efficiency of the simulation.
The collision forces and torques are determined via a discrete element method (DEM) that assumes a soft contact between overlapping rigid particles (cf. Cundall & Strack 1979). In our case, the normal and tangential collision components are given by a linear spring-dashpot model, similar to Costa et al. (2015) and Biegert et al. (2017). Following van der Hoef et al. (2006), the spring and damping coefficients of the normal collision model, and , are determined via the dry coefficient of restitution dry , a material parameter that is here chosen to be 0.97 (Vowinckel et al. 2021), and the collision time . The latter is chosen according to the findings in Rettinger & Rüde (2020) as = 4¯Δ /Δ , where¯is an average particle diameter, and ensures an adequate temporal resolution of the collision. As shown in Thornton et al. (2013), the spring and damping coefficient of the tangential model are related to the ones of the normal direction via the Poisson's ratio , such that = and = √ , with = 2(1 − )/(2 − ). The magnitude of the tangential collision force is limited by the Coulomb friction, determined as a product of the friction coefficient and the absolute value of the normal collision force. In the present simulations, we use = 0.22 and = 0.15 as reported in Joseph & Hunt (2004). The external force is given as the gravitational and buoyancy forces due to the gravitational acceleration , i.e. ext , = ( − ) , .

Fluid-particle coupling
To establish the coupling between the fluid and the granular phase in an accurate manner, we follow Rettinger & Rüde (2020) and distinguish between resolved and unresolved hydrodynamic forces to compute hyd , and hyd , . For the resolved part, we use the LBM-specific momentum exchange method as proposed by Aidun et al. (1998) to apply an explicit mapping of the particles onto the computational grid. This is achieved by flagging cells with their centers contained inside of particles as solid, effectively removing them from the fluid domain (cf. figure 1). This results in a sharp interface between the fluid and solid phase, along which no-slip boundary conditions for the fluid are applied. Here, we use the central linear interpolation (CLI) scheme of Ginzburg et al. (2008) that allows for second-order accurate results by including information about the exact surface position. The momentum exchanged locally with the particle due to its no-slip boundary condition is then integrated over the whole particle surface, as in Wen et al. (2014). Following Ladd (1994, this measure determines the resolved part of the fluid-particle interaction force fp , and torque fp , acting on this particle, which are averaged over two consecutive fluid time steps for improved stability. Solid cells that are no longer occupied by the particle due to its motion are converted back to fluid cells. Additionally, the otherwise missing PDF information is restored in these cells with an approach similar to Dorschner et al. (2015), using density and pressure tensor information from surrounding fluid cells and the particle's velocity.
As shown in Rettinger & Rüde (2020), this approach is able to reliably and accurately predict the resolved part of the fluid-particle interactions of single spheres. For two approaching particles, however, the mesh resolution of the narrow gap between the particles' surfaces is usually too coarse to fully resolve the strong lubrication interaction originating from the fluid that is being squeezed out of the gap of size . For those cases, a lubrication correction model must be applied that accounts for these unresolved forces (Nguyen & Ladd 2002;Biegert et al. 2017). Thus, the total hydrodynamic interaction force and torque on a particle is here computed as (2.4) These lubrication correction forces and torques explicitly account for the pair-wise lubrication forces and torques due to relative normal, tangential translational, and tangential rotational velocities, and are given in Rettinger & Rüde (2020). As suggested by validation studies therein, the normal and tangential lubrication corrections are only active for < 2Δ /3 and < Δ /2, respectively. As these corrections scale as lub,cor , ∝ −1 and lub,cor , ∝ ln( ), they would grow to infinity for vanishing gap sizes. Hence, a calibrated lower limit of lub ,min = (0.001 + 0.000035 , /Δ ) , /2 is applied in their calculation.

Simulation description
In this section, we detail the set up of the simulation, including the generation of the sediment beds, the physical parameterization, the description of the computational setup, and, finally, the evaluation of relevant rheological quantities.

Setup description
The general scenario is to consider linear shear flows with a constant shear rate = /ℎ across sediment beds of polydisperse, spherical particles (cf. figure 2), where is the velocity of the moving top wall, ℎ = − ℎ is the clear fluid height, is the vertical extent of the domain, and ℎ is the height of the sediment. To this end, we generate a grain size distribution with diameter values for particles by sampling from a log-normal distribution, defined by the parameters LN and 2 LN . Those parameters are related to the desired mean and variance 2 of the distribution via  Note, that we decided to use the arithmetic mean diameter for the parameterization instead of the median diameter ,50 as it is also well-defined for bidisperse grain size distributions. As will be detailed in §3.2, we target a numerical resolution of the mean diameter of¯/Δ = 20. Especially for large variances, care must be taken to maintain a reasonable numerical resolution for all particle sizes including its smallest values. Hence, we dismiss diameter values below 10 cells to guarantee a reasonable resolution of the flow field around the particles. The statistical properties of the polydisperse sediments including the ratio of largest to smallest diameter in the bed, given by ,max = max , and ,min = min , , can be found in table 1. Note that was chosen below 20 for strong polydispersity to compensate for the lower limit of admissible diameters and to obtain¯/Δ ≈ 20. We also use a log-normal distribution, albeit with a much smaller variance, for the monodisperse case as encountered in experimental studies (Boyer et al. 2011;Aussillous et al. 2013) to prevent an artificially close packing observable in perfectly mono-sized sphere beds.
Subsequently, the initial sediment beds for the main simulations of a fully coupled fluid-particle system are created by a precursor simulation without fluid. A constant density is assigned to the particles. Initially, they are placed inside a tall domain, with a uniform spacing in all directions that prevents potentially large overlaps, and given a random velocity. Due to gravity, they then settle due to gravity on a plate of size × = 51.2¯× 25.6¯= 1024 × 512 cells, where and are the streamwise and spanwise extent, respectively, of the horizontally periodic computational domain. The precursor simulations are run until all particles have come to rest to yield the initial bed height ℎ 0 for the main simulation. This state is typically achieved after some minutes of simulation time on a regular workstation. We noticed that this precursor simulation requires the same physical parameters, such as gravitational acceleration and submerged weight, as in the main simulation to prevent large accelerations followed by abrupt position changes in the initial phase of the main simulation. Since ℎ 0 can only be roughly estimated a priori, an iterative procedure is applied to find the right amount of particles necessary to achieve comparable bed heights among the different runs. In all cases, the bed is generated to obtain an initial bed height of approximately ℎ 0 = 340, i.e. ℎ 0 /¯= 17 (cf. table 1). This requires around 26000 particles for the monodisperse case to around 14500 particles for the strongly polydisperse setup. A visualization of the generated sediment beds and the diameter distribution for all four cases can be seen in figure 3.

Physical parameterization
The main simulation is executed in a cuboidal domain of size × × = 1024 × 512 × 480 cells. The domain is completely filled with a viscous fluid, defined by the kinematic viscosity and density . Periodic boundary conditions are applied in streamwise ( ) and spanwise ( ) direction, while no-slip boundaries are applied at the particle surface as well as the top and bottom planes bounding the vertical direction ( ). The top plane is moving in -direction with a constant velocity = 0.03 in lattice units. The sphere packing is initialized by the results from the precursor simulations to prescribe ℎ 0 . We fix all particles with a vertical center position smaller than 3/4¯throughout the simulation to form a bottom roughness. This measure prevents artificial slipping of the complete bed over the bottom plane (Jain et al. 2017;Biegert et al. 2017). A linear shear profile is assigned to the fluid above the sediment bed as an initial condition (cf. figure 3).
Apart from the density ratio / , we characterize the sediment mobility by the Shields parameter Θ: where = is the shear stress, and is the magnitude of the gravitational acceleration. Additionally, we define a particle Reynolds number =¯/ using = √︁ / . For those non-dimensional parameters, we choose Θ = 0.5, = 0.76, and / = 1.5 in all simulations to have comparable results. The value of the Shields parameter is well above the expected threshold for incipient motion, given as Θ ≈ 0.12 by Ouriemi et al. (2007), to ensure an adequate mobility of the particles. This results in a bulk Reynolds number based on channel properties = ℎ /(2 ) of around 14 and a Stokes number, =¯2 / , of around 0.85, which makes the simulations fall into the viscous regime (Bagnold 1954). Due to the low Reynolds number, we obtain a laminar Couette-like flow profile in the bulk region above the bed, where is constant. Finally, we define the reference time scale as ref =¯/ . We explicitly note that the set of physical parameters of the simulations is determined using the initial values of the bed and the fluid height, since ℎ becomes a result of the simulation and varies over time when the sediment bed dilates under shear as will be detailed in §3.3.
To ensure an accurate resolution  Specifically in the present work, each simulation run is executed for 48h on 7680 processes on the SuperMUC-NG supercomputer at LRZ in Garching, Germany. The resulting 2.5 × 10 8 grid cells, simulated for around 9 × 10 6 time steps in each case, make the studies at hand one of the largest and computationally most costly simulation campaigns of polydisperse sediment beds reported in literature.
Movies of the simulations are provided as supplementary material.

Evaluation procedure for simulation data
Since the goal of the present study is to investigate the rheological behavior of sediment beds in the framework of the ( )-rheology, we have to obtain the values for , = / , and = / . These quantities can be determined from vertical profiles of and (Houssais et al. 2016;Vowinckel et al. 2021). From the numerical simulations, we obtain high fidelity data of individual particle positions and velocities, as well as flow velocities as a function of time and space. To process the data for robust rheological interpretations, we apply spatial and temporal averaging.
As a first step, we perform spatial averaging and analyze it over time to determine the initialization period needed to obtain a statistically stationary state. This measure ensures that transient effects such as the dilation of the granular packing under shear and the initial sorting of the polydisperse grains are excluded from the statistical analysis (cf. Appendix A). We subdivide the domain into binned averaging volumes of size 0 = × × Δ , stacked vertically upon each other. In order to obtain the vertical particle volume fraction profile at a specific time , we make use of the particle diameter and its center coordinates. The horizontal planes between the stacked 0 slice each sphere into several sphere segments, whose volume can be determined analytically. We then add up all the volumes of the sphere segments within a 0 and divide this accumulated particle volume by the total averaging volume | 0 | to obtain the particle volume fraction ( , ), where the discrete vertical center coordinate of the respective 0 . As a next step, we apply a central moving average of width 10Δ , which corresponds to half of the mean particle diameter. This measure is needed to even out the layering at the sub-particle scale that introduces fluctuations within horizontally averaged profiles (Vowinckel et al. 2021). From these vertical profiles and with linear interpolation, we can evaluate the bed height ℎ ( ) given as the vertical position, for which (ℎ , ) = 0.1 (Kidanemariam & Uhlmann 2014). Note that other authors have used different threshold values for this definition (Houssais et al. 2016;Biegert et al. 2017), but due to the sharp gradient of the profile at the interface region, the actual value to determine ℎ does not have an impact on our analysis of rheological quantities. The temporal evolution of the bed height due to the movement of the top particle layer is illustrated in figure 4. It can be seen that when increasing the polydispersity of the bed, fluctuations in ℎ become larger and also, on average, the bed expands more.
Based on these evaluations, we define an instant of time that marks the beginning of our averaging time, 0 . As mentioned above, this is done to exclude the initial dilation phase of the sediment bed and, in particular, possible morphological effects due to vertical grain size segregation for the polydisperse cases, see Appendix A. Hence, no significant changes in the rheological quantities nor the local particle size distributions are observed during the evaluation period. The temporal averaging windows for the different cases are stated in table 1 and visualized in figure 4 as gray shaded areas. The slightly different end times originate from the different total run time of the simulations.
These considerations finally allow us to obtain the time-averaged particle volume fraction as where the angular brackets indicate averaging in time as implied by the subscript . Similarly, we evaluate the time-averaged bed height and state it in table 2. Analogously, we perform the spatial and temporal averaging of the streamwise fluid velocity . There, we define an indicator function Γ being 1 in the fluid and 0 otherwise that separates the fluid from the particle phase to compute so-called intrinsic spatial averages (Vowinckel et al. 2017(Vowinckel et al. , 2019b: where the subscript of the angular brackets now indicates spatial averaging. This is again followed by a central moving average. Temporal averaging as in Eq. (3.4) finally yields , , the vertical fluid profile consecutively averaged over space and time. We note that we observed temporal fluctuations in the instantaneous flow profiles within the bulk of the sediment bed, i.e. where the fluid and particle velocities are very small. Those fluctuations presumably originate from ongoing sorting effects inside the bed that appear over long time spans (Ferdowsi et al.

2017
). As such, longer simulation times would be desirable to increase the temporal averaging window and obtain a more robust statistical steady state. It was shown by Vowinckel et al. (2021), however, that unsteady effects are negligible when analyzing the rheological properties in the viscous regime.
We obtain the local shear rate as the spatial derivative of , . Owing to the spatial heterogeneity of our polydisperse sediment beds that may still be subject to ongoing sorting, we decided to use the absolute value of the local shear rate, i.e. | |, as a robust measure to compute the rheological quantities (Madraki et al. 2017). The actual shear stress is extracted from the bulk region of the flow, where it is constant due to the linear flow profile. The normalized shear stress values of all cases are reported in table 2, which are close to the target Shields number of 0.5. The granular pressure, on the other hand, is obtained from via which is the total submerged weight of the sediment bed. The complete data sets for all four simulation cases can be found in the supplementary data. Looking at the particle volume fraction profile, a layering is visible near the bottom plane (figure 5a), which is due to the ordered structure induced by the spheres mounted to the bottom plane. Therefore, we discard the data from the lower parts of the bed, i.e. where < 5¯, to exclude potential artefacts induced by the boundary condition of the bottom roughness.
We can directly obtain the maximum solid volume fraction from the particle volume fraction profile. To this end, we evaluate its average in the bulk region of the bed, i.e.
Its value for the different setups is given in table 2. As expected, increases with polydispersity since the voids between larger particles can be filled by smaller particles. The maximum packing fractions are close to the values commonly reported in literature for random close sphere packings with log-normal size distributions (Brouwers 2014;Farr 2013).
For brevity, we will omit the indication of the averaging operator and use instead of to denote the averaged particle volume fraction for the remainder of the work.  (2005) to show that the rheology of the fluid-particle mixture is governed by . Based on these considerations, Boyer et al. (2011) proposed the following empirical correlations as a rheological model, which became known as the ( )-rheology and reads in its most general form

Rheology of monodisperse sediment beds
The macroscopic friction coefficient, thus, has the two contributions and ℎ from frictionalcontact-based and hydrodynamic stresses, respectively. The expression of was originally proposed by Jop et al. (2005) and Cassar et al. (2005) while studying submarine granular flow down an inclined plane. Notably, the parameter represents the value of for which = ( 2 + 1 )/2, i.e. the average of 1 and 2 . This parameter can therefore be understood as the transition from a frictional dominated to a more suspended regime where binary particle collisions prevail and the role of hydrodynamic stress becomes increasingly important. The parameters 1 and are particle properties that represent the minimum friction and maximum particle volume fraction, respectively, for → 0, i.e. the jamming point of the dense suspension when the granular flow ceases. According to Cassar et al. (2005), 2 is the maximum value for the friction coefficient at higher shear rates, whereas, this value serves as the threshold that distinguishes the two contributions from particle contact and hydrodynamic interactions in the framework of Boyer et al. (2011). The coefficients = 1, = 5/2 can be determined from the analytical solution for effective viscosities of dilute suspensions originating from Einstein (1905), and is a parameter that has been determined empirically by best fit to experimental data (Morris & Boulay 1999;Boyer et al. 2011).
For the sake of the arguments that follow, we decided to deviate from the commonly encountered notation of , which has previously been denoted as 0 (e.g. In contrast, Tapia et al. (2019) investigated a region of ∈ [3 × 10 −4 , 10 −1 ] to address the effect of particle roughness on the rheology of dense suspensions. For that reason, they used slightly roughened (SR) and highly roughened (HR) spheres in their experiments. Instead of fitting the complete Boyer model (4.1), these authors suggested a simplified scaling, which only contains the √ term close to the jamming transition and used this approach to determine the friction factor at the jamming point by extrapolating their data. This approach worked very well for the given range of , but it also required a fitting of the coefficient that was, thus, found to be different from the Einstein formulation. Following the reasoning given in Tapia  a constant which implies 1 = 2 . This effectively removes the second term of from (4.1) and, thus, is not required for this analysis. A summary of the values that have been reported in literature and discussed in the preceding paragraphs is given in table 3. Note that the particles used in all of these experimental studies were monodisperse spheres.

Comparison to simulation results
In an effort to compare our simulation results against experimental data of pressure-imposed rheometry, we evaluate our data following the procedure described in §3.3 to extract all rheological quantities as vertical profiles through the sediment bed (cf. figure 5). Combining the data from these profiles, we are able to investigate and as a function of within the range ∈ [10 −9 , 10 3 ]. This analysis is shown in figure 6 for the monodisperse case. In the upper panel of this figure, the macroscopic friction factor is given as a function of the viscous number . For comparison, we plot our data together with the experimentally obtained data from Boyer et al. . Consequently, the simulation data is well predicted by the parameterized models (4.1) for > 10 −5 . This range is in agreement with the values used in these experimental studies to calibrate the coefficients 1 , 2 , and . For lower values of , our data underestimates the two correlations, which confirms the creep regime reported by Houssais et al. (2016) and visible in their data. In this regime, the plotted parameterizations of the model predict that levels off to a constant value, whereas the available data shows another significant shift towards a lower level of .
The simulation results for ( )/ match well with the experimental data of Boyer et al. .589. The latter shows some significant scatter, originating from the five distinct experiments varying the Shields numbers. Excellent agreement between our data and the rheology model is observed for the range ∈ [10 −9 , 1], which contains the range of viscous numbers used in Boyer et al.
(2011) to parameterize the model. For larger , the simulation data exhibits smaller values than either of the models. In this range, we observe a more rapid decrease of from to 0. This region corresponds to the interface between the densely packed sediment bed and free flow region. The deviations reflect the difficulty to use the empirical correlation of Boyer et al. (2011) in the extrapolated region of a more dilute regime (Vowinckel et al. 2021). By comparing the two parameterizations, we see that the parameter in (4.2) controls the viscous number range of this transition region. We note that the value of , used for the normalization of our simulation data, is 0.631 and thus larger than the ones from other studies. As already noted §3.3, our value of is close to the one reported for a random sphere packing which can be expected since it is obtained from the bulk region of the sediment bed, i.e., the region of vanishingly low shear rates and, consequently, small viscous numbers. This is in contrast to other studies (Boyer et al. 2011;Vowinckel et al. 2021), where stronger shearing was applied that led to a notable dilation of the suspension and, thus, a decrease in . Furthermore, Singh et al. (2018) observed a strong influence of the inter-particle friction coefficient on for sheared systems and found values of that are similar to ours for a friction coefficient of = 0.15. To focus on the general behavior of the ( ) relation rather than the limiting value, which is therefore different in our simulation but also in existing studies, we always present and analyze the normalized values in this work. This also effectively removes the dependence on from the the ( ) model (4.2). In summary, our data of the monodisperse case agrees well with existing experimental data and previously derived parameterizations of the rheology model. This overall confirms the validity of our simulation approach for densely packed sediment beds in shear flow and enables further predictive simulations. These studies will feature polydisperse setups for direct comparison with the monodisperse models. Furthermore, we observe a systematic shift in towards lower values for < 10 −5 , also present in the experimental data of Houssais et al. (2016). This range, however, was not addressed by Boyer et al. (2011) nor Houssais et al. (2016 and is thus not contained in the existing rheological model. In the following section, we will evaluate and enhance the parameterization of the empirical coefficients in (4.1) for the effects of polydispersity by focusing on the collisional and hydrodynamic regime for ∈ [10 −5 , 10 2 ]. We then proceed in §6 to study the creep regime in more detail and propose an extended model that is able to capture the observed behavior.

Simulation results
We now apply the same analysis as for the monodisperse case in §4.3 for the additional three setups of polydisperse sediment beds summarized in table 1 that reflect different degrees of  polydispersity as indicated by the variance of the grain size distribution. This analysis again yields and as a function of and is shown in figure 7. Similar to figure 6, the left column shows the macroscopic friction factor from our data together with the model parameterizations from Boyer et al. (2011) andHoussais et al. (2016). For increasing polydispersity, we observe a decrease of within the range ∈ [10 −5 , 10 0 ]. Note that and are plotted on logarithmic scales, i.e. even small deviations that become visible in this range are large in actual values, as can be seen in the respective insets. All cases reproduce the creep regime for < 10 −6 , as already observed for the monodisperse case. This effect becomes slightly more pronounced with increasing polydispersity.
The right column of figure 7 shows our data for the particle volume fraction over normalized by , and model parameterizations from Boyer et al. (2011) and Morris & Boulay (1999). There, the drop from to 0 occurs at lower values of when the polydispersity is increased, which results in a shift by up to one order of magnitude in for poly-100 compared to mono. An interesting feature emerges for values of around ≈ 10 −4 that can be seen most prominently for the poly-100 case where values larger than are observable. We found this to be a result of vertical sorting of the polydisperse sediment, where finer sediments from the topmost sediment layer translate to and accumulate in a lower layer, thereby increasing the particle volume fraction in this region.
Summarizing, increasing the polydispersity of the sediment bed while keeping all other physical parameters constant has a distinct effect on and as a function of . As a result, the agreement between the simulation data and the existing model parameterizations by Boyer et al. (2011), Houssais et al. (2016), and Morris & Boulay (1999 deteriorates with increasing polydispersity. In the following, we will enhance the parameterization of the rheological model in (4.1) and (4.2) for the effects of polydispersity by focusing on the frictional and hydrodynamic regime for ∈ [10 −5 , 10 2 ]. For now, we exclude the creep regime for the remainder of this section to provide a consistent comparison with the analyses of Boyer et al. (2011) andHoussais et al. (2016). However, we will study this regime in more detail in the subsequent section §6.

Effect of polydispersity on model parameterization
In order to improve the parameterization of equations (4.1) and (4.2), we evaluate the parameters 1 , 2 , , and determined from fits of our simulation results to reveal trends as a function of increasing polydispersity. To this end, we apply a fit of (4.1) and (4.2) to our data. We follow the reasoning of Boyer et al. (2011) and determine 1 , 2 , and as free parameters, while keeping = 5/2 and = 1 to recover the Einstein relation for the effective viscosity of dilute suspensions. Similar to Houssais et al. (2016), we apply the fit over the range ∈ [10 −5 , 10 2 ] and exclude the values for lower to focus on the regimes dominated by frictional and hydrodynamic stresses. Owing to the large value range over several orders of magnitude, we fit ln( ) to instead of directly. The resulting coefficients are reported in table 4, and the corresponding plots are additionally presented in figure 7. We explicitly note that is extracted from our simulation results as a quantity of the individual sediment bed and is not fitted here.
Comparing the case mono to Boyer et al. (2011), our values for 2 and are almost identical, and also agrees very well, but we found a value for 1 that is closer to the results of Houssais et al. (2016). This could be attributed to the material parameters that enter our particle contact algorithm described in §2.2, such as the restitution coefficient and friction coefficient, which are parameters that are not reported by neither one of these experimental studies.
For increasing polydispersity, the friction coefficients 1 and 2 decrease, while increases. Additionally, changes in the four cases as well, although the values remain on a very low level for all cases. A significant shift was detected from = 0.0042 to = 0.0006 for the cases poly-10 and poly-50, respectively, whereas remains on this lower level for poly-100, Owing to the large range of , it is challenging to precisely determine the exact value of via curve fitting.
In the case of , the fitted curves reproduce the position and extent of the drop from to 0 particle volume fraction very well. This is achieved by increasing for larger polydispersities, resulting in significantly larger values than given by Morris &Boulay (1999) andBoyer et al. (2011). Slight deviations of the simulation data from the fitted correlations can still be seen for ≈ 10 2 where the curves predict values larger than present in the data.
Generally, the fitted curves plotted in figure 7 show a very good agreement with the simulation results for the here considered range of viscous numbers. This confirms our assumption that an adequate parameterization of the existing models for ( ) and ( ) allows for an extension that takes polydispersity into account. In a next step, we attempt to formalize the observed trends in the obtained coefficients as functions of polydispersity.

Model parameterization as a function of polydispersity
From the fits to the four different simulated cases, we find that the coefficients entering (4.1) and (4.2) depend on the polydispersity of the sediment bed. The parameters 1 and 2 decrease when the polydispersity is increased, whereas increases. Even though seemingly decreases with increasing polydispersity, we refrain from interpreting these values as an actual trend due to the aforementioned difficulties in its determination. Based on these findings, we aim to extend the existing rheological model to incorporate polydispersity in a general way and without individual calibration or fitting. As such, it becomes readily applicable in macroscopic simulations and can significantly improve the predictions of the rheology of polydisperse sediment beds.
To this end, we have to select a parameter that characterizes polydispersity in a concise way. A set of possible parameters can be found in table 1, namely the variance of the underlying log-normal distribution as well as the the diameter ratio ,max / ,min . It is also reported in table 2 that these parameters directly influence the maximum particle volume fraction that indicates the jamming condition. Here, we choose to be the characteristic parameter as it is already present in the existing rheological framework as a key parameter. This choice of the governing parameter is in line with recent work by Pednekar et al. (2018) and the quantity can be obtained in a robust manner from either the vertical profile of the particle volume fraction or from ( ) as → 0. For an a priori determination of , a reasonable estimation can be obtained by assuming a perfect log-normal distribution and making use of available packing fraction predictors (e.g. Brouwers 2014; Farr 2013). Previous studies on dry granular flows have suggested to account for polydispersity by using the weighted arithmetic mean of the particle diameter in the definition of the inertial number (Tripathi & Khakhar 2011). Since this geometric quantity does not appear in the definitions of ( )-rheology framework, we identified as the more suitable measure to account for polydispersity of dense suspensions in a quantitative manner. In figure 8, the fitted coefficients are plotted as a function of . In a next step, a functional expression for each parameter is determined which describes the dependence on . For the three parameters with a clear trend, we assume a linear dependence on . This is the strongest assumption we can justify based on the number of data points available. For , we refrain from further assumptions and use the average of the fitted values, while also reporting its standard deviation. A sensitivity study revealed that the dependence on the exact value of is only weak, so that solely its order of magnitude, which is captured well by the average, has a significant effect. This justifies the model simplification and keeps the number of coefficients to a minimum. Applying a linear regression, the resulting correlations for each parameter are given as The above relations are plotted as orange solid lines in figure 8 as well, exhibiting a reasonable agreement to the values determined by individual fits. For a quantitative comparison, we assess the predictive power of the rheology model (4.1) and (4.2) to reproduce our observed simulation results using the parameters 1 , 2 , and proposed by Boyer et al. (2011), Houssais et al. (2016 and Morris & Boulay (1999), as well as the ones found by the individual fits performed in §5.2 and compare it against the prediction using the parameterization given by the calibrated expressions (5.1)-(5.4). To this end, we compute the 2 value as measure to quantify the agreement between observations and a prediction model where¯is the average value of all observations. The maximum 2 = 1, thus, indicates perfect agreement between the model prediction and the observations, whereas smaller values mean lower agreement.
The 2 values are reported in table 5, where again we use the logarithmized data to compute 2 for due to its large value range. Note that we evaluated the 2 for the range of ∈ [10 −5 , 10 2 ], which corresponds to the value range used for fitting and excludes the creep regime. For ( ), the parameterizations from Boyer et al. (2011) andHoussais et al. (2016) offer a fairly good predictive quality for the monodisperse case and then deviate for increasing polydispersity, which is in line with our previous observations. This is improved when applying the fitted coefficients which produces an almost perfect agreement in all four cases. Our expressions for 1 , 2 , and , (5.1)-(5.3), yield a performance very similar to the fitted parameters. In particular, this shows that the results are rather insensitive to the actual choice of as the values differ by one order of magnitude in the case of poly-100, which can be seen as an additional justification for assuming a constant . The same findings regarding the predictive quality can be reported for the particle volume fraction . The individual fits and the correlation for , (5.4), yield very good agreement for all the cases, whereas the parameterization by Boyer et al. (2011) and Morris & Boulay (1999), i.e.
From these results, we conclude that our approach of including the effect of polydispersity via a functional dependence of the coefficients on successfully improves the macroscopic rheology models. Since the maximum particle volume fraction already appears in the original model, this strategy can readily be integrated and applied in macroscopic modeling approaches.
For ( ), however, the region of small , and accordingly small , values can not be captured via the present formulation of (4.1). As such, the applicability would be limited to cases with > 10 −5 . To solve this issue, the model for ( ) has to be extended to explicitly account for the creep regime as will be detailed in the next section. The model for ( ), on the other hand, correctly predicts a constant value of for these small viscous numbers and is thus already applicable to this regime without further modifications.

Evaluation of the creep regime
The creep regime is characterized as a slow deformation of granular material under very low shear rates. In terms of the ( )-rheology, this becomes evident by a macroscopic friction factor that does not level off to a constant value in the frictional regime, but decreases to even smaller values for lower and lower viscous numbers. Assessing this regime is challenging, because it requires very low viscous numbers. In fact, to the knowledge of the authors, the only experimental campaign that was able to investigate the rheology of the creep regime for granular flows immersed in a viscous shearing fluid is the study of Houssais et al. (2016), who reported values down to = 10 −9 . However, their results are subject to a substantial amount of scatter in this range due to the general difficulty of measuring such small and in an experimental apparatus that cannot be fully shielded from external disturbances and may touch the sensor accuracy of the measurement instruments. Additionally, this study was carried out in an annular flume that introduces some artifacts due to the curved side walls. In our simulations of a straight horizontal domain with no side walls being present and the ability to control and evaluate the setup very accurately, these experimental imperfections are not an issue. Despite the differences in the experimental setup of Houssais et al. (2016) and our numerical simulation, we confirm the observation of the creep regime in our simulation data, as seen in figure 7, albeit with less scatter. This is true not only for the monodisperse case, that yields very good agreement with the experimental data of Houssais et al. (2016) across the entire range of (figures 11a and 6), but also for all other cases considered (figures 11b-d).
Houssais et al. (2016) perceived creep as localized, intermittent particle motion for which a description with temporally averaged quantities like and might be less appropriate. To gain more insight into the dynamics of the creep regime and its mechanisms, we turn to the instantaneous but still spatially-averaged profiles of . These are visualized over time in figure 9 for all four simulated cases. Note that the displayed vertical region is restricted to ∈ [5, 15]t o better focus on the creep regime. Furthermore, we plot the viscous number in terms of log 10 due to its large value range. In all cases, we observe a short start-up phase which is followed by a statistically stationary state with temporal as well as vertical fluctuations. These steady fluctuations agree qualitatively well with the ones reported for hard particles by Bouzid et al. (2015), who carried out two-dimensional simulations of sheared dry systems in the quasi-static limit. This observation is in line with the particle properties used in our study, where the restitution coefficient and the particle friction were chosen to reflect silica grains. Similar to the results by Bouzid et al. (2015), no burst-like behavior can be observed in figure 9. On the contrary, Bouzid et al. (2015) observed such intermittent motion only for soft particles with restitution coefficients as low as 0.1, which could then be better described by a non-local rheology (e.g. Kamrin & Koval 2012).
Recently, Gillissen & Ness (2020) showed that temporal fluctuations of , rather than its average, characterize the creep regime for inhomogeneous flow conditions. These fluctuations are seen as the reason why the ( )-rheology by Boyer et al. (2011), derived for homogeneous conditions, fails to capture the creep regime. Even though our considered setup is a homogeneous shear flow, we also observe significant fluctuations in this region of the bed. We, therefore, follow the same argument and evaluate the vertical root-mean-square profile rms . It is based on the deviations of the vertical instantaneous profiles from the temporally averaged one, evaluated over the same time span as the temporal average (excluding the initial start-up phase, cf. table 1).
This analysis of the vertical profiles of and rms is shown in figure 10 for the four simulated cases. We observe that for viscous numbers above 10 −6 (mono) to around 10 −5 (poly-100), the fluctuations are smaller than the average . This is in agreement with results reported by Gillissen & Ness (2020) for homogeneous shear, and thus similar flow conditions. Furthermore,  this range corresponds to the viscous numbers, for which the existing ( )-rheology was found to agree well with our simulation data, see §5. Turning towards the creep regime, corresponding to the lower layers of the bed, the fluctuations exceed the averaged value by around two orders of magnitude. This was not observed by Gillissen & Ness (2020) for the case of homogeneous shear flow, as they could not access such small viscous numbers, so that the focus of this study was on inhomogeneous, and rather distinct, flow conditions of a Kolmogorov flow. Interestingly, our evaluation also shows that the fluctuations surpass the average at larger viscous numbers of around 10 as well. This coincides with the bed load transport layer at the fluid-sediment interface and is the region where the particles move along the bed's surface in an intermittent fashion, as they temporarily get trapped between particles and then proceed to slide or roll over them. While the magnitude of these fluctuations thus might provide additional insight into the mechanisms of the creep regime, we note that the development of such rheological models is still an active field of research (Gillissen & Ness 2020). In particular, information about these fluctuations is usually not available in two-phase models and would require additional closure relations to be applicable there. Instead, we focus on the steady-state rheology and aim to include the creep regime as an extension to the existing ( )-rheology in the next sections.

Extension of model to creep
Since the data by Boyer et al. (2011) did not access such low viscous numbers, the description of this regime is, hence, lacking in the ( )-rheology. To this end, we follow the reasoning of Cassar et al. (2005) and Jop et al. (2005), and define a creep regime in addition to the frictional and hydrodynamic regime. Similarly to the frictional regime, this brings a lower and an upper  (4.1), we have shifted the lower limit of the macroscopic friction from 1 to 0 , whereas 1 becomes the upper limit of the creeping regime that centers around the viscous number of the creep regime, i.e., . The proposed extension (6.1) recovers the original formulation (4.1) by choosing = 0 or 0 = 1 . We explicitly note that we here aim to model the rheological behavior for very small, but non-zero viscous numbers, i.e., → 0. This quasi-static, but still dynamic, regime might thus be different from the static case at = 0 (Perrin et al. 2019).

Testing the extended model for the creep regime
Similar to §5.2, we apply curve fitting to find appropriate values for the newly introduced coefficients 0 ∈ [0, 1] and for all simulations conducted. To this end, we extend the range of to the full range observed in the simulations, i.e. ∈ [10 −9 , 10 2 ]. Since the extended formulation (6.1) is meant as an extension of the classical ( )-rheology (4.1), we keep the values of the previously determined coefficients 1 , 2 and as reported in table 4. This also effectively prevents possible overfitting.
The results are shown in figure 11, together with the existing parameterizations of the original model and the fits from §5.2. The obtained coefficients are given in the respective subcaption of the figures. In all cases, the fit of the extended model (green line) is able to follow the shift to  the creep regime and reproduces our simulation data very well, especially for the extended range ∈ [10 −9 , 10 −5 ]. We also note that the curves of the extended model and the fit from §5.2 (orange line) collapse for > 10 −4 , where the extension term for the creep regime effectively evaluates to 1 and thus reduces to the original model. Analyzing the trend of the values determined for the two new parameters 0 and , we again notice a decrease in the friction coefficient 0 with increasing polydispersity. This decrease, however, is less significant than before for 1 and 2 and a difference of only around 10% can be seen between the monodisperse case and the one with strongest polydispersity. Generally, 0 is about three times smaller than 1 . Determining the parameter faces similar challenges as discussed for before which thus shows no clear trend with polydispersity. It is obvious, however, that its value averages out around 10 −6 , which is more than three orders of magnitude smaller than and confirms the physical meaning of discussed above to describe the average value of for the creep regime.
Due to the observed marginal sensitivity of 0 and on the polydispersity, and the general difficulty of measurements for the creep regime, we do not attempt to express a functional dependence on as in the previous section. In order to obtain a general parameterization of the creep-extended model, we instead propose to use the following expressions, evaluated as the 2 ln( ( )) 2 ( )/  Table 6: 2 values of different parameterizations of the rheology model, (4.1) and (4.2), evaluated with respect to the simulated data for ∈ [10 −9 , 10 2 ]. The polydispersity extension developed in §5 features fits with coefficients from table 4, and the correlations (5.1)-(5.4). The creep extension, (6.1), in the current section re-uses these coefficients or correlations, respectively, and adds the coefficients from figure 11 for the individual fits or the correlations from (6.2)-(6.3).
average of the fitted coefficients: 0 = 0.082 ± 0.004, (6.2) = 1.30 × 10 −6 ± 2.91 × 10 −7 . (6.3) We evaluate the performance of our creep-extended rheology model by computing the 2 for the different empirical correlations over the entire range of ∈ [10 −9 , 10 3 ]. For that, we compare (i) (4.1) with the parameters of Boyer et al.  et al. (2016), but also to the previously developed polydisperse model from §5.3, the creep-extended rheology outperforms all other available correlations. The fact that we observe an almost perfect match for both, the fit and the correlations, confirms the validity of our approach to account for polydispersity.
For completeness, we also show the 2 values for ( ) over the extended range of , in contrast to the limited range used in table 5. From there, we see that the creep regime does not influence the predictive performance of the polydispersity-extended ( ) model from §5, since it is the region of constant particle volume fraction and thus already covered by the model (4.2). Overall, the parameterization of the creep-extended rheological model via the proposed correlations yields 2 values larger than 0.992 for all here considered cases for both, and , and without any further calibration. This is a significant improvement to the previous rheology model and its parameterizations.

Conclusion
In this work, we studied the rheological properties of polydisperse, densely packed sediment beds in a laminar shear flow through particle-resolved direct numerical simulations. This was achieved by large-scale 3D simulation domains using an efficiently coupled lattice Boltzmanndiscrete element method to fully resolve all relevant scales in space and time. In particular, particle collisions are modeled by linear spring-damper models in normal and tangential directions, with a Coulomb-like friction model. Additionally, a lubrication model is applied for short-range hydrodynamic interactions. Four different sediment beds were created in a precursor simulation ranging from monodisperse to strongly polydisperse with a maximum to minimum diameter ratio close to 10. As a key feature, the non-uniformity of the sediment yields increasing values for the maximum packing fraction. The beds consisted of up to 26000 particles, and the flow conditions were chosen to obtain several layers of mobile particles. As such, the present simulations are one of the most extensive numerical studies on mobile polydisperse sediment beds.
From the simulation results, we obtained depth-resolved spatially and temporally averaged profiles of rheological quantities. These enabled us to study the impact of polydispersity on the scaling of the macroscopic friction coefficient and the particle volume fraction as a function of the viscous number , i.e., the ( )-rheology. We compared our results to previous experimental studies of dense suspensions of neutrally buoyant spheres and sheared sediment beds of inertial particles and found excellent agreement for the monodisperse case. Owing to the wide value range of the viscous number, ∈ [10 −9 , 10 3 ], and the highly-resolved data, we were able to enhance the ( )-rheology and its parameterization for the effects of polydispersity and creeping flow. The effect of polydispersity has so far not been investigated for continuous grain-size distributions, and we addressed this issue by focusing on the frictional and hydrodynamic regimes. Based on our systematic simulation campaign, we derived an improved parameterization of the rheological model of Boyer et al. (2011) that explicitly accounts for polydispersity. This was achieved by expressing the two coefficients 1 and 2 , and the free parameter as functions of . The parameter is already present in the original rheological model and is here determined as the maximum observable packing fraction for a log-normal grain size distribution with a given variance, which determines the degree of polydispersity.
The effect of creep has so far been reported in Houssais et al. (2016) only, but this regime was excluded from the discussion of the rheology in this study. Our results confirm the existence of a creeping regime that is distinctively different from the well-known frictional and hydrodynamic regimes at higher viscous numbers (Boyer et al. 2011). For vanishing shear, the macroscopic friction levels off to a quasi-static, creeping state that yields values of , which are substantially lower than the frictional regime would suggest. This observation gave rise to the idea to enhance the ( )-rheology to explicitly account for the creep regime following the argument of Jop et al. (2005). This was done at the cost of introducing two additional parameters. However, we remark that these new parameters are physically based quantities related to particle properties as they express the quasi-static friction for the creeping state and the characteristic viscous number that describes the transition from the frictional to the creeping regime. These two parameters were determined by fitting the extended empirical correlation to our simulation data, and we found them to be less dependent on the maximum particle volume fraction. Compared to the frictional regime, the friction coefficient of the creeping regime is reduced by a factor of three.
Finally, our study demonstrates that particle properties that enter the ( )-rheology framework may change the entire system's rheological properties. Since the scaling laws obtained so far involve several idealizations and particular choices for the sediment material used, more work will be needed to explore the effects of different particle and flow properties on the rheological behavior of sheared sediment beds. This highlights another benefit of our simulation approach, where such changes can be made with ease, allowing for efficient parametric studies.
Supplementary data. Supplementary material and movies are available online.
Funding. B. V. gratefully acknowledges the support through the German Research Foundation (DFG) grant VO2413/2-1. U. R. gratefully acknowledges financial support by the Federal Ministry of Education and Research (BMBF) through the SKAMPY project, grant 01 ICH 15003 A, and by the Bavarian State Ministry of Science and the Arts through the Competence Network for Scientific High Performance Computing in Bavaria (KONWIHR).

Declaration of Interests.
The authors report no conflict of interest.

Appendix A. Vertical size segregation
For polydisperse sediment beds that are exposed to shear stress, it is known that a vertical size segregation sets in (e.g. Ferdowsi et al. 2017). Consequently, larger particles move to the top of the bed while smaller particles descend to lower sediment layers. A similar phenomenon, the brazil nut effect, can be observed in dry granular beds subjected to vibrations (Rosato et al. 1987).
We study the dynamics of this vertical sorting by assessing the composition of the topmost layers of the bed. To this end, we define that particles with a vertical center of mass position above ℎ top = 15.5¯belong to the bed's top region, which is roughly 2¯below the average sediment bed height ℎ , cf. table 2. We then sort these top topmost particles according to their diameters into bins of size¯/5. Evaluating the size distribution over time, we are able to investigate the size-based segregation in this top layer. This evaluation is shown in figure 12 for equally spaced time steps throughout the complete simulation, i.e., ∈ [0, 12000] ref . Since such an effect is not present in the monodisperse case, we exclude it from these discussions.
In all cases, we see a qualitatively similar behavior. The smaller size fractions, relative to the overall diameter distribution, decreases in number over time. These particles, thus, move to lower layers of the bed and the smallest particles almost vanish completely from the top layers. This process is initially very pronounced but then slows down gradually. At the same time, the number of larger particles increases in the upper layer, although the absolute change is significantly weaker than for the smaller ones. All these changes in the composition primarily happen during the initial stage of the simulation, so that a steady state develops after > 6000 ref . This indicates that the fast segregation process, as described by Ferdowsi et al. (2017), is already completed. Therefore, we do not expect further strong morphological changes during the second half of the simulation from which we obtain the data for our evaluations, cf. table 2. Since the present study focuses on sheared polydisperse sediments under well-developed conditions, this initial run-up phase was excluded from the statistical analysis presented in §4 - §6.