## 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 Reference Frey and Church2011). The transport typically occurs along a slope or by a fluid flow shearing the sediment (Jerolmack & Daniels Reference Jerolmack and Daniels2019) and can lead to bedform evolution, such as ripples and dunes, even for laminar flow conditions (Lajeunesse *et al.* Reference Lajeunesse, Malverti, Lancien, Armstrong, Métivier, Coleman, Smith, Davies, Cantelli and Parker2010). 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.* Reference Aussillous, Chauchat, Pailha, Médale and Guazzelli2013; Houssais *et al.* Reference Houssais, Ortiz, Durian and Jerolmack2016; Kidanemariam Reference Kidanemariam2016; Vowinckel *et al.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021). 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 Reference Morris and Boulay1999; Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011). 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.* Reference Gallier, Lemaire, Peters and Lobry2014; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018; Vowinckel *et al.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021).

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 Reference Morris and Boulay1999). While Morris & Boulay (Reference Morris and Boulay1999) 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.* Reference Boyer, Guazzelli and Pouliquen2011; Dagois-Bohy *et al.* Reference Dagois-Bohy, Hormozi, Guazzelli and Pouliquen2015). For laminar viscous flows, i.e. a Stokes number $St=\rho _p\dot {\gamma }d^2_p/\eta _f$ smaller than 10 (Bagnold Reference Bagnold1954; Ness & Sun Reference Ness and Sun2016), where $\rho _p$ is the particle density and $d_p$ is the characteristic particle diameter, this measure allowed Boyer *et al.* (Reference Boyer, Guazzelli and Pouliquen2011) 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 (Reference Guazzelli and Pouliquen2018) 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.* Reference Aussillous, Chauchat, Pailha, Médale and Guazzelli2013; Maurin, Chauchat & Frey Reference Maurin, Chauchat and Frey2016; Vowinckel *et al.* Reference Vowinckel, Biegert, Luzzatto-Fegiz and Meiburg2019*a*). This analogy is important for two-phase fluid sediment transport modelling (Jenkins & Hanes Reference Jenkins and Hanes1998; Hsu, Jenkins & Liu Reference Hsu, Jenkins and Liu2004), where the fluid–particle mixture is treated as two separated continua with interconnected conservation laws of mass and momentum (Ouriemi, Aussillous & Guazzelli Reference Ouriemi, Aussillous and Guazzelli2009). The empirical correlations of the $\mu (J)$-rheology can provide the constitutive equations needed to close this set of equations (Chauchat *et al.* Reference Chauchat, Cheng, Nagel, Bonamy and Hsu2017; Lee & Huang Reference Lee and Huang2018; Lee Reference Lee2021). 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.* (Reference Boyer, Guazzelli and Pouliquen2011) 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.* (Reference Revil-Baudard, Chauchat, Hurther and Barraud2015), who investigated sheet-flow processes under turbulent flow conditions, that these correlations need adjustments for more dilute systems, whereas Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* (Reference Barker, Schaeffer, Bohórquez and Gray2015), who then proposed an extension to tackle this problem (Barker & Gray Reference Barker and Gray2017).

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 Reference Denn and Morris2014; Pähtz *et al.* Reference Pähtz, Durán, De Klerk, Govender and Trulsson2019). 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.* (Reference Boyer, Guazzelli and Pouliquen2011) reported $\mu _1=0.32$ and $\phi _m=0.585$ for the monodisperse case, but it has been shown by Tapia, Pouliquen & Guazzelli (Reference Tapia, Pouliquen and Guazzelli2019) 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 Reference Guazzelli and Pouliquen2018). 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 Reference Biegert, Vowinckel and Meiburg2017), 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 Reference Chang and Powell1994; Gondret & Petit Reference Gondret and Petit1997). 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.* Reference Boyer, Guazzelli and Pouliquen2011) to $\phi _m=0.64$. Consequently, models for $\phi _m$ in bidisperse volume-imposed rheometry were proposed by Dörr, Sadiki & Mehdizadeh (Reference Dörr, Sadiki and Mehdizadeh2013) and Mwasame, Wagner & Beris (Reference Mwasame, Wagner and Beris2016) that can also be applied to polydisperse systems (Pednekar, Chun & Morris Reference Pednekar, Chun and Morris2018).

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 Reference Trulsson, Andreotti and Claudin2012; Ness & Sun Reference Ness and Sun2016) 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.* (Reference Amarsid, Delenne, Mutabaruka, Monerie, Perales and Radjai2017) 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.* Reference Bauer2020*a*) 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 (Reference Rettinger and Rüde2017) and Rettinger & Rüde (Reference Rettinger and Rüde2020). This extends our pore-resolved simulations of fluid flow through porous media (Fattahi *et al.* Reference Fattahi, Waluga, Wohlmuth, Rüde, Manhart and Helmig2016; Gil *et al.* Reference Gil, Galache, Godenschwager and Rüde2017; Rybak *et al.* Reference Rybak, Schwarzmeier, Eggenweiler and Rüde2021), and is in line with previous erosion studies using a similar methodology (Derksen Reference Derksen2011; Rettinger *et al.* Reference Rettinger, Godenschwager, Eibl, Preclik, Schruff, Frings and Rüde2017). We follow the approach by Vowinckel *et al.* (Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021) 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 Reference Kuhnle1993; Frings Reference Frings2008). The rather large disparity of the grain sizes is achieved using the efficient parallelization scheme of Eibl & Rüde (Reference Eibl and Rüde2018). 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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), including the classical empirical correlations of the $\mu (J)$-rheology (Boyer *et al.* Reference Boyer, Guazzelli and Pouliquen2011). 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 (Reference Rettinger and Rüde2020). 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.* Reference Bauer2020*a*), 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.

