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

(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 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

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

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

(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. (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üde2020b). The details of our simulation approach are presented in Bauer et al. (Reference Bauer2020a). 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.

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

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

(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. (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 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 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.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:

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

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. (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).

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. (Reference Boyer, Guazzelli and Pouliquen2011) and Houssais et al. (Reference Houssais, Ortiz, Durian and Jerolmack2016). Additionally, curves of (4.1) and (4.2) are shown, parameterized as 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) (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. (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.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. (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.

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

(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 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$.

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

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

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.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 7a), but also for all other cases considered (figure 7c,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).

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

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 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:

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

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

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.