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

Abstract

Classical scaling relationships for rheological quantities such as the $\mu (J)$-rheology have become increasingly popular for closures of two-phase flow modelling. However, these frameworks have been derived for monodisperse particles. We aim to extend these considerations to sediment transport modelling by using a more realistic sediment composition. We investigate the rheological behaviour 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 10, 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.

1. 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 us to characterize sediment transport in laminar flows in terms of the rheology to investigate the fluid–particle mixture's deformation behaviour 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 with data previously obtained in rheometer studies with dense suspensions of neutrally buoyant particles (e.g. Morris & Boulay 1999; Boyer, Guazzelli & Pouliquen 2011). For these classical rheological investigations, a shear rate $\dot {\gamma }$ is applied to a dense granular material suspended in a fluid with viscosity $\eta _f$ to investigate the total shear stress $\tau$ acting on the fluid–particle mixture in the shearing direction and the imposed particle pressure $p_p$ in the wall-normal direction. The total shear comprises hydrodynamic and frictional interparticle stresses, with the latter becoming more important with increasing particle volume fraction $\phi$ (Gallier et al. 2014; Guazzelli & Pouliquen 2018; Vowinckel et al. 2021).

In this regard, two types of rheometer set-ups are possible. On the one hand, the volume-imposed 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, $\eta _s=\tau /\eta _f\dot {\gamma }$ and $\eta _n=p_p/\eta _f\dot {\gamma }$, respectively, and to derive empirical correlations for these two quantities as functions of $\phi$. 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 $St=\rho _p\dot {\gamma }d^2_p/\eta _f$ smaller than 10 (Bagnold 1954; Ness & Sun 2016), where $\rho _p$ is the particle density and $d_p$ is the characteristic particle diameter, this measure allowed Boyer et al. (2011) to define a macroscopic friction coefficient $\mu =\tau /p_p$ that depends on the viscous number $J=\eta _f\dot {\gamma }/p_p$. Based on this, the authors were able to propose empirical correlations for $\mu (J)$ and $\phi (J)$ that distinguish between stress contributions from particle contact and hydrodynamic interactions. This framework has become known as the $\mu (J)$-rheology. In this article, we will follow the nomenclature of Guazzelli & Pouliquen (2018) and use the symbol $J$ rather than $I_v$ 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 $p_p$ at some depth in the sediment bed is equal to the submerged weight of the overlying grains (Aussillous et al. 2013; Maurin, Chauchat & Frey 2016; Vowinckel et al. 2019a). This analogy is important for two-phase fluid sediment transport modelling (Jenkins & Hanes 1998; Hsu, Jenkins & Liu 2004), where the fluid–particle mixture is treated as two separated continua with interconnected conservation laws of mass and momentum (Ouriemi, Aussillous & Guazzelli 2009). The empirical correlations of the $\mu (J)$-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 $\mu (J)$ and $\phi (J)$ 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<\phi <0.58$ and $J>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 $J\approx 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 annular flume employed in this study. Hence, for cases where the modelled flow conditions exceed the range of the calibration data, the $\mu (J)$-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 $\mu (J)$-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 $\mu (J)$-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 $J$ and large $\phi$, the frictional interparticle 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 $\mu _1$ and maximum particle volume fraction $\phi _m$. For example, Boyer et al. (2011) reported $\mu _1=0.32$ and $\phi _m=0.585$ for the monodisperse case, but it has been shown by Tapia, Pouliquen & Guazzelli (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 $\mu (J)$-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, Vowinckel & Meiburg 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 mm in diameter. This calls for an extension of the $\mu (J)$-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 with the monodisperse case (Chang & Powell 1994; Gondret & Petit 1997). In these studies, the non-uniformity of the bidisperse grains was up to $d_{p,{max}}/d_{p,{min}}=13.75$, where $d_{p,{max}}$ and $d_{p,{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 $\phi _m=0.585$ for the monodisperse case (Boyer et al. 2011) to $\phi _m=0.64$. Consequently, models for $\phi _m$ in bidisperse volume-imposed rheometry were proposed by Dörr, Sadiki & Mehdizadeh (2013) and Mwasame, Wagner & Beris (2016) that can also be applied to polydisperse systems (Pednekar, Chun & Morris 2018).

As a next step, two-dimensional (2-D) discrete element method (DEM) simulations with grains of continuous polydispersity have been carried out where the fluid drag was approximated by Stokes drag and lubrication (Trulsson, Andreotti & Claudin 2012; Ness & Sun 2016) and the variation of the grain size was kept constant at $d_{p,{max}}/d_{p,{min}}=3.0$ and $1.4$, respectively. A recent study by Amarsid et al. (2017) extended these considerations to a lattice Boltzmann method (LBM)–DEM for simulations in two dimensions for $d_{p,max}/d_{p,min}=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, three-dimensional (3-D) simulations with a systematic focus on the degree of polydispersity in pressure-imposed rheometry or even sheared sediment beds have not been considered yet. The present study addresses this issue.

We employ the open-source simulation framework waLBerla (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. $St<10$. To this end, we utilize the combined LBM–DEM 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. 2021), 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 $J\approx 10^{-9}$. The investigated sediment beds have a non-uniformity of $d_{p,{max}}/d_{p,{min}}$ up to a factor of 10, 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 $\mu (J)$-rheology to account for the sediment polydispersity by linking the non-uniformity to the critical volume fraction $\phi _m$ and propose a straightforward extension to creeping flow conditions that recovers the original $\mu (J)$-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 set-up in § 3. We then infer the pressure-imposed rheology and validate our simulation approach in § 4 by comparing the monodisperse case with the experimental data of Boyer et al. (2011) and Houssais et al. (2016), including the classical empirical correlations of the $\mu (J)$-rheology (Boyer et al. 2011). Finally, we utilize the data from our simulation campaign to present extensions of the $\mu (J)$-rheology for polydispersity and creeping flow in § 5 and § 6, respectively.

2. Numerical method

For the numerical studies presented here, we couple the LBM for fluid flow with a DEM to account for particle interactions of polydisperse, spherical grains. This approach has proved 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 open-source high-performance framework waLBerla (cf. Bauer et al. 2020a), and its implementation can be found in the official software repository (https://walberla.net/). A sketch of the numerical scheme is presented in figure 1.

Figure 1. Schematic representation of the coupled LBM–DEM approach for fully resolved particulate flow simulations. The orange circles depict two colliding spheres, $i$ and $j$. The underlying uniform grid is used for the LBM, which simulates the fluid flow inside the fluid (light blue) cells. The solid (light brown) cells, whose centres are contained inside the particles, do not carry fluid information.

2.1. LBM

The LBM is a relatively recent approach for the simulation of viscous fluid flow. It describes the evolution of particle distribution functions (p.d.f.s) on a uniform computational grid and thereby fulfils 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 $D3Q19$ two-relaxation-time model of Ginzburg, Verhaeghe & d'Humieres (2008). The relaxation times, connected via the parameter $\varLambda =3/16$, determine the kinematic fluid viscosity $\nu _f$ and allow for accurate flow simulations. The local fluid pressure $p_f$ and velocity $\boldsymbol {u}_f$ are obtained via zeroth- and first-order moments of the p.d.f.s 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 $\Delta x = 1$, the time step size $\Delta t = 1$ and a reference fluid density of $\rho _f = 1$. These will be used in the remainder of this work.

2.2. DEM

The motion of a spherical particle $i$ can be described by the Newton–Euler equations

(2.1)$$\begin{gather} m_{p,i} \frac{\text{d} \boldsymbol{u}_{p,i}}{\text{d} t} = \boldsymbol{F}_{p,i} = \boldsymbol{F}_{p,i}^{col} + \boldsymbol{F}_{p,i}^{hyd} + \boldsymbol{F}_{p,i}^{ext}, \end{gather}$$

(2.2)$$\begin{gather}I_{p,i} \frac{\text{d} \boldsymbol{\omega}_{p,i}}{\text{d} t} = \boldsymbol{T}_{p,i} = \boldsymbol{T}_{p,i}^{col} + \boldsymbol{T}_{p,i}^{hyd}. \end{gather}$$

Here, $m_{p,i} = \rho _{p} V_{p,i}$ is the mass of the particle of density $\rho _{p}$ and volume $V_{p,i}$, and ${I_{p,i} = (m_{p,i} d_{p,i}^2)/10}$ is the moment of inertia for a sphere of diameter $d_{p,i}$. The temporal change of the particle's translational velocity is thus given by the acting forces $\boldsymbol {F}_{p,i}$, with contributions from the collisions $\boldsymbol {F}_{p,i}^{col}$, the hydrodynamic interactions $\boldsymbol {F}_{p,i}^{hyd}$ and external sources $\boldsymbol {F}_{p,i}^{ext}$. Similarly, the angular velocity changes according to the acting torque $\boldsymbol {T}_{p,i}$, 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 $\Delta t_p = \Delta t/10$. Consequently, 10 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 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, $k_n$ and $d_n$, are determined via the dry coefficient of restitution $e_{dry}$, a material parameter that is here chosen to be $0.97$ (Vowinckel et al. 2021), and the collision time $T_c$. The latter is chosen according to the findings in Rettinger & Rüde (2020) as $T_c = 4 \bar {d}_p \Delta t /\Delta x$, where $\bar {d}_p$ is an average particle diameter, and ensures an adequate temporal resolution of the collision. As shown in Thornton, Cummins & Cleary (2013), the spring and damping coefficient of the tangential model are related to the ones of the normal direction via Poisson's ratio $\nu _p$, such that $k_t = \kappa _p k_n$ and $d_t = \sqrt {\kappa _p}d_n$, with $\kappa _p = 2(1-\nu _p)/(2-\nu _p)$. The magnitude of the tangential collision force is limited by the Coulomb friction, determined as a product of the friction coefficient $\mu _p$ and the absolute value of the normal collision force. In the present simulations, we use $\nu _p = 0.22$ and $\mu _p=0.15$ as reported in Joseph & Hunt (2004).

The external force is given as the gravitational and buoyancy forces due to the gravitational acceleration $\boldsymbol {g}$, i.e. $\boldsymbol {F}_{p,i}^{ext} = (\rho _p-\rho _f) V_{p,i} \boldsymbol {g}$.

2.3. 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 $\boldsymbol {F}_{p,i}^{hyd}$ and $\boldsymbol {T}_{p,i}^{hyd}$. For the resolved part, we use the LBM-specific momentum exchange method as proposed by Aidun, Lu & Ding (1998) to apply an explicit mapping of the particles onto the computational grid. This is achieved by flagging cells with their centres 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 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 $\boldsymbol {F}_{p,i}^{fp}$ and torque $\boldsymbol {T}_{p,i}^{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 p.d.f. 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 $\delta _n$. 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 $i$ is here computed as

(2.3)$$\begin{gather} \boldsymbol{F}_{p,i}^{hyd} = \boldsymbol{F}_{p,i}^{fp} + \boldsymbol{F}_{p,i}^{lub,cor}, \end{gather}$$

(2.4)$$\begin{gather}\boldsymbol{T}_{p,i}^{hyd} = \boldsymbol{T}_{p,i}^{fp} + \boldsymbol{T}_{p,i}^{lub,cor}. \end{gather}$$

These lubrication correction forces and torques explicitly account for the pairwise 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 $\delta _n<2\Delta x/3$ and $\delta _n<\Delta x/2$, respectively. As these corrections scale as $\boldsymbol {F}_{p,i}^{lub,cor} \propto \delta _n^{-1}$ and $\boldsymbol {T}_{p,i}^{lub,cor} \propto \ln (\delta _n)$, they would grow to infinity for vanishing gap sizes. Hence, a calibrated lower limit of ${\delta _{n,{min}}^{lub}=(0.001 + 0.000035 d_{p,i}/\Delta x)\,d_{p,i}/2}$ is applied in their calculation.

3. 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 set-up and, finally, the evaluation of relevant rheological quantities.

3.1. Set-up description

The general scenario is to consider linear shear flows with a constant shear rate $\dot {\gamma }=U_w/h_f$ across sediment beds of polydisperse, spherical particles (cf. figure 2), where $U_w$ is the velocity of the moving top wall, $h_f=L_z-h_b$ is the clear fluid height, $L_z$ is the vertical extent of the domain and $h_b$ is the height of the sediment. To this end, we generate a grain size distribution with diameter values for $N_p$ particles by sampling from a log-normal distribution, defined by the parameters $\mu _{LN}$ and $\sigma _{LN}^2$. Those parameters are related to the desired mean $\mu _X$ and variance $\sigma _X^2$ of the distribution via

(3.1a,b)\begin{equation} \mu_{LN} = \ln\left(\frac{\mu_X^2}{\sqrt{\mu_X^2 + \sigma_X^2}}\right) \quad \text{and}\quad \sigma_{LN}^2 = \ln\left(1+ \frac{\sigma_X^2}{\mu_X^2}\right), \end{equation}

which yields the mean diameter

(3.2)\begin{equation} \bar{d}_p = \frac{1}{N_p} \sum_{i=1}^{N_p} d_{p,i}. \end{equation}

Note, that we decided to use the arithmetic mean diameter for the parameterization instead of the median diameter $d_{p,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 $\bar {d}_p/\Delta x=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.

Figure 2. Sketch of physical set-up as a side view, including a slice of the initial flow field above the sediment bed.

The statistical properties of the polydisperse sediments including the ratio of largest to smallest diameter in the bed, given by $d_{p,{max}} = \max _i d_{p,i}$ and $d_{p,{min}} = \min _i d_{p,i}$, can be found in table 1. Note that $\mu _X$ was chosen below $20$ for strong polydispersity to compensate for the lower limit of admissible diameters and to obtain $\bar {d}_p / \Delta x \approx 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.

Table 1. Parameters and properties of the different sediment beds, where length scales are expressed in lattice units.

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 $\rho _p$ 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 on a plate of size $L_x\times L_y=51.2\bar {d}_p \times 25.6 \bar {d}_p = 1024\times 512$ cells, where $L_x$ and $L_y$ 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 $h_b^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 $h_b^0$ can only be roughly estimated a priori, an iterative procedure is applied to find the right number of particles $N_p$ 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 $h_b^0 = 340$, i.e. $h_b^0 / \bar {d}_p = 17$ (cf. table 1). This requires around $26\,000$ particles for the monodisperse case to around $14\,500$ particles for the strongly polydisperse set-up. A visualization of the generated sediment beds and the diameter distribution for all four cases can be seen in figure 3.

Figure 3. The four different set-ups and their diameter distribution (from (a,b) to (g,h): mono; poly-10; poly-50; poly-100). Colouring of particles is according to the diameter with a logarithmic colour scale. See table 1 for detailed information about bed configurations. Along the diameter distribution, the cumulative distribution function (c.d.f.) based on a kernel density estimate is provided.

3.2. Physical parameterization

The main simulation is executed in a cuboidal domain of size $L_x\times L_y \times L_z = 1024 \times 512 \times 480$ cells. The domain is completely filled with a viscous fluid, defined by the kinematic viscosity $\nu _f$ and density $\rho _f$. Periodic boundary conditions are applied in the streamwise ($x$) and spanwise ($y$) directions, while no-slip boundaries are applied at the particle surface as well as the top and bottom planes bounding the vertical direction ($z$). The top plane is moving in the $x$-direction with a constant velocity $U_w = 0.03$ in lattice units. The sphere packing is initialized by the results from the precursor simulations to prescribe $h_{b}^0$. We fix all particles with a vertical centre position smaller than $3/4\bar {d}_p$ throughout the simulation to form a bottom roughness. This measure prevents artificial slipping of the complete bed over the bottom plane (Biegert et al. 2017; Jain, Vowinckel & Fröhlich 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 $\rho _p/\rho _f$, we characterize the sediment mobility by the Shields parameter $\varTheta$ as follows:

(3.3)\begin{equation} \varTheta = \frac{\tau}{g (\rho_p-\rho_f) \bar{d}_p}, \end{equation}

where $\tau = \rho _f \nu _f \dot {\gamma }$ is the shear stress and $g$ is the magnitude of the gravitational acceleration. Additionally, we define a particle Reynolds number $Re_p = u_\tau \bar {d}_p/\nu _f$ using $u_\tau = \sqrt {\tau /\rho _f}$.

For those non-dimensional parameters, we choose $\varTheta = 0.5$, $Re_p = 0.76$ and $\rho _p/\rho _f = 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 $\varTheta _c \approx 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, $Re_b = U_w h_f/(2 \nu _f)$, of around $14$ and a Stokes number, $St = \rho _p \bar {d}_p^2 \dot {\gamma }/\eta _f$, 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 $\tau$ is constant. Finally, we define the reference time scale as $t_{ref} = \bar {d}_p / U_w$. 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 $h_b$ 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 of fluid–particle interaction, a numerical resolution of approximately $20$ cells per mean diameter is chosen in all simulations, i.e. $\bar {d}_p/\Delta x \approx 20$ (Costa et al. 2015; Biegert et al. 2017; Rettinger & Rüde 2017, 2020). Since such a high resolution inherently renders the present numerical simulations computationally challenging, a performance-optimized implementation of the numerical methods as well as efficient communication routines must be applied to stay within adequate runtimes without exhausting computational resources (Eibl & Rüde 2018; Bauer, Köstler & Rüde 2020b). The details of our simulation approach are presented in Bauer et al. (2020a). The approach has successfully been applied in previous large-scale studies of particle-resolved simulations (e.g. Rettinger et al. 2017; Götz et al. 2010), where its excellent performance on high-performance computing clusters has been demonstrated. Specifically in the present work, each simulation run is executed for 48 hours on $7680$ processes on the SuperMUC-NG supercomputer at LRZ in Garching, Germany. The resulting $2.5\times 10^8$ grid cells, simulated for around $9\times 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 the literature.

Movies of the simulations are provided as supplementary material, available at https://doi.org/10.1017/jfm.2021.870, together with relevant simulation results.

3.3. Evaluation procedure for simulation data

Since the goal of the present study is to investigate the rheological behaviour of sediment beds in the framework of the $\mu (J)$-rheology, we have to obtain the values for $p_p$, $\mu =\tau /p_p$ and $J=\eta _f \dot {\gamma }/p_p$. These quantities can be determined from vertical profiles of $\dot {\gamma }$ and $\phi$ (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 analyse 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 $V_0=L_x \times L_y \times \Delta x$, stacked vertically upon each other. In order to obtain the vertical particle volume fraction profile at a specific time $t$, we make use of the particle diameter and its centre coordinates. The horizontal planes between the stacked $V_0$ slice each sphere into several sphere segments, whose volume $V_s$ can be determined analytically. We then add up all the volumes of the sphere segments within a $V_0$ and divide this accumulated particle volume $\sum V_s$ by the total averaging volume $|V_0|$ to obtain the particle volume fraction $\phi (z,t)$, where $z$ the discrete vertical centre coordinate of the respective $V_0$. As a next step, we apply a central moving average of width $10\Delta x$, which corresponds to half of the mean particle diameter. This measure is needed to even out the layering at the subparticle 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 $h_b(t)$ given as the vertical position, for which $\phi (h_b,t) = 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 $h_b$ 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 $h_b$ become larger and also, on average, the bed expands more.

Figure 4. Bed height $h_b$ as a function of time extracted from the instantaneous vertical volume fraction profiles for all simulation set-ups. The grey area depicts the region used for temporal averaging. For cases (a) mono, (b) poly-10, (c) poly-50 and (d) poly-100.

Based on these evaluations, we define an instant of time that marks the beginning of our averaging time, $t_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 grey 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

(3.4)\begin{equation} \langle \phi \rangle_{t}(z) = \frac{1}{t_1-t_0} \int_{t_0}^{t_1} \phi(z,t) \,\text{d}t, \end{equation}

where the angular brackets indicate averaging in time as implied by the subscript $t$. Similarly, we evaluate the time-averaged bed height and state it in table 2.

Table 2. Sediment bed and flow quantities extracted from the simulation data, together with duration of the time-averaging period.

Analogously, we perform the spatial and temporal averaging of the streamwise fluid velocity $u_f$. We define an indicator function $\varGamma$ 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, 2019b)

(3.5)\begin{equation} \langle u_f \rangle_{V}(z,t) = \frac{1}{\int_{V_0} \varGamma \,\text{d}V} \int_{V_0} \varGamma u_f(x,y,z,t) \,\text{d}V, \end{equation}

where the subscript $V$ of the angular brackets now indicates spatial averaging. This is again followed by a central moving average. Temporal averaging as in (3.4) finally yields $\langle u_f \rangle _{V,t}$, 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 analysing the rheological properties in the viscous regime.

We obtain the local shear rate as the spatial derivative of $\langle u_f \rangle _{V,t}$. 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. $\lvert \dot {\gamma }\rvert$, as a robust measure to compute the rheological quantities (Madraki et al. 2017). The actual shear stress $\tau$ 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 $\langle \phi \rangle _t$ via

(3.6)\begin{equation} p_p(z) = \left(\rho_p-\rho_f\right) g \int_{z}^{\infty} \langle\phi\rangle_t(z') \,\text{d}z'. \end{equation}

This definition is in line with the one proposed by the two-phase model of Aussillous et al. (2013) and successfully used in the analysis of Vowinckel et al. (2021). Note that we do not introduce an artificial confining pressure $P_0$ at the top wall as suggested by Houssais et al. (2016), because our simulation data yield full information of vertically resolved porosity profiles across the entire depth of the channel. These data allow for a straightforward computation of the vertical profiles of $\mu$ and $J$. The final profiles of the relevant quantities are exemplified in figure 5 by showing the results for the monodisperse case. In this figure, the granular pressure is normalized by $P_{tot} = (\rho _p-\rho _f) g \int _{0}^{\infty } \langle \phi \rangle _t(z') \,\text {d}z'$, 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 $z < 5\bar {d}_p$, to exclude potential artefacts induced by the boundary condition of the bottom roughness.

Figure 5. Spatially and temporally averaged profiles of different quantities for the monodisperse case. The dashed horizontal line represents the bed height. The solid horizontal line is at $z=5\bar {d}_p$ and all profile data below is discarded in further analysis. The data of the vertical profiles for all simulation runs are provided as supplemental material.

We can directly obtain the maximum solid volume fraction $\phi _m$ from the particle volume fraction profile. To this end, we evaluate its average in the bulk region of the bed, i.e.

(3.7)\begin{equation} \phi_m = \frac{1}{5\bar{d}_p} \int_{5 \bar{d}_p}^{10 \bar{d}_p} \langle\phi\rangle_t(z) \,\text{d}z. \end{equation}

Its value for the different set-ups is given in table 2. As expected, $\phi _m$ 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 the literature for random close sphere packings with log-normal size distributions (Farr 2013; Brouwers 2014).

For brevity, we will omit the indication of the averaging operator and use $\phi$ instead of $\langle \phi \rangle _t$ to denote the averaged particle volume fraction for the remainder of the work.

4. Rheology of monodisperse sediment beds

4.1. Rheological model for dense suspensions

The rheology of monodisperse, neutrally buoyant, spherical particles in a viscous fluid has been assessed experimentally by shearing walls that impose a constant volume on the fluid–particle mixture (e.g Krieger & Dougherty 1959; Morris & Boulay 1999; Stickel & Powell 2005; Guazzelli & Morris 2011). This approach is commonly referred to as volume-imposed rheometry. The scenario has been extended to a pressure-imposed rheometry, where a constant confining pressure is applied on the top wall that remains movable in the vertical direction. This measure allows us to investigate the dilation/consolidation of a granular suspension under varying shear (e.g. Boyer et al. 2011; Dagois-Bohy et al. 2015; Tapia et al. 2019). As already laid out in the introduction, this scenario bears a straightforward analogy to the shearing of sediment beds. Hence, the pressure-imposed rheometry and the corresponding empirical correlations derived from the rheological experiments to predict the macroscopic friction and the particle volume fraction as functions of the viscous number $J=\eta _f \dot {\gamma }/p_p$ are the focus of this work.

Using their experimental apparatus, Boyer et al. (2011) followed the argument of Cassar, Nicolas & Pouliquen (2005) to show that the rheology of the fluid–particle mixture is governed by $J$. Based on these considerations, Boyer et al. (2011) proposed the following empirical correlations as a rheological model, which became known as the $\mu (J)$-rheology and reads in its most general form:

(4.1)$$\begin{gather} \mu(J) = \underbrace{\mu_1 + \frac{\mu_2-\mu_1}{1+J_f/J}}_{\mu^f(J)}+ \underbrace{a_\mu J^{1/2} + b_\mu J}_{\mu^h(J)}, \end{gather}$$

(4.2)$$\begin{gather}\phi(J) = \frac{\phi_m}{1+(K_n J)^{1/2}}. \end{gather}$$

The macroscopic friction coefficient, thus, has the two contributions, $\mu ^f$ and $\mu ^h$, from frictional-contact-based and hydrodynamic stresses, respectively. The expression of $\mu ^f$ was originally proposed by Jop, Forterre & Pouliquen (2005) and Cassar et al. (2005) while studying submarine granular flow down an inclined plane. Notably, the parameter $J_f$ represents the value of $J$ for which $\mu ^f = (\mu _2 + \mu _1)/2$, i.e. the average of $\mu _1$ and $\mu _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 $\mu _1$ and $\phi _m$ are particle properties that represent the minimum friction and maximum particle volume fraction, respectively, for $J\rightarrow 0$, i.e. the jamming point of the dense suspension when the granular flow ceases. According to Cassar et al. (2005), $\mu _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 $a_\mu =1$, $b_\mu =5/2 \phi _m$ can be determined from the analytical solution for effective viscosities of dilute suspensions originating from Einstein (1905), and $K_n$ 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 $J_f$, which has previously been denoted as $I_0$ (e.g. Cassar et al. 2005; Boyer et al. 2011; Houssais et al. 2016) or $J_0$ (e.g. Guazzelli & Pouliquen 2018; Vowinckel et al. 2021).

4.2. Existing model parameterizations

In the work of Boyer et al. (2011), viscous numbers in the range $J \in [10^{-6},10^{-1}]$ were investigated. Since $\lim _{J\rightarrow 0}\mu ^f(J) = \mu _1$ and $\lim _{J\rightarrow 0}\phi (J) = \phi _m$, the parameters $\mu _1=0.32$ and $\phi _m = 0.585$ were obtained within the lower limit of $J$. Additionally, the parameters $\mu _2=0.7$ and $J_f=0.005$ were determined by fitting to the experimental data. The coefficient $K_n$ was evaluated as $K_n=1$ by Boyer et al. (2011), whereas Morris & Boulay (1999) found a value of $K_n=0.75$ in their experiments on shear-induced particle migration.

Recently, further experimental studies of an annular flume set-up with monodisperse spheres were reported by Houssais et al. (2016) and Tapia et al. (2019), which differ most notably in the range of measured $J$ values. In Houssais et al. (2016), a sediment bed of monodisperse spheres was sheared by a laminar Couette flow to obtain values of $J \in [10^{-9},10]$, which extended the data range to significantly lower $J$. This study revealed a novel regime for $\mu$, labelled as the creep regime and it is discussed in more detail in § 6. To provide a comparison with (4.2), Houssais et al. (2016) decided to exclude these low $J$-values from their analysis to obtain fitted coefficients for the region $J \in [3\times 10^{-5},2]$ that show very good agreement with the results of Boyer et al. (2011).

In contrast, Tapia et al. (2019) investigated a region of $J \in [3\times 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 $\sqrt J$ 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 $J$, but it also required a fitting of the coefficient $a_\mu$ that was, thus, found to be different from the Einstein formulation. Following the reasoning given in Tapia et al. (2019), they assumed a constant $\mu ^f$ which implies $\mu _1 = \mu _2$. This effectively removes the second term of $\mu ^f$ from (4.1) and, thus, $J_f$ is not required for this analysis.

A summary of the values that have been reported in the 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.

Table 3. Summary of previous work in the context of the $\mu (J)$-rheology, (4.1), together with the reported coefficients. The symbol (*) indicates the values for $J$ and $J_f$ were adapted to match our definition of the viscous number and the symbol (${\dagger}$) indicates the range used for fitting.

4.3. Comparison with 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 $\mu$ and $\phi$ as functions of $J$ within the range $J \in [10^{-9},10^3]$. This analysis is shown in figure 6 for the monodisperse case. In panel (a) of this figure, the macroscopic friction factor $\mu$ is given as a function of the viscous number $J$. For comparison, we plot our data together with the experimentally obtained data from Boyer et al. (2011) and Houssais et al. (2016), as well as the proposed parameterization therein of the $\mu (J)$ model (4.1) as summarized in table 3. Panel (b) of the same figure shows our data for the particle volume fraction $\phi$ over $J$ normalized by $\phi _m$, and the predictions using (4.2) with the coefficient $K_n$ from Boyer et al. (2011) and from Morris & Boulay (1999).

Figure 6. Rheological quantities ((a), $\mu$; (b), $\phi$) as a function of viscous number $J$ for a monodisperse sediment bed. Data from the present monodisperse simulation is shown, together with experimental data from Boyer et al. (2011) and Houssais et al. (2016). Additionally, curves of (4.1) and (4.2) are shown, parameterized as proposed by Boyer et al. (2011), Houssais et al. (2016) and Morris & Boulay (1999) (cf. table 2).

Comparing our simulation results of $\mu (J)$ with the existing experimental data shows a very good agreement, in particular with the data from Houssais et al. (2016) over the complete range of $J$. Consequently, the simulation data is well predicted by the parameterized models (4.1) for $J > 10^{-5}$. This range is in agreement with the values used in these experimental studies to calibrate the coefficients $\mu _1$, $\mu _2$ and $J_f$. For lower values of $J$, 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 $\mu$ levels off to a constant value, whereas the available data shows another significant shift towards a lower level of $\mu$.

The simulation results for $\phi (J)/\phi _m$ match well with the experimental data of Boyer et al. (2011), normalized by $\phi _m=0.585$, and Houssais et al. (2016), normalized by $\phi _m=0.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 $J \in [10^{-9},1]$, which contains the range of viscous numbers used in Boyer et al. (2011) to parameterize the model. For larger $J$, the simulation data exhibits smaller $\phi$ values than either of the models. In this range, we observe a more rapid decrease of $\phi$ from $\phi _m$ to 0. This region corresponds to the interface between the densely packed sediment bed and free flow region. The deviations reflect the difficulty in using 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 $K_n$ in (4.2) controls the viscous number range of this transition region. We note that the value of $\phi _m$, used for the normalization of our simulation data, is 0.631 and thus larger than the ones from other studies. As already noted in § 3.3, our value of $\phi _m$ 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 $\phi _m$. Furthermore, Singh et al. (2018) observed a strong influence of the interparticle friction coefficient $\mu _p$ on $\phi _m$ for sheared systems and found values of $\phi _m$ that are similar to ours for a friction coefficient of $\mu _p=0.15$. To focus on the general behaviour of the $\phi (J)$ relation rather than the limiting value, which is therefore different in our simulation but also in existing studies, we always present and analyse the normalized $\phi$ values in this work. This also effectively removes the dependence on $\phi _m$ from the $\phi (J)$ 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 set-ups for direct comparison with the monodisperse models. Furthermore, we observe a systematic shift in $\mu$ towards lower values for $J<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 $J \in [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 behaviour.

5. Rheological model for polydisperse sediment beds

5.1. Simulation results

We now apply the same analysis as for the monodisperse case in § 4.3 for the additional three set-ups 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 $\mu$ and $\phi$ as functions of $J$ and is shown in figure 7.

Figure 7. Rheological quantities ((a,c,e,g), $\mu$; (b,d,f,h), $\phi$) as a function of viscous number $J$ for the four different set-ups (from (a,b) to (g,h): mono, poly-10, poly-50, poly-100). Colour and style as in figure 6. Additionally, best fits as explained in § 5.2 are given as orange curves. The insets in (a,c,e,g) magnify the region of small viscous numbers, using a linear axis for $\mu$.

Similar to figure 6, figure 7(a,c,e,g) shows the macroscopic friction factor $\mu$ from our data together with the model parameterizations from Boyer et al. (2011) and Houssais et al. (2016). For increasing polydispersity, we observe a decrease of $\mu$ within the range $J \in [10^{-5},10^0]$. Note that $\mu$ and $J$ 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 $J < 10^{-6}$, as already observed for the monodisperse case. This effect becomes slightly more pronounced with increasing polydispersity.

Figure 7(b,d,f,h) shows our data for the particle volume fraction $\phi$ over $J$ normalized by $\phi _m$, and model parameterizations from Boyer et al. (2011) and Morris & Boulay (1999). There, the drop from $\phi _m$ to 0 occurs at lower values of $J$ when the polydispersity is increased, which results in a shift by up to one order of magnitude in $J$ for poly-100 compared with mono. An interesting feature emerges for values of $\phi$ around $J\approx 10^{-4}$ that can be seen most prominently for the poly-100 case where values larger than $\phi _m$ 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 $\mu$ and $\phi$ as functions of $J$. 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 $J \in [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) and Houssais et al. (2016). However, we will study this regime in more detail in the subsequent section, § 6.

5.2. Effect of polydispersity on model parameterization

In order to improve the parameterization of (4.1) and (4.2), we evaluate the parameters $\mu _1, \mu _2, J_f$ and $K_n$ 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 $\mu _1$, $\mu _2$ and $J_f$ as free parameters, while keeping $a_\mu =5/2\phi _m$ and $b_\mu =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 $J \in [10^{-5},10^2]$ and exclude the values for lower $J$ 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 (\mu )$ to $J$ instead of $\mu$ directly. The resulting coefficients are reported in table 4, and the corresponding plots are additionally presented in figure 7. We explicitly note that $\phi _m$ is extracted from our simulation results as a quantity of the individual sediment bed and is not fitted here.

Table 4. Coefficients applied for the equations of the $\mu (J)$- and $\phi (J)$-rheology for the curves shown in figure 7, with $\phi _m$ from table 2. The fits are obtained using data of $J\in [10^{-5},10^2]$.

Comparing the case mono with Boyer et al. (2011), our values for $\mu _2$ and $J_f$ are almost identical, and $K_n$ also agrees very well, but we found a value for $\mu _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 either one of these experimental studies.

For increasing polydispersity, the friction coefficients $\mu _1$ and $\mu _2$ decrease, while $K_n$ increases. Additionally, $J_f$ 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 $J_f=0.0042$ to $J_f=0.0006$ for the cases poly-10 and poly-50, respectively, whereas $J_f$ remains on this lower level for poly-100. Owing to the large range of $J$, it is challenging to precisely determine the exact value of $J_f$ via curve fitting.

In the case of $\phi$, the fitted curves reproduce the position and extent of the drop from $\phi _m$ to 0 particle volume fraction very well. This is achieved by increasing $K_n$ for larger polydispersities, resulting in significantly larger values than given by Morris & Boulay (1999) and Boyer et al. (2011). Slight deviations of the simulation data from the fitted correlations can still be seen for $J \approx 10^2$ where the curves predict values larger than those present in the data.

Generally, the fitted curves plotted in figure 7 show a very good agreement with the simulation results for the range of viscous numbers considered here. This confirms our assumption that an adequate parameterization of the existing models for $\mu (J)$ and $\phi (J)$ 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.

5.3. 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 $\mu _1$ and $\mu _2$ decrease when the polydispersity is increased, whereas $K_n$ increases. Even though $J_f$ 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 diameter ratio $d_{p,{max}}/d_{p,{min}}$. It is also reported in table 2 that these parameters directly influence the maximum particle volume fraction $\phi _m$ that indicates the jamming condition. Here, we choose $\phi _m$ 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 $\phi (J)$ as $J\rightarrow 0$. For an a priori determination of $\phi _m$, a reasonable estimation can be obtained by assuming a perfect log-normal distribution and making use of available packing fraction predictors (e.g. Farr 2013; Brouwers 2014). Alternatively, for our specific case of a truncated log-normal distribution, the relation

(5.1)\begin{equation} \phi_m = 0.000633\, \sigma_X^2 + 0.635 \end{equation}

can be used to approximate $\phi _m$ as a function of the variance $\sigma _X^2$. Note, however, that this empirical correlation is specific to the grain size distributions used in the present study and different types of sediments may require adjustments of this empirical correlation. Previous studies on dry granular flows have suggested accounting 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 the $\mu (J)$-rheology framework, we identified $\phi _m$ 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 $\phi _m$.

Figure 8. Fitted coefficients (blue) from table 4 as a function of parameter $\phi _m$, which is used to describe polydispersity. Additionally, the correlations (5.2)–(5.5) are included as orange lines.

In a next step, a functional expression for each parameter is determined which describes the dependence on $\phi _m$. For the three parameters with a clear trend, we assume a linear dependence on $\phi _m$. This is the strongest assumption we can justify based on the number of data points available. For $J_f$, 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 $J_f$ 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

(5.2)$$\begin{gather} \mu_1 ={-}0.980 \, \phi_m +0.872, \end{gather}$$

(5.3)$$\begin{gather}\mu_2 ={-}6.157 \, \phi_m +4.554, \end{gather}$$

(5.4)$$\begin{gather}J_f = 0.0027 \pm 0.0024, \end{gather}$$

(5.5)$$\begin{gather}K_n = 60.490 \, \phi_m -37.008. \end{gather}$$

The above relations are plotted as orange solid lines in figure 8 as well, exhibiting a reasonable agreement with 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 $\mu _1$, $\mu _2$, $J_f$ and $K_n$ 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 them against the prediction using the parameterization given by the calibrated expressions (5.2)–(5.5). To this end, we compute the $R^2$ value as a measure to quantify the agreement between observations $o$ and a prediction model $m$ as

(5.6)\begin{equation} R^2 = 1- \frac{\sum_i (o_i - m_i)^2}{\sum_i (o_i - \bar{o})^2}, \end{equation}

where $\bar {o}$ is the average value of all observations. The maximum $R^2=1$, thus, indicates perfect agreement between the model prediction and the observations, whereas smaller values mean lower agreement.

The $R^2$ values are reported in table 5, where again we use the logarithmized data to compute $R^2$ for $\mu$ due to its large value range. Note that we evaluated the $R^2$ for the range of $J \in [10^{-5},10^2]$, which corresponds to the value range used for fitting and excludes the creep regime. For $\mu (J)$, the parameterizations from Boyer et al. (2011) and Houssais 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 $\mu _1$, $\mu _2$ and $J_f$, (5.2)–(5.4), yield a performance very similar to the fitted parameters. In particular, this shows that the results are rather insensitive to the actual choice of $J_f$ 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 $J_f$. The same findings regarding the predictive quality can be reported for the particle volume fraction $\phi$. The individual fits and the correlation for $K_n$, (5.5), yield very good agreement for all the cases, whereas the parameterization by Boyer et al. (2011) and Morris & Boulay (1999), i.e. $K_n= 1$ and $K_n= 0.75$, respectively, are not as accurate.

Table 5. The $R^2$ values for different parameterizations of the rheology model, (4.1) and (4.2), evaluated with respect to the simulated data for $J \in [10^{-5},10^2]$, thus excluding the creep regime. Present contributions consist of individual fits for each case with coefficients from table 4, and the novel correlations (5.2)–(5.5) taking into account polydispersity. Bold indicates the most relevant data.

From these results, we conclude that our approach of including the effect of polydispersity via a functional dependence of the coefficients on $\phi _m$ 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 modelling approaches. Owing to its simplicity, the proposed correlations should be used with care for different bed compositions as they might yield $\phi _m$-values outside the studied range that may lead to unrealistic values for $\mu _1$, $\mu _2$ or $K_n$.

For $\mu (J)$, however, the region of small $J$, and accordingly small $\mu$, values cannot be captured via the present formulation of (4.1). As such, the applicability would be limited to cases with $J>10^{-5}$. To solve this issue, the model for $\mu (J)$ has to be extended to explicitly account for the creep regime as will be detailed in the next section. The model for $\phi (J)$, on the other hand, correctly predicts a constant value of $\phi _m$ for these small viscous numbers and is thus already applicable to this regime without further modifications.

6. Rheological model for creep regime

6.1. 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 $\mu (J)$-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 $J = 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 $J$ and $\mu$ 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 artefacts due to the curved sidewalls. In our simulations of a straight horizontal domain with no sidewalls being present and the ability to control and evaluate the set-up very accurately, these experimental imperfections are not an issue. Despite the differences in the experimental set-up 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 $J$ (figures 6 and 7a), but also for all other cases considered (figure 7c,e,g).

Houssais et al. (2016) perceived creep as localized, intermittent particle motion for which a description with temporally averaged quantities like $J$ and $\mu$ 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 $J$. These are visualized over time in figure 9 for all four simulated cases. Note that the displayed vertical region is restricted to $z\in [5,15]\bar {d}_p$ to better focus on the creep regime. Furthermore, we plot the viscous number in terms of $\log _{10}J$ 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 2-D 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 behaviour 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).

Figure 9. Temporal evolution of vertical $J$-profiles. Due to the range of values, we plot $\log _{10}(J)$ to indicate the order of magnitude and choose the colour scale to focus on very low viscous numbers. For cases (a) mono, (b) poly-10, (c) poly-50, (d) poly-100.

Recently, Gillissen & Ness (2020) showed that temporal fluctuations of $J$, rather than its average, characterize the creep regime for inhomogeneous flow conditions. These fluctuations are seen as the reason why the $\mu (J)$-rheology by Boyer et al. (2011), derived for homogeneous conditions, fails to capture the creep regime. Even though our considered set-up 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 (r.m.s.) profile $J_{rms}$. It is based on the deviations of the vertical instantaneous $J$ 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 $J$ and $J_{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 $J$. 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 $\mu (J)$-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.

Figure 10. Vertical profiles of the time-averaged $J$ and the r.m.s. value of its fluctuations, evaluated over the same time span $[t_0,t_1]$ as the temporal averaging given in table 1. For cases (a) mono, (b) poly-10, (c) poly-50 and (d) poly-100.

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 $\mu (J)$-rheology in the next sections.

6.2. 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 $\mu (J)$-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 limit of macroscopic friction, so that there remains a smooth transition in-between the different regimes. This yields the following extension of (4.1) to adequately capture the creep regime in the rheological framework:

(6.1)\begin{equation} \mu(J) = \underbrace{\mu_0 + \frac{\mu_1-\mu_0}{1+J_c/J}}_{\mu^c} + \underbrace{\frac{\mu_2-\mu_1}{1+J_f/J}}_{\mu^f}+ \underbrace{\frac{5}{2} \phi_m J^{1/2} + J}_{\mu^h}. \end{equation}

In comparison with the original model of Boyer et al. (2011), (4.1), we have shifted the lower limit of the macroscopic friction from $\mu _1$ to $\mu _0$, whereas $\mu _1$ becomes the upper limit of the creeping regime that centres around the viscous number of the creep regime, i.e. $J_c$. The proposed extension (6.1) recovers the original formulation (4.1) by choosing $J_c=0$ or $\mu _0 = \mu _1$. We explicitly note that we here aim to model the rheological behaviour for very small, but non-zero viscous numbers, i.e. $J\rightarrow 0$. This quasi-static, but still dynamic, regime might thus be different from the static case at $J=0$ (Perrin et al. 2019).

6.3. 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 $\mu _0 \in [0,1]$ and $J_c$ for all simulations conducted. To this end, we extend the range of $J$ to the full range observed in the simulations, i.e. $J \in [10^{-9},10^2]$. Since the extended formulation (6.1) is meant as an extension of the classical $\mu (J)$-rheology (4.1), we keep the values of the previously determined coefficients $\mu _1, \mu _2$ and $J_f$ 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 caption of this figure. 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 $J\in [10^{-9},10^{-5}]$. We also note that the curves of the extended model and the fit from § 5.2 (orange line) collapse for $J > 10^{-4}$, where the extension term $\mu ^c$ for the creep regime effectively evaluates to $\mu _1$ and thus reduces to the original model. Analysing the trend of the values determined for the two new parameters $\mu _0$ and $J_c$, we again notice a decrease in the friction coefficient $\mu _0$ with increasing polydispersity. This decrease, however, is less significant than before for $\mu _1$ and $\mu _2$ and a difference of only around 10 % can be seen between the monodisperse case and the one with strongest polydispersity. Generally, $\mu _0$ is approximately three times smaller than $\mu _1$. Determining the parameter $J_c$ faces similar challenges as discussed for $J_f$ 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 $J_f$ and confirms the physical meaning of $J_c$ discussed above to describe the average value of $J$ for the creep regime.

Figure 11. Macroscopic friction factor $\mu$ as function of viscous number $J$. The legend is as in figure 7. In addition, the fit of the extended model, (6.1), is shown in green. The insets show a magnified view for low values of $J$ using a linear $y$-axis. For cases (a) mono ($\mu _0 = 0.087, J_c = 1.39\times 10^{-6}$), (b) poly-10 ($\mu _0 = 0.082, J_c = 8.01\times 10^{-7}$), (c) poly-50 ($\mu _0 = 0.080, J_c = 1.44\times 10^{-6}$) and (d) poly-100 ($\mu _0 = 0.076, J_c = 1.55\times 10^{-6}$).

Due to the observed marginal sensitivity of $\mu _0$ and $J_c$ on the polydispersity, and the general difficulty of measurements for the creep regime, we do not attempt to express a functional dependence on $\phi _m$ 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 average of the fitted coefficients:

(6.2)$$\begin{gather} \mu_0 = 0.082 \pm 0.004 , \end{gather}$$

(6.3)$$\begin{gather}J_c = 1.30\times 10^{{-}6} \pm 2.91\times 10^{{-}7}. \end{gather}$$

We evaluate the performance of our creep-extended rheology model by computing the $R^2$ for the different empirical correlations over the entire range of $J\in [10^{-9},10^3]$. For that, we compare: (i) (4.1) with the parameters of Boyer et al. (2011); (ii) (4.1) with the parameters of Houssais et al. (2016); (iii) (6.1) with the parameters given in table 4 and figure 11; and (iv) (6.1) with the parameters given by correlations (5.2)–(5.5) and (6.2)–(6.3). The resulting $R^2$ values are given in table 6. In comparison with the existing model parameterizations of Boyer et al. (2011) and Houssais 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.

Table 6. The $R^2$ values of different parameterizations of the rheology model, (4.1) and (4.2), evaluated with respect to the simulated data for $J \in [10^{-9},10^2]$. The polydispersity extension developed in § 5 features fits with coefficients from table 4, and the correlations (5.2)–(5.5). The creep extension, (6.1), in the current section reuses these coefficients or correlations, respectively, and adds the coefficients from figure 11 for the individual fits or the correlations from (6.2)–(6.3). Bold indicates the most relevant data.

For completeness, we also show the $R^2$ values for $\phi (J)$ over the extended range of $J$, 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 $\phi (J)$ 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 $R^2$ values larger than 0.992 for all here considered cases for both, $\mu$ and $\phi$, and without any further calibration. This is a significant improvement on the previous rheology model and its parameterizations.

7. 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 3-D simulation domains using an efficiently coupled LBM–DEM to fully resolve all relevant scales in space and time. In particular, particle collisions are modelled 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 $26\,000$ 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 $\mu$ and the particle volume fraction $\phi$ as a function of the viscous number $J$, i.e. the $\mu (J)$-rheology. We compared our results with 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, $J\in [10^{-9},10^3]$, and the highly resolved data, we were able to enhance the $\mu (J)$-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 $\mu _1$ and $\mu _2$, and the free parameter $K_n$ as functions of $\phi _m$. The parameter $\phi _m$ 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 $\mu$, which are substantially lower than the frictional regime would suggest. This observation gave rise to the idea to enhance the $\mu (J)$-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 with 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 $\mu (J)$-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 behaviour of sheared sediment beds. Additionally, strongly inhomogeneous particle size distributions might require a more local description of the bed composition to accurately describe the rheology of the sediment. This highlights another benefit of our simulation approach, where such changes can be made with ease, allowing for efficient parametric studies.