### 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.* (Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017). For the present studies, we employ the $D3Q19$ two-relaxation-time model of Ginzburg, Verhaeghe & d'Humieres (Reference Ginzburg, Verhaeghe and d'Humieres2008). 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

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 Reference Wachs2019) 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 Reference Cundall and Strack1979). In our case, the normal and tangential collision components are given by a linear spring–dashpot model, similar to Costa *et al.* (Reference Costa, Boersma, Westerweel and Breugem2015) and Biegert *et al.* (Reference Biegert, Vowinckel and Meiburg2017). Following van der Hoef *et al.* (Reference van der Hoef, Ye, van Sint Annaland, Andrews, Sundaresan and Kuipers2006), 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.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021), and the collision time $T_c$. The latter is chosen according to the findings in Rettinger & Rüde (Reference Rettinger and Rüde2020) 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 (Reference Thornton, Cummins and Cleary2013), 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 (Reference Joseph and Hunt2004).

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 (Reference Rettinger and Rüde2020) 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 (Reference Aidun, Lu and Ding1998) 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.* (Reference Ginzburg, Verhaeghe and d'Humieres2008) 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.* (Reference Wen, Zhang, Tu, Wang and Fang2014). Following Ladd (Reference Ladd1994), 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.* (Reference Dorschner, Chikatamarla, Bösch and Karlin2015), using density and pressure tensor information from surrounding fluid cells and the particle's velocity.

As shown in Rettinger & Rüde (Reference Rettinger and Rüde2020), 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 Reference Nguyen and Ladd2002; Biegert *et al.* Reference Biegert, Vowinckel and Meiburg2017). Thus, the total hydrodynamic interaction force and torque on a particle $i$ is here computed as

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 (Reference Rettinger and Rüde2020). 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

*a*,

*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

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.

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.* Reference Boyer, Guazzelli and Pouliquen2011; Aussillous *et al.* Reference Aussillous, Chauchat, Pailha, Médale and Guazzelli2013) 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 $\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.

### 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.* Reference Biegert, Vowinckel and Meiburg2017; Jain, Vowinckel & Fröhlich Reference Jain, Vowinckel and Fröhlich2017). 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:

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.* (Reference Ouriemi, Aussillous, Medale, Peysson and Guazzelli2007), 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 Reference Bagnold1954). 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.* Reference Costa, Boersma, Westerweel and Breugem2015; Biegert *et al.* Reference Biegert, Vowinckel and Meiburg2017; Rettinger & Rüde Reference Rettinger and Rüde2017, Reference Rettinger and Rüde2020). 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 Reference Eibl and Rüde2018; Bauer, Köstler & Rüde Reference Bauer, Köstler and Rüde2020*b*). The details of our simulation approach are presented in Bauer *et al.* (Reference Bauer2020*a*). The approach has successfully been applied in previous large-scale studies of particle-resolved simulations (e.g. Rettinger *et al.* Reference Rettinger, Godenschwager, Eibl, Preclik, Schruff, Frings and Rüde2017; Götz *et al.* Reference Götz, Iglberger, Stürmer and Rüde2010), 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.* Reference Houssais, Ortiz, Durian and Jerolmack2016; Vowinckel *et al.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021). 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.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021).

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 Reference Kidanemariam and Uhlmann2014). Note that other authors have used different threshold values for this definition (Houssais *et al.* Reference Houssais, Ortiz, Durian and Jerolmack2016; Biegert *et al.* Reference Biegert, Vowinckel and Meiburg2017), 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.

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

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.

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.* Reference Vowinckel, Nikora, Kempe and Fröhlich2017, Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019*b*)

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.* Reference Ferdowsi, Ortiz, Houssais and Jerolmack2017). 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.* (Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021), 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.* Reference Madraki, Hormozi, Ovarlez, Guazzelli and Pouliquen2017). 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

This definition is in line with the one proposed by the two-phase model of Aussillous *et al.* (Reference Aussillous, Chauchat, Pailha, Médale and Guazzelli2013) and successfully used in the analysis of Vowinckel *et al.* (Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021). Note that we do not introduce an artificial confining pressure $P_0$ at the top wall as suggested by Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), 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 5*a*), 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.

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.

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 Reference Farr2013; Brouwers Reference Brouwers2014).

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 Reference Krieger and Dougherty1959; Morris & Boulay Reference Morris and Boulay1999; Stickel & Powell Reference Stickel and Powell2005; Guazzelli & Morris Reference Guazzelli and Morris2011). 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.* Reference Boyer, Guazzelli and Pouliquen2011; Dagois-Bohy *et al.* Reference Dagois-Bohy, Hormozi, Guazzelli and Pouliquen2015; Tapia *et al.* Reference Tapia, Pouliquen and Guazzelli2019). 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.* (Reference Boyer, Guazzelli and Pouliquen2011) followed the argument of Cassar, Nicolas & Pouliquen (Reference Cassar, Nicolas and Pouliquen2005) to show that the rheology of the fluid–particle mixture is governed by $J$. Based on these considerations, Boyer *et al.* (Reference Boyer, Guazzelli and Pouliquen2011) proposed the following empirical correlations as a rheological model, which became known as the $\mu (J)$-rheology and reads in its most general form:

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 (Reference Jop, Forterre and Pouliquen2005) and Cassar *et al.* (Reference Cassar, Nicolas and Pouliquen2005) 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.* (Reference Cassar, Nicolas and Pouliquen2005), $\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.* (Reference Boyer, Guazzelli and Pouliquen2011). 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 (Reference Einstein1905), and $K_n$ is a parameter that has been determined empirically by best fit to experimental data (Morris & Boulay Reference Morris and Boulay1999; Boyer *et al.* Reference Boyer, Guazzelli and Pouliquen2011).

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.* Reference Cassar, Nicolas and Pouliquen2005; Boyer *et al.* Reference Boyer, Guazzelli and Pouliquen2011; Houssais *et al.* Reference Houssais, Ortiz, Durian and Jerolmack2016) or $J_0$ (e.g. Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018; Vowinckel *et al.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021).

### 4.2. Existing model parameterizations

In the work of Boyer *et al.* (Reference Boyer, Guazzelli and Pouliquen2011), 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.* (Reference Boyer, Guazzelli and Pouliquen2011), whereas Morris & Boulay (Reference Morris and Boulay1999) 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) and Tapia *et al.* (Reference Tapia, Pouliquen and Guazzelli2019), which differ most notably in the range of measured $J$ values. In Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* (Reference Boyer, Guazzelli and Pouliquen2011).

In contrast, Tapia *et al.* (Reference Tapia, Pouliquen and Guazzelli2019) 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.* (Reference Tapia, Pouliquen and Guazzelli2019), 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.

### 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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), 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.* (Reference Boyer, Guazzelli and Pouliquen2011) and from Morris & Boulay (Reference Morris and Boulay1999).

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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* (Reference Boyer, Guazzelli and Pouliquen2011), normalized by $\phi _m=0.585$, and Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), 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.* (Reference Boyer, Guazzelli and Pouliquen2011) 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.* (Reference Boyer, Guazzelli and Pouliquen2011) in the extrapolated region of a more dilute regime (Vowinckel *et al.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021). 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.* Reference Boyer, Guazzelli and Pouliquen2011; Vowinckel *et al.* Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021), 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.* (Reference Singh, Mari, Denn and Morris2018) 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016). This range, however, was not addressed by Boyer *et al.* (Reference Boyer, Guazzelli and Pouliquen2011) nor Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.

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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016). 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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Morris & Boulay (Reference Morris and Boulay1999). 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.* (Reference Boyer, Guazzelli and Pouliquen2011), Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) and Morris & Boulay (Reference Morris and Boulay1999) 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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016). 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.* (Reference Boyer, Guazzelli and Pouliquen2011) 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), 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.

Comparing the case mono with Boyer *et al.* (Reference Boyer, Guazzelli and Pouliquen2011), 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016). 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 (Reference Morris and Boulay1999) and Boyer *et al.* (Reference Boyer, Guazzelli and Pouliquen2011). 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.* (Reference Pednekar, Chun and Morris2018) 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 Reference Farr2013; Brouwers Reference Brouwers2014). Alternatively, for our specific case of a truncated log-normal distribution, the relation

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 Reference Tripathi and Khakhar2011). 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$.

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

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.* (Reference Boyer, Guazzelli and Pouliquen2011), Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) and Morris & Boulay (Reference Morris and Boulay1999), 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

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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Morris & Boulay (Reference Morris and Boulay1999), i.e. $K_n= 1$ and $K_n= 0.75$, respectively, are not as accurate.

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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) across the entire range of $J$ (figures 6 and 7*a*), but also for all other cases considered (figure 7*c*,*e*,*g*).

Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* (Reference Bouzid, Izzet, Trulsson, Clément, Claudin and Andreotti2015), 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.* (Reference Bouzid, Izzet, Trulsson, Clément, Claudin and Andreotti2015), no burst-like behaviour can be observed in figure 9. On the contrary, Bouzid *et al.* (Reference Bouzid, Izzet, Trulsson, Clément, Claudin and Andreotti2015) 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 Reference Kamrin and Koval2012).

Recently, Gillissen & Ness (Reference Gillissen and Ness2020) 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.* (Reference Boyer, Guazzelli and Pouliquen2011), 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 (Reference Gillissen and Ness2020) 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 (Reference Gillissen and Ness2020) 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 Reference Gillissen and Ness2020). 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.* (Reference Boyer, Guazzelli and Pouliquen2011) 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.* (Reference Cassar, Nicolas and Pouliquen2005) and Jop *et al.* (Reference Jop, Forterre and Pouliquen2005), 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:

In comparison with the original model of Boyer *et al.* (Reference Boyer, Guazzelli and Pouliquen2011), (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.* Reference Perrin, Clavaud, Wyart, Metzger and Forterre2019).

### 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.

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:

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.* (Reference Boyer, Guazzelli and Pouliquen2011); (ii) (4.1) with the parameters of Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016); (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.* (Reference Boyer, Guazzelli and Pouliquen2011) and Houssais *et al.* (Reference Houssais, Ortiz, Durian and Jerolmack2016), 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 $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.* (Reference Boyer, Guazzelli and Pouliquen2011) 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.* (Reference Houssais, Ortiz, Durian and Jerolmack2016) 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.* Reference Boyer, Guazzelli and Pouliquen2011). 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.* (Reference Jop, Forterre and Pouliquen2005). 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.

## Supplementary material

Supplementary material is available at https://doi.org/10.1017/jfm.2021.870.

## Acknowledgements

The authors thank M. Houssais for sharing her data and gratefully acknowledge the Erlangen Regional Computing Center (www.rrze.fau.de/) as well as the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on their supercomputers. We thank three anonymous reviewers for their valuable comments that helped to improve the manuscript.

## 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.

## Data availability statement

The complete simulation data used for the analysis is made available on Zenodo via https://doi.org/10.5281/zenodo.5529671.

## Author contributions

B.V. conceived the original idea. C.R. performed the simulations and implemented the data analysis. C.R. and B.V. contributed equally in analysing the data, developing the model extension, reaching conclusion and in writing the paper. S.E. implemented functionalities essential for polydisperse set-ups, assisted in the planning of the post-processing steps, and helped shape the research. U.R. supervised the project. All authors reviewed the final manuscript.

## 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.* Reference Ferdowsi, Ortiz, Houssais and Jerolmack2017). 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.* Reference Rosato, Strandburg, Prinz and Swendsen1987).

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 centre of mass position above $h^{top} = 15.5 \bar {d}_p$ belong to the bed's top region, which is roughly $2\bar {d}_p$ below the average sediment bed height $\langle h_b \rangle _t$ (cf. table 2). We then sort these $N_p^{top}$ topmost particles according to their diameters into bins of size $\bar {d}_p / 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. $t \in [0,12\,000]\, t_{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 behaviour. 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 $t > 6000\,t_{ref}$. This indicates that the fast segregation process, as described by Ferdowsi *et al.* (Reference Ferdowsi, Ortiz, Houssais and Jerolmack2017), 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. We additionally verified that all lower layers are hardly affected by the segregation process since they exhibit a diameter distribution that is close to the one of the entire sediment bed, given in figure 3. This allows us to regard the majority of the bed as a homogeneous packing, simplifying the analysis.