4. Extending models of mixture viscosity

To improve the mixture rheological model based on the characterization of each phase, we propose to split the apparent viscosity of the mixture phase, as the sum of the fluid and granular contributions, following Boyer et al. (2011) [12]. In particular, using previous definition of the different viscosities, one has

$$ \Sigma_{xy}^{\text{tot}} \equiv \eta_{\text{eq}}^{m} \dot{\gamma} = \Sigma_{xy}^{f} + \Sigma_{xy}^{p} \equiv \eta_{\text{eq}}^{f} \dot{\gamma} + \mu_{\text{eq}}^{p} P^{p}, $$

leading to

$$ \eta_{\text{eq}}^{m} = \eta_{\text{eq}}^{f} + \mu_{\text{eq}}^{p} P^{p} / \dot{\gamma}. $$

Equivalently, the effective friction is

$$ \mu_{\text{eq}}^{m} = \eta_{\text{eq}}^{f} \dot{\gamma} / P^{p} + \mu_{\text{eq}}^{p}, $$

using the granular pressure $P^p$ as the only relevant pressure for a Coulomb-type model of the mixture phase. Then, using the correlations obtained previously for the equivalent fluid phase (see Fig. 8) and for the equivalent granular phase, respectively, one proposes

$$ \frac{\eta_{\text{eq}}^{M}}{\eta_{f}} = \underbrace{(1 - \phi)^{-\beta}}_{\eta_{\text{eq}}^{f} / \eta_{f}} + \underbrace{\mu_{\text{eq}}^{p}(\phi) \mathcal{F}(\phi, \text{Re}_{p})}_{\eta_{\text{eq}}^{p} / \eta_{f}}, \qquad (55) $$

with $\beta = \beta_0(1 + l_g^2 \rho_f \dot{\gamma} / \eta_f) = \beta_0(1 + (l_g/d)^2 \text{Re}_p)$ as obtained in Sec. IV B 1 and $\mu_{\text{eq}}^p(\phi)$ from Eq. (39) with $\phi_c = 0.615$, $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.75$, and $\sqrt{K_0} = 0.18$ as obtained in Sec. IV B 2. In Eq. (55), $\mathcal{F}$ models $P^p / (\eta_f \dot{\gamma})$, which corresponds to the friction-to-viscosity relationship, or equivalently the dimensionless granular pressure. It shall be noted that $\mathcal{F} \equiv \mathcal{J}^{-1}$ as long as $\dot{\gamma} = \dot{\gamma}^p$. It will be assumed to be mostly the case in the following for simplicity. Moreover, if the flow is purely viscous then $\mathcal{F} \equiv \mathcal{J}^{-1} = K^{-1}$. For inertial situations, $\mathcal{F}$ could be more complex as $\mathcal{J} \neq K$ and then it is no longer the relevant dimensionless parameter to characterize the granular phase [see Eqs. (40) and (41)]. Then, $\mathcal{F}$ shall be explicitly written as a function of both $\phi$ and the local $\text{Re}_l^p$, as long as $\dot{\gamma} = \dot{\gamma}^p$ [see Eq. (41)]. It then depends on the local value of $\text{Re}_l^p$, and we propose here to simplify by giving an estimation based on the control parameter $\text{Re}_p$, as

$$ \mathcal{F}(\phi, \text{Re}_p) \equiv \mathcal{J}^{-1} = (1 + \alpha r \text{Re}_p) \left( \frac{b\phi}{\phi_c - \phi} \right)^2. \qquad (56) $$

The associated model for the effective viscosity is labeled as $\eta_{\text{eq}}^F$. In Fig. 16, we compare the effective viscosity obtained from the present IBM-DEM simulations $\eta_{\text{eq}}^m$ with $\eta_{\text{eq}}^F$ [Fig. 16(c)]. For comparison, the models of Refs. [49] and [12], labeled as $\eta_{\text{eq}}^{KD}$ and $\eta_{\text{eq}}^{BGP}$, respectively, are shown in Figs. 16(a) and 16(b), respectively. For clarity, all the models used in Fig. 16 are summarized in Table II. The model of Ref. [12] corresponds to $\text{Re}_p = 0$ in Eq. (56) and a different fluid contribution, namely, $\eta_{\text{eq}}^{BGP} / \eta_f = 1 + \frac{5}{2}\phi(1 - \phi/\phi_c)^{-1} + \mu_{\text{eq}}^p(\phi)\mathcal{F}(\phi, \text{Re}_p = 0)$. Note that equivalent models have been recently proposed leading to similar results as the latter model [27,57]. Disregarding too large $\phi$, i.e., large viscosity, such approach clearly improves the collapse of the different cases ($\text{Re}_p, \theta$) considered in this study, when using model (56) (as dispersion of the data clearly diminishes in [Fig. 16(c)]). This is not necessarily the case for other models as in [Fig. 16(a)] or [Fig. 16(b)]. Nevertheless, the present model [Fig. 16(c)] seems to fail at perfectly capturing the expected trend at small $\phi$ (solid black line), even if all data collapse on a same curve.

To understand previous observations, we report $(\eta_{\text{eq}}^p / \eta_f) / \mu_{\text{eq}}^p \equiv P^p / (\eta_f \dot{\gamma})$, i.e., $\mathcal{F}$ as a function of $\phi$ in Fig. 15. It is shown to depend only weakly on $\text{Re}_p$. Yet, $\mathcal{J}^{-1} = K^{-1}$, i.e., $\mathcal{F}(\phi, \text{Re}_p = 0)$ as modeled by Eq. (56), fails at predicting the evolution of $\mathcal{F}$ at small $\phi$ (blue line). This suggests that the local inertia of the grains $\text{Re}_l^p$ strongly depends on the state variable $\phi$ and more weakly on $\text{Re}_p$.
<table><caption>TABLE II. Rheological models assessed in Fig. 16. Here, $\mu_{eq}^p(\phi)$ is given by (39) with $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.75$ and $\sqrt{K_0} = 0.18$ and $\beta = \beta_0 + a r^{2/3} \text{Re}_p \text{Ga}^{-4/3}$ with $\beta_0 = 0.565$ and $a = 0.85$.</caption><thead><tr><th>Name</th><th>Reference</th><th>Apparent viscosity of the fluid phase</th><th>Apparent viscosity of the granular phase</th><th>Apparent viscosity of the mixture phase</th><th></th></tr></thead><tbody><tr><td>KD</td><td>[49]</td><td>N/A</td><td>N/A</td><td>$\frac{\eta_{eq}^{KD}}{\eta_f} = (1 - \frac{\phi}{\phi_m})^{-\frac{5}{2}\phi_m}$</td><td>$\phi_m = 0.62$</td></tr><tr><td>BGP</td><td>[12]</td><td>$\frac{\eta_{eq}^f}{\eta_f} = 1 + \frac{5}{2}\phi(1 - \frac{\phi}{\phi_c})^{-1}$</td><td>$\frac{\eta_{eq}^p}{\eta_f} = \mu_{eq}^p(\phi)(\frac{\phi}{\phi_c - \phi})^2$</td><td>$\frac{\eta_{eq}^{BGP}}{\eta_f} = \frac{\eta_{eq}^f}{\eta_f} + \frac{\eta_{eq}^p}{\eta_f}$</td><td>$\phi_c = 0.615$</td></tr><tr><td>F</td><td></td><td>$\frac{\eta_{eq}^f}{\eta_f} = (1 - \phi)^{-\beta}$</td><td>$\frac{\eta_{eq}^p}{\eta_f} = \mu_{eq}^p(\phi)(1 + \alpha r \text{Re}_p)(\frac{b\phi}{\phi_c - \phi})^2$</td><td>$\frac{\eta_{eq}^F}{\eta_f} = \frac{\eta_{eq}^f}{\eta_f} + \frac{\eta_{eq}^p}{\eta_f}$</td><td>$\phi_c = 0.615$<br>$\alpha = 0.635$</td></tr><tr><td>FLBC</td><td></td><td>$\frac{\eta_{eq}^f}{\eta_f} = (1 - \phi)^{-\beta}$</td><td>$\frac{\eta_{eq}^p}{\eta_f} = \mu_{eq}^p(\phi)(\frac{b\phi}{\phi_c - \phi})^2 + (\frac{\phi}{\phi_c - \phi})^2$</td><td>$\frac{\eta_{eq}^{FLBC}}{\eta_f} = \frac{\eta_{eq}^f}{\eta_f} + \frac{\eta_{eq}^p}{\eta_f}$</td><td>$\phi_c = 0.615$<br>$b = 0.7$</td></tr></tbody></table>

As $r \text{Re}_p^p = \mathcal{I}^2/\mathcal{J}$, then beyond $\mathcal{K}$, the ratio $\mathcal{I}^2/\mathcal{J}$ is also mostly a function of $\phi$. We obtain from our data

$$ \mathcal{F}(\phi) \equiv \mathcal{J}^{-1} \equiv \mathcal{K}^{-1} \left(1 + \alpha \frac{\mathcal{I}^2}{\mathcal{J}}\right) = \left(\frac{b\phi}{\phi_c - \phi}\right)^2 + \left(\frac{\phi}{\phi_c - \phi}\right). \quad (57) $$

Note that the last term in Eq. (57) fits the term $\alpha \mathcal{I}^2 \mathcal{J}^{-1} \mathcal{K}^{-1}$. This new model is shown in Fig. 15 (green solid line) to predict more accurately the trend of $\mathcal{F}$ on the entire range of $\phi$. It shall be noted that the influence of the average length scale $h_g/d$, which could affect the nature of the diffuse interface between the top of the granular bed and the pure fluid layer, does not seem to significantly modify the relevance of model (57) as long as $\phi > 0.05$ (see Fig. 22 in Appendix B).

Using the latter model for $\mathcal{F}$ (57) leads to $\eta_{eq}^{FLBC}$, shown in Fig. 16(d). This is shown to improve the prediction of the mixture viscosity at moderate $\phi$, i.e., for $\eta_{eq}^m/\eta_f \lesssim 10$, to capture the diffuse interface. It can also be used to predict the effective friction $\mu_{eq}^m$ as shown in Fig. 13 (solid red line).

V. SUMMARY OF THE MAIN FINDINGS

The bedload transport of a granular bed by a laminar Couette flow, has been studied numerically using a fully resolved IBM-DEM approach, referred here to as the microscale description. Bedload

<figure><img src="image_17.png" /><figcaption>FIG. 15. $(\eta_{eq}^p/\eta_f)/\mu_{eq}^p \equiv \mathcal{F}$ as a function of $\phi$ for all the cases considered in the present work: blue solid line corresponds to model (56) for $\text{Re}_p = 0$ and green solid line is Eq. (57) (see text for details).</figcaption></figure>
<table><caption>TABLE II. Rheological models assessed in Fig. 16. Here, $\mu_{eq}^p(\phi)$ is given by (39) with $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.75$ and $\sqrt{K_0} = 0.18$ and $\beta = \beta_0 + a r^{2/3} \text{Re}_p \text{Ga}^{-4/3}$ with $\beta_0 = 0.565$ and $a = 0.85$.</caption><thead><tr><th>Name</th><th>Reference</th><th>Apparent viscosity of the fluid phase</th><th>Apparent viscosity of the granular phase</th><th>Apparent viscosity of the mixture phase</th><th></th></tr></thead><tbody><tr><td>KD</td><td>[49]</td><td>N/A</td><td>N/A</td><td>$\frac{\eta_{eq}^{KD}}{\eta_f} = (1 - \frac{\phi}{\phi_m})^{-\frac{5}{2}\phi_m}$</td><td>$\phi_m = 0.62$</td></tr><tr><td>BGP</td><td>[12]</td><td>$\frac{\eta_{eq}^f}{\eta_f} = 1 + \frac{5}{2}\phi(1 - \frac{\phi}{\phi_c})^{-1}$</td><td>$\frac{\eta_{eq}^p}{\eta_f} = \mu_{eq}^p(\phi)(\frac{\phi}{\phi_c - \phi})^2$</td><td>$\frac{\eta_{eq}^{BGP}}{\eta_f} = \frac{\eta_{eq}^f}{\eta_f} + \frac{\eta_{eq}^p}{\eta_f}$</td><td>$\phi_c = 0.615$</td></tr><tr><td>F</td><td></td><td>$\frac{\eta_{eq}^f}{\eta_f} = (1 - \phi)^{-\beta}$</td><td>$\frac{\eta_{eq}^p}{\eta_f} = \mu_{eq}^p(\phi)(1 + \alpha r \text{Re}_p)(\frac{b\phi}{\phi_c - \phi})^2$</td><td>$\frac{\eta_{eq}^F}{\eta_f} = \frac{\eta_{eq}^f}{\eta_f} + \frac{\eta_{eq}^p}{\eta_f}$</td><td>$\phi_c = 0.615$<br>$\alpha = 0.635$</td></tr><tr><td>FLBC</td><td></td><td>$\frac{\eta_{eq}^f}{\eta_f} = (1 - \phi)^{-\beta}$</td><td>$\frac{\eta_{eq}^p}{\eta_f} = \mu_{eq}^p(\phi)(\frac{b\phi}{\phi_c - \phi})^2 + (\frac{\phi}{\phi_c - \phi})$</td><td>$\frac{\eta_{eq}^{FLBC}}{\eta_f} = \frac{\eta_{eq}^f}{\eta_f} + \frac{\eta_{eq}^p}{\eta_f}$</td><td>$\phi_c = 0.615$<br>$b = 0.7$</td></tr></tbody></table>

As $r \text{Re}_p^p = \mathcal{I}^2/\mathcal{J}$, then beyond $\mathcal{K}$, the ratio $\mathcal{I}^2/\mathcal{J}$ is also mostly a function of $\phi$. We obtain from our data

$$ \mathcal{F}(\phi) \equiv \mathcal{J}^{-1} \equiv \mathcal{K}^{-1} \left(1 + \alpha \frac{\mathcal{I}^2}{\mathcal{J}}\right) = \left(\frac{b\phi}{\phi_c - \phi}\right)^2 + \left(\frac{\phi}{\phi_c - \phi}\right). \quad (57) $$

Note that the last term in Eq. (57) fits the term $\alpha \mathcal{I}^2 \mathcal{J}^{-1} \mathcal{K}^{-1}$. This new model is shown in Fig. 15 (green solid line) to predict more accurately the trend of $\mathcal{F}$ on the entire range of $\phi$. It shall be noted that the influence of the average length scale $h_g/d$, which could affect the nature of the diffuse interface between the top of the granular bed and the pure fluid layer, does not seem to significantly modify the relevance of model (57) as long as $\phi > 0.05$ (see Fig. 22 in Appendix B).

Using the latter model for $\mathcal{F}$ (57) leads to $\eta_{eq}^{FLBC}$, shown in Fig. 16(d). This is shown to improve the prediction of the mixture viscosity at moderate $\phi$, i.e., for $\eta_{eq}^m/\eta_f \lesssim 10$, to capture the diffuse interface. It can also be used to predict the effective friction $\mu_{eq}^m$ as shown in Fig. 13 (solid red line).

## V. SUMMARY OF THE MAIN FINDINGS

The bedload transport of a granular bed by a laminar Couette flow, has been studied numerically using a fully resolved IBM-DEM approach, referred here to as the microscale description. Bedload

<figure><img src="image_17.png" /><figcaption>FIG. 15. $(\eta_{eq}^p/\eta_f)/\mu_{eq}^p \equiv \mathcal{F}$ as a function of $\phi$ for all the cases considered in the present work: blue solid line corresponds to model (56) for $\text{Re}_p = 0$ and green solid line is Eq. (57) (see text for details).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 16. Comparison of the predicted apparent viscosity of the mixture versus that obtained from present simulations $\eta_{eq}^m/\eta_f$: (a) Krieger and Dougherty (1959) $\eta_{eq}^{KD}/\eta_f$; (b) Boyer *et al.* (2011) $\eta_{eq}^{BGP}/\eta_f$; (c) model (55) and (56) $\eta_{eq}^F/\eta_f$; (d) model (55)–(57) $\eta_{eq}^{FLBC}/\eta_f$: see Table I for detail of symbols. The two dashed lines represents the discrepancy ratio of 2 and 0.5, respectively.</figcaption></figure>

transport is known for decades to be strongly dependent on the Shields number $\theta$. Above a critical Shields number $\theta_c$ at which the granular bed starts moving, increasing $\theta$ leads to more significant granular transport. As $\theta \to \theta_c$ from above, this therefore tends to a singular sharp interface of the granular bed with only a thin layer of the order of the grain diameter moving. However, when $\theta$ increases, the interface between the static granular bed and the pure fluid phase becomes more diffuse and is characterized by a thicker moving layer. Accordingly, the vertical distribution of the solid fraction $\phi$ evolves from a nearly singular one at the bed surface when $\theta \approx \theta_c$ towards a smoother shape when $\theta$ increases. Yet, in bedload transport, the upper bed surface hardly gets rid of some sharp transition or discontinuity.

Accordingly, the strategy to model bedload transport at a scale larger than the grain, referred to as mesoscale, remains uncertain, while this is obviously unmissable when dealing with large scale systems. To unify physical processes of the granular suspension system from weak transport to intense transport, we have proposed here to provide a mixture approach, for which the entire system, including the granular phase and the fluid phase, is solved using a single equivalent fluid. For that purpose, an upscaling—or averaging—method from microscale simulation results has been adopted to provide an independent description of each phase at the mesoscale prior describing the mixture phase. This approach leads to a diffuse interface at the upper surface of the granular bed whatever the Shields number $\theta$.

A. Volume-average homogenization remains relevant even in a system presenting a significant shear at the micrograin scale

When upscaling, separation between the grain size $d$ and the averaging scale $h_g$ is required for the validity of the asymptotic development at the mesoscale, which reads $h_g \gg d/2$. We have shown that this constraint has to be released to capture the high velocity gradient at the bed surface, as for instance shown for the velocity profiles in Fig. 20 with $h_g$ in the range $[d/10, d]$. Obviously, the invariance of the flow in the horizontal plane and in time allowed us to decrease $h_g$ to capture gradients, but only down to a certain point. For too small $h_g$, mesoscale quantities becomes noisy
<figure><img src="image_18.png" /><figcaption>FIG. 16. Comparison of the predicted apparent viscosity of the mixture versus that obtained from present simulations $\eta_{eq}^m/\eta_f$: (a) Krieger and Dougherty (1959) $\eta_{eq}^{KD}/\eta_f$; (b) Boyer *et al.* (2011) $\eta_{eq}^{BGP}/\eta_f$; (c) model (55) and (56) $\eta_{eq}^F/\eta_f$; (d) model (55)–(57) $\eta_{eq}^{FLBC}/\eta_f$: see Table I for detail of symbols. The two dashed lines represents the discrepancy ratio of 2 and 0.5, respectively.</figcaption></figure>

transport is known for decades to be strongly dependent on the Shields number $\theta$. Above a critical Shields number $\theta_c$ at which the granular bed starts moving, increasing $\theta$ leads to more significant granular transport. As $\theta \to \theta_c$ from above, this therefore tends to a singular sharp interface of the granular bed with only a thin layer of the order of the grain diameter moving. However, when $\theta$ increases, the interface between the static granular bed and the pure fluid phase becomes more diffuse and is characterized by a thicker moving layer. Accordingly, the vertical distribution of the solid fraction $\phi$ evolves from a nearly singular one at the bed surface when $\theta \approx \theta_c$ towards a smoother shape when $\theta$ increases. Yet, in bedload transport, the upper bed surface hardly gets rid of some sharp transition or discontinuity.

Accordingly, the strategy to model bedload transport at a scale larger than the grain, referred to as mesoscale, remains uncertain, while this is obviously unmissable when dealing with large scale systems. To unify physical processes of the granular suspension system from weak transport to intense transport, we have proposed here to provide a mixture approach, for which the entire system, including the granular phase and the fluid phase, is solved using a single equivalent fluid. For that purpose, an upscaling—or averaging—method from microscale simulation results has been adopted to provide an independent description of each phase at the mesoscale prior describing the mixture phase. This approach leads to a diffuse interface at the upper surface of the granular bed whatever the Shields number $\theta$.

A. Volume-average homogenization remains relevant even in a system presenting a significant shear at the micrograin scale

When upscaling, separation between the grain size $d$ and the averaging scale $h_g$ is required for the validity of the asymptotic development at the mesoscale, which reads $h_g \gg d/2$. We have shown that this constraint has to be released to capture the high velocity gradient at the bed surface, as for instance shown for the velocity profiles in Fig. 20 with $h_g$ in the range $[d/10, d]$. Obviously, the invariance of the flow in the horizontal plane and in time allowed us to decrease $h_g$ to capture gradients, but only down to a certain point. For too small $h_g$, mesoscale quantities becomes noisy
and strong oscillations along the shear direction in the granular phase appear. This is actually the limit of validity of the mesoscale approach: microscopic structures—granular arrangement—are not filtered out by this mesoscale, breaking the isotropic condition for equations validity at a scale smaller than the mesoscale. This highlights that the grain size and the macroscopic length scale of deformation are of the same order of magnitude, in bedload transport.

The influence of $h_g$ on velocity profiles, and particularly the granular profiles, have been shown to be quite significant, i.e., on the length scale of the diffuse mesoscale interface (Fig. 20). Yet, and more surprisingly, the apparent viscosity of the granular phase seems to be relatively independent of $h_g$, as shown in Fig. 21 (besides an increase of the dispersion of the data when $h_g$ is decreased). Diffusing the bed interface with volume average thus keeps the right momentum balance within the diffuse interface (Sec. IV A). Then, choosing setting $h_g = d$ gives a reasonable description of the flow at the mesoscale, and has therefore been used to discuss the rheological models for the different phases. This is supported by results discussed in Appendix B.

## B. Weak inertia $Re_p$ and shear localization modify mixture rheological models

Regarding the second question, we have shown that the effective viscosity $\eta_{eq}^f$ of the equivalent fluid evolves with the granular volume fraction $\phi$ within the bed, with in particular $\eta_{eq}^f$ increasing with $\phi$. This is in line with previous studies on homogeneous suspension flows. However, $\eta_{eq}^f$ is found to be smaller than that found for homogeneous neutrally buoyant suspensions. Moreover, we observe a dependency of this effective viscosity with $\theta$ and $Re_p$. These have been discussed as a consequence of the very specific nature of the mixture flow in bedload application in Sec. IV B 1. Weight of the grain remains crucial to describe such system and the associated dissipation in the equivalent fluid phase. In particular, it leads to both sharpening of bed interface in the limit of relatively small $\theta$ and a vertical agitation of the grain influencing fluid dissipation at finite inertia $Re_p$. Combinations of these two effects in bedload hardly allow to reach assumptions required to obtain more classical viscosity model for homogeneous neutrally buoyant suspension. A model extended from homogeneous and buoyant systems has therefore been proposed to describe the fluid-phase effective viscosity Eqs. (51) and (52) for such configuration including its specific dynamics. This model was shown to capture $\eta_{eq}^f(\phi)$ in the range of parameters covered here, and therefore including both sharp interface effect and weak inertia. For finite inertia $Re_p$, it was shown that this model combined $Re_p$ and a characteristic dimensionless length scale $l_g/d$ to be associated with fluid structure induced by vertical agitation of the granular phase, by analogy to eddy viscosity in turbulent models.

However, the granular phase is shown to be well described by the classical $\mu$-I rheology of the granular effective friction, now extended to the $\mu$-K rheology for immersed granular flows. The parameters which characterize this rheology, are found to be independent of $\theta$ and $Re_p$.

Finally, in line with Ref. [12], we discuss the model of mixture viscosity as the sum of that of the fluid and granular phases, respectively. Accordingly, as long as the model for each individual phases found before are used, the main input on the mixture viscosity modeling concerns the friction-to-viscosity law $F$ (55) required to set the granular viscosity law from a granular friction law. Assuming no difference between fluid shear and granular shear, this terms corresponds to the so-called viscous number $J$ in granular rheology, i.e., $F \equiv J^{-1}$. Based on that assumption, we found that inertia has to be accounted for to extract its relevant dependency with $\phi$ close to the bed surface (Sec. IV B 4 and Fig. 15). This leads to an extra order $\phi$ term in the mixture viscosity in the small $\phi$ limit (i.e., close to the bed surface) emanating from the granular phase (i.e., contacts), unlike homogeneous neutrally buoyant suspension, for which this term only comes from the fluid contribution, and is usually referred to as equivalent to the Einstein's viscosity. Surprisingly the effective viscosity of the mixture phase in the limit of small $\phi$ becomes in the present configuration closer to the Einstein's viscosity than the equivalent fluid phase was. Moreover, the newly proposed
and strong oscillations along the shear direction in the granular phase appear. This is actually the limit of validity of the mesoscale approach: microscopic structures—granular arrangement—are not filtered out by this mesoscale, breaking the isotropic condition for equations validity at a scale smaller than the mesoscale. This highlights that the grain size and the macroscopic length scale of deformation are of the same order of magnitude, in bedload transport.

The influence of $h_g$ on velocity profiles, and particularly the granular profiles, have been shown to be quite significant, i.e., on the length scale of the diffuse mesoscale interface (Fig. 20). Yet, and more surprisingly, the apparent viscosity of the granular phase seems to be relatively independent of $h_g$, as shown in Fig. 21 (besides an increase of the dispersion of the data when $h_g$ is decreased). Diffusing the bed interface with volume average thus keeps the right momentum balance within the diffuse interface (Sec. IV A). Then, choosing setting $h_g = d$ gives a reasonable description of the flow at the mesoscale, and has therefore been used to discuss the rheological models for the different phases. This is supported by results discussed in Appendix B.

## B. Weak inertia $Re_p$ and shear localization modify mixture rheological models

Regarding the second question, we have shown that the effective viscosity $\eta_{eq}^f$ of the equivalent fluid evolves with the granular volume fraction $\phi$ within the bed, with in particular $\eta_{eq}^f$ increasing with $\phi$. This is in line with previous studies on homogeneous suspension flows. However, $\eta_{eq}^f$ is found to be smaller than that found for homogeneous neutrally buoyant suspensions. Moreover, we observe a dependency of this effective viscosity with $\theta$ and $Re_p$. These have been discussed as a consequence of the very specific nature of the mixture flow in bedload application in Sec. IV B 1. Weight of the grain remains crucial to describe such system and the associated dissipation in the equivalent fluid phase. In particular, it leads to both sharpening of bed interface in the limit of relatively small $\theta$ and a vertical agitation of the grain influencing fluid dissipation at finite inertia $Re_p$. Combinations of these two effects in bedload hardly allow to reach assumptions required to obtain more classical viscosity model for homogeneous neutrally buoyant suspension. A model extended from homogeneous and buoyant systems has therefore been proposed to describe the fluid-phase effective viscosity Eqs. (51) and (52) for such configuration including its specific dynamics. This model was shown to capture $\eta_{eq}^f(\phi)$ in the range of parameters covered here, and therefore including both sharp interface effect and weak inertia. For finite inertia $Re_p$, it was shown that this model combined $Re_p$ and a characteristic dimensionless length scale $l_g/d$ to be associated with fluid structure induced by vertical agitation of the granular phase, by analogy to eddy viscosity in turbulent models.

However, the granular phase is shown to be well described by the classical $\mu$-I rheology of the granular effective friction, now extended to the $\mu$-K rheology for immersed granular flows. The parameters which characterize this rheology, are found to be independent of $\theta$ and $Re_p$.

Finally, in line with Ref. [12], we discuss the model of mixture viscosity as the sum of that of the fluid and granular phases, respectively. Accordingly, as long as the model for each individual phases found before are used, the main input on the mixture viscosity modeling concerns the friction-to-viscosity law $F$ (55) required to set the granular viscosity law from a granular friction law. Assuming no difference between fluid shear and granular shear, this terms corresponds to the so-called viscous number $J$ in granular rheology, i.e., $F \equiv J^{-1}$. Based on that assumption, we found that inertia has to be accounted for to extract its relevant dependency with $\phi$ close to the bed surface (Sec. IV B 4 and Fig. 15). This leads to an extra order $\phi$ term in the mixture viscosity in the small $\phi$ limit (i.e., close to the bed surface) emanating from the granular phase (i.e., contacts), unlike homogeneous neutrally buoyant suspension, for which this term only comes from the fluid contribution, and is usually referred to as equivalent to the Einstein's viscosity. Surprisingly the effective viscosity of the mixture phase in the limit of small $\phi$ becomes in the present configuration closer to the Einstein's viscosity than the equivalent fluid phase was. Moreover, the newly proposed
<figure><img src="image_19.png" /><figcaption>FIG. 17. Vertical profiles of streamwise granular velocity (filled symbols), fluid velocity (open symbols), and solid volume fraction (solid lines). (●) case A; (■) case B; (▲) case C; (◆) case D.</figcaption></figure>

model for the mixture viscosity is found to capture nicely mixture viscosity laws obtained from micro-scale simulations in the range of $\theta$ and $Re_p$ considered here.

## ACKNOWLEDGMENTS

We thank Prof. Alessandro Bottaro for helpful discussions and suggestions. This work was supported by the Chaire Bioskins of the IdEx of the University of Toulouse, of which A. Bottaro is the holder. We also thank Annaig Pedrono for her support regarding the solvers JADIM and GraDyM. This work was performed using HPC resources from CALMIP (Grant 2018-[P1027]).

## APPENDIX A: EFFECT OF THE SPATIAL RESOLUTION AND THE DOMAIN SIZE ON THE MESOSCALE RESULTS

Here, we assess the effect of the spatial resolution, the domain size and the “thickness” of the immersed boundary $\Delta_\alpha$ (see Sec. II A for definition). Four cases are considered and presented in Table III. The first one is the numerical setup used in the main text of the present work, the three others correspond to a larger domain, a higher spatial resolution and a thinner immersed boundary $\Delta_\alpha$. Note that for all these cases, we set as characteristic size of the weighting function $h_g = d$ (see Sec. III B).

Figures 17 and 18 present the mesoscopic velocity profiles and constitutive laws, respectively, for the various phases. The results of the four cases are roughly identical. As for the velocity profiles, results are superimposed except at the top of the moving bedload layer. In this region, the granular

<table><caption>TABLE III. Summary of the different cases considered in Appendix A. Here, $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22).</caption><thead><tr><th>Case</th><th>Spatial resolution</th><th>Computational domain size</th><th>Thickness of the immersed boundary</th><th>Comment</th></tr><tr><td></td><td><i>d</i> / Δ<i>x</i></td><td><i>L</i><sub>x</sub> × <i>L</i><sub>y</sub> × <i>L</i><sub>z</sub></td><td>Δ<sub>α</sub></td><td></td></tr></thead><tbody><tr><td>A (●)</td><td>10</td><td>4<i>d</i> × 20<i>d</i> × 4<i>d</i></td><td>1.3Δ<i>x</i> = 0.13<i>d</i></td><td>reference</td></tr><tr><td>B (■)</td><td>10</td><td>8<i>d</i> × 20<i>d</i> × 8<i>d</i></td><td>1.3Δ<i>x</i> = 0.13<i>d</i></td><td>larger domain</td></tr><tr><td>C (▲)</td><td>20</td><td>4<i>d</i> × 20<i>d</i> × 4<i>d</i></td><td>2.6Δ<i>x</i> = 0.13<i>d</i></td><td>refined grid</td></tr><tr><td>D (◆)</td><td>20</td><td>4<i>d</i> × 20<i>d</i> × 4<i>d</i></td><td>1.3Δ<i>x</i> = 0.065<i>d</i></td><td>reduced IBM thickness</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>FIG. 17. Vertical profiles of streamwise granular velocity (filled symbols), fluid velocity (open symbols), and solid volume fraction (solid lines). (●) case A; (■) case B; (▲) case C; (◆) case D.</figcaption></figure>

model for the mixture viscosity is found to capture nicely mixture viscosity laws obtained from micro-scale simulations in the range of $\theta$ and $Re_p$ considered here.

## ACKNOWLEDGMENTS

We thank Prof. Alessandro Bottaro for helpful discussions and suggestions. This work was supported by the Chaire Bioskins of the IdEx of the University of Toulouse, of which A. Bottaro is the holder. We also thank Annaig Pedrono for her support regarding the solvers JADIM and GraDyM. This work was performed using HPC resources from CALMIP (Grant 2018-[P1027]).

## APPENDIX A: EFFECT OF THE SPATIAL RESOLUTION AND THE DOMAIN SIZE ON THE MESOSCALE RESULTS

Here, we assess the effect of the spatial resolution, the domain size and the “thickness” of the immersed boundary $\Delta_\alpha$ (see Sec. II A for definition). Four cases are considered and presented in Table III. The first one is the numerical setup used in the main text of the present work, the three others correspond to a larger domain, a higher spatial resolution and a thinner immersed boundary $\Delta_\alpha$. Note that for all these cases, we set as characteristic size of the weighting function $h_g = d$ (see Sec. III B).

Figures 17 and 18 present the mesoscopic velocity profiles and constitutive laws, respectively, for the various phases. The results of the four cases are roughly identical. As for the velocity profiles, results are superimposed except at the top of the moving bedload layer. In this region, the granular

<table><caption>TABLE III. Summary of the different cases considered in Appendix A. Here, $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22).</caption><thead><tr><th>Case</th><th>Spatial resolution</th><th>Computational domain size</th><th>Thickness of the immersed boundary</th><th>Comment</th></tr><tr><th></th><th>d/Δx</th><th>L<sub>x</sub> × L<sub>y</sub> × L<sub>z</sub></th><th>Δ<sub>α</sub></th><th></th></tr></thead><tbody><tr><td>A (●)</td><td>10</td><td>4d × 20d × 4d</td><td>1.3Δx = 0.13d</td><td>reference</td></tr><tr><td>B (■)</td><td>10</td><td>8d × 20d × 8d</td><td>1.3Δx = 0.13d</td><td>larger domain</td></tr><tr><td>C (▲)</td><td>20</td><td>4d × 20d × 4d</td><td>2.6Δx = 0.13d</td><td>refined grid</td></tr><tr><td>D (◆)</td><td>20</td><td>4d × 20d × 4d</td><td>1.3Δx = 0.065d</td><td>reduced IBM thickness</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>FIG. 18. Effect of the spatial resolution and size of the domain on the constitutive laws of the equivalent (a) fluid phase, (b), (c) granular phase, and (d), (e) mixture phase: (black) case A; (magenta) case B; (red) case C; (blue) case D.</figcaption></figure>

velocity is slightly faster for cases A and B than for cases C and D. The maximum difference is $0.11\bar{\Gamma}d$ which corresponds to about 19% of the maximum velocity of case A.

As for the constitutive laws, the results are mostly superimposed. The apparent viscosity of the fluid phase slightly depends on mesh refinement [see Fig. 18(a) for $\phi \le 0.5$. At higher solid volume fraction, variations induced by the averaging procedure occur and interpretation of the results becomes more difficult. A best-fit of the results using the model (32), namely $\eta_f^{eq}/\eta_f = (1-\phi)^{-\beta}$, gives $\beta = 1.71$ and 1.62 for cases A and B, respectively, and 1.73 and 1.71 for cases C and D, respectively (note that here we only use data for $\phi \le 0.5$).

As for the equivalent granular phase [Figs. 18(b) and 18(c)], fluctuations are observed at small $K$ where the bed is static. Results are roughly similar between the four cases. The same conclusion can be drawn for the equivalent mixture phase [Fig. 18(d)].

Figure 19 compares the vertical distribution of volumetric forces acting on the equivalent granular phase for cases A and B. The results are mostly superimposed except in the static zone $y \le h_{\text{bed}} - 4d$. In this zone, the amplitude of the oscillations is observed to decrease as the domain's size is increased. More quantitatively, when the size is doubled in both horizontal directions, the oscillations amplitude is reduced by about 30% in frame (a) and 60% in frame (b). This supports the argument that these oscillations are due to the somewhat limited number of static grains on which we apply the coarse-graining approach. As they remain at the exact same location, the time averaging of these correlated samples leads to such oscillations. We verified that this effect is not problematic for the results presented in the main text, since it corresponds to very low values of $K \le 10^{-5}$ which are not considered here, and a very limited range of the solid volume fraction $\phi \ge 0.601$, close to that of jamming.

Overall, Figs. 17–19 allow us to conclude that the results presented in the main text are not significantly influenced by the size of the domain, the spatial resolution, and the “thickness” of the immersed boundary. In addition, they provide an explanation about the oscillations in the static region of the granular bed observed in the vertical profiles.
<figure><img src="image_20.png" /><figcaption>FIG. 18. Effect of the spatial resolution and size of the domain on the constitutive laws of the equivalent (a) fluid phase, (b), (c) granular phase, and (d), (e) mixture phase: (black) case A; (magenta) case B; (red) case C; (blue) case D.</figcaption></figure>

velocity is slightly faster for cases A and B than for cases C and D. The maximum difference is $0.11\bar{\Gamma}d$ which corresponds to about 19% of the maximum velocity of case A.

As for the constitutive laws, the results are mostly superimposed. The apparent viscosity of the fluid phase slightly depends on mesh refinement [see Fig. 18(a) for $\phi \le 0.5$. At higher solid volume fraction, variations induced by the averaging procedure occur and interpretation of the results becomes more difficult. A best-fit of the results using the model (32), namely $\eta_f^{eq}/\eta_f = (1-\phi)^{-\beta}$, gives $\beta = 1.71$ and 1.62 for cases A and B, respectively, and 1.73 and 1.71 for cases C and D, respectively (note that here we only use data for $\phi \le 0.5$).

As for the equivalent granular phase [Figs. 18(b) and 18(c)], fluctuations are observed at small $K$ where the bed is static. Results are roughly similar between the four cases. The same conclusion can be drawn for the equivalent mixture phase [Fig. 18(d)].

Figure 19 compares the vertical distribution of volumetric forces acting on the equivalent granular phase for cases A and B. The results are mostly superimposed except in the static zone $y \le h_{\text{bed}} - 4d$. In this zone, the amplitude of the oscillations is observed to decrease as the domain's size is increased. More quantitatively, when the size is doubled in both horizontal directions, the oscillations amplitude is reduced by about 30% in frame (a) and 60% in frame (b). This supports the argument that these oscillations are due to the somewhat limited number of static grains on which we apply the coarse-graining approach. As they remain at the exact same location, the time averaging of these correlated samples leads to such oscillations. We verified that this effect is not problematic for the results presented in the main text, since it corresponds to very low values of $K \le 10^{-5}$ which are not considered here, and a very limited range of the solid volume fraction $\phi \ge 0.601$, close to that of jamming.

Overall, Figs. 17–19 allow us to conclude that the results presented in the main text are not significantly influenced by the size of the domain, the spatial resolution, and the “thickness” of the immersed boundary. In addition, they provide an explanation about the oscillations in the static region of the granular bed observed in the vertical profiles.
<figure><img src="image_21.png" /><figcaption>FIG. 19. Comparison of vertical distribution of all the terms of the momentum balance Eqs. (42) and (44) of the equivalent granular phase in the (a) vertical direction, (b) streamwise direction, for the cases A (pale) and B (bright). Same legend as Figs. 5(a) and 6(a). See Eqs. (46) and (48).</figcaption></figure>

## APPENDIX B: INFLUENCE OF THE CHARACTERISTIC LENGTH OF THE WEIGHTING FUNCTION $h_g$ ON THE MESOSCOPIC RESULTS

In this Appendix, we assess the dependency of the mesoscopic results on the characteristic length of the weighting function $h_g$. Here, $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III.

The equations of the averaged phases (18) and (19) are supposed to be valid only if $h_g \gg d$. Conversely, one needs a $h_g$ not too large to capture the local variations of the mesoscopic quantities. Therefore, it is of crucial importance to assess the effect of $h_g$ on the results.

Figure 20 presents the streamwise velocity profiles for both phases and for the solid volume fraction for $h_g = d$, $d/2$, $d/4$ and $d/10$. Note that $h_g = d/10$ here corresponds to $h_g = \Delta x$, that is the smallest value of $h_g$ which can be considered here.

The velocity profiles strongly depend on the value of $h_g$. In particular, the velocity of the granular phase in the moving bedload layer is quite sensitive to the specific value of $h_g$. For instance, at $y = h_{\text{bed}} + d$, the granular velocity decreases by about 40% from $h_g/d = 1/10$ to $h_g/d = 1$. The velocity difference between the fluid phase and the granular phase is about $4V_s$ for $h_g/d = 1$ and almost zero for $h_g/d = 1/10$ [see the inset of Fig. 20(a)]. The influence of $h_g$ on the solid volume fraction $\phi$ is even more dramatic, as shown in Fig. 20(b). For $h_g/d < 1/2$, the fluctuations of $\phi$ are of the same order of magnitude as $\phi$ itself. Note that the observed oscillations have a period of $d$.

<figure><img src="image_22.png" /><figcaption>FIG. 20. (a) Velocity profile of the fluid phase (dashed) and the granular phase (solid) for various of $h_g = \{d/10, d/4, d/2, d\}$. Velocity is scaled by $V_s$ with $V_s = (\rho_p - \rho_f)gd^2/(18\eta_f)$ the Stokes velocity. (b) Solid volume fraction for the same values of $h_g$. Here, $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIG. 19. Comparison of vertical distribution of all the terms of the momentum balance Eqs. (42) and (44) of the equivalent granular phase in the (a) vertical direction, (b) streamwise direction, for the cases A (pale) and B (bright). Same legend as Figs. 5(a) and 6(a). See Eqs. (46) and (48).</figcaption></figure>

## APPENDIX B: INFLUENCE OF THE CHARACTERISTIC LENGTH OF THE WEIGHTING FUNCTION $h_g$ ON THE MESOSCOPIC RESULTS

In this Appendix, we assess the dependency of the mesoscopic results on the characteristic length of the weighting function $h_g$. Here, $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III.

The equations of the averaged phases (18) and (19) are supposed to be valid only if $h_g \gg d$. Conversely, one needs a $h_g$ not too large to capture the local variations of the mesoscopic quantities. Therefore, it is of crucial importance to assess the effect of $h_g$ on the results.

Figure 20 presents the streamwise velocity profiles for both phases and for the solid volume fraction for $h_g = d$, $d/2$, $d/4$ and $d/10$. Note that $h_g = d/10$ here corresponds to $h_g = \Delta x$, that is the smallest value of $h_g$ which can be considered here.

The velocity profiles strongly depend on the value of $h_g$. In particular, the velocity of the granular phase in the moving bedload layer is quite sensitive to the specific value of $h_g$. For instance, at $y = h_{bed} + d$, the granular velocity decreases by about 40% from $h_g/d = 1/10$ to $h_g/d = 1$. The velocity difference between the fluid phase and the granular phase is about $4V_s$ for $h_g/d = 1$ and almost zero for $h_g/d = 1/10$ [see the inset of Fig. 20(a)]. The influence of $h_g$ on the solid volume fraction $\phi$ is even more dramatic, as shown in Fig. 20(b). For $h_g/d < 1/2$, the fluctuations of $\phi$ are of the same order of magnitude as $\phi$ itself. Note that the observed oscillations have a period of $d$.

<figure><img src="image_22.png" /><figcaption>FIG. 20. (a) Velocity profile of the fluid phase (dashed) and the granular phase (solid) for various of $h_g = \{d/10, d/4, d/2, d\}$. Velocity is scaled by $V_s$ with $V_s = (\rho_p - \rho_f)gd^2/(18\eta_f)$ the Stokes velocity. (b) Solid volume fraction for the same values of $h_g$. Here, $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIG. 21. Apparent rheology of the mesoscopic fluid phase (top row), granular phase (second and third rows) and mixture phase (bottom row): (a, blue) $h_g = d/2$, (b, magenta) $h_g = d/4$; (c, red) $h_g = d/10$. Results for $h_g = d$ are also plotted as reference (black). For the apparent viscosity of the fluid phase, the model $\eta_{\text{eq}}^f/\eta_f = (1-\phi)^{-1.71}$ is also plotted (solid black line). Here, $\text{Re}_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III.</figcaption></figure>

Figure 21 presents the constitutive laws of the fluid phase $\eta_f(\phi)$, the granular phase $\mu_p(K) - \phi(K)$ and the mixture phase $\eta_m(\phi)$ for $h_g = d/2$, $d/4$ and $d/10$, respectively (see Sec. IV B for more details and definitions). Results obtained for $h_g = d$ are also plotted on each graph as reference.

One can see that, as in Fig. 20(b), fluctuations and dispersion appear when $h_g$ decreases. Nevertheless, the dispersed data seems to remained centered around the results obtained at large $h_g$.

Note also that the rheology of the mixture phase for $\phi < 0.3$ is well captured even for $h_g = d/10$. This part of the rheology correspond to the moving bedload layer of particles. In this flow region, temporal averaging is efficient in reducing spatial fluctuations; furthermore, a significant part of the total viscosity comes from the fluid which was shown to be less sensitive to the value of $h_g$. For $\phi = 0.3$ and $h_g = d$, we obtain a total viscosity of $2.67\eta_f$ and an apparent viscosity of the fluid phase of $1.84\eta_f$, that is almost 70% of the total viscosity comes from the fluid contribution.
<figure><img src="image_23.png" /><figcaption>FIG. 21. Apparent rheology of the mesoscopic fluid phase (top row), granular phase (second and third rows) and mixture phase (bottom row): (a, blue) $h_g = d/2$, (b, magenta) $h_g = d/4$; (c, red) $h_g = d/10$. Results for $h_g = d$ are also plotted as reference (black). For the apparent viscosity of the fluid phase, the model $\eta_{\text{eq}}^f/\eta_f = (1-\phi)^{-1.71}$ is also plotted (solid black line). Here, $\text{Re}_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III.</figcaption></figure>

Figure 21 presents the constitutive laws of the fluid phase $\eta_f(\phi)$, the granular phase $\mu_p(K) - \phi(K)$ and the mixture phase $\eta_m(\phi)$ for $h_g = d/2$, $d/4$ and $d/10$, respectively (see Sec. IV B for more details and definitions). Results obtained for $h_g = d$ are also plotted on each graph as reference.

One can see that, as in Fig. 20(b), fluctuations and dispersion appear when $h_g$ decreases. Nevertheless, the dispersed data seems to remained centered around the results obtained at large $h_g$.

Note also that the rheology of the mixture phase for $\phi < 0.3$ is well captured even for $h_g = d/10$. This part of the rheology correspond to the moving bedload layer of particles. In this flow region, temporal averaging is efficient in reducing spatial fluctuations; furthermore, a significant part of the total viscosity comes from the fluid which was shown to be less sensitive to the value of $h_g$. For $\phi = 0.3$ and $h_g = d$, we obtain a total viscosity of $2.67\eta_f$ and an apparent viscosity of the fluid phase of $1.84\eta_f$, that is almost 70% of the total viscosity comes from the fluid contribution.
<figure><img src="image_24.png" /><figcaption>FIG. 22. $\mathcal{F} \equiv P^p/(\eta_f \dot{\gamma})$ as a function of $\phi$ for $h_g = d$ (black symbols), $h_g = d/2$ (blue symbols), $h_g = d/4$ (magenta symbols), $h_g = d/10$ (red symbols). As $d/4$ and $d/10$ are shown to be dispersed (see Fig. 21), $h_g = d$ and $h_g = d/2$ symbols are slightly larger for clarity. Here, $\text{Re}_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III. Dash line and solid are reported from Fig. 15, i.e., from Ref. [12] and the extended model (57).</figcaption></figure>

Similarly, the influence of $h_g$ on the friction-to-viscosity parameter $\mathcal{F}$ is shown in Fig. 22. Clearly, $h_g$ affects the obtained results at small $\phi$, which indicates that the extended model (57) includes the influence of the transition from the granular bed to the pure fluid layer. However, this separation is mostly visible for $\phi < 0.1$, even lower. Then, it clearly shows that the model (57) is, at least, required to capture the rapid transition observed in $0.05 < \phi < 0.4$ whatever $h_g$.

## APPENDIX C: DRAG FORCE ON THE PARTICLES

In this Appendix, we present the results for the drag force acting on the particles within the granular bed. To this end, we assume that the term $n \langle f_1 \rangle_x^p$ present in Eqs. (44) and (48) contains the drag force only, that is we neglect the added-mass force. This is relevant in the present work since we verified that the local acceleration of the apparent fluid and granular phase, namely $D_f \langle \mathbf{u} \rangle^f / Dt$ and $D_p \langle \mathbf{u} \rangle^p / Dt$ were negligible.

Richardson and Zaki (1954) [58] proposed a drag law for a suspension of mono-dispersed spherical particles in an infinite medium which is valid for arbitrary particle Reynolds numbers, which reads

$$ \mathbf{F}_D = \frac{(\rho_p - \rho_f)\phi g}{V_t(1-\phi)^{l-2}} (\langle \mathbf{u} \rangle^f - \langle \mathbf{u} \rangle^p), \quad (C1) $$

where $V_t$ is the terminal velocity of the particle and $l$ is an exponent which depends on the particle Reynolds number. In the present range of $\text{Re}_p$, one can assume that (i) the terminal velocity follows the Stokes law, i.e., $V_t = V_s = (\rho_p - \rho_f)gd^2/(18\eta_f)$; (ii) the value of $l$ is roughly constant and we set $l = 4.65$. The mean drag force, after Ref. [58], reads

$$ \underbrace{\frac{\mathbf{F}_D}{\eta_f \dot{\Gamma}/d}}_{(F^*)} = \underbrace{\frac{18\phi}{(1-\phi)^{2.65}}}_{(\Phi^*)} \underbrace{\frac{(\langle \mathbf{u} \rangle^f - \langle \mathbf{u} \rangle^p)}{\dot{\Gamma}d}}_{(\Delta U^*)}. \quad (C2) $$

Figure 23 assesses the validity of Eq. (C2) on one configuration, namely $\text{Re}_p = 1$, $\theta = 0.67$, $Ga = 1.22$, by showing either $\Delta U^*$ versus $F^*/\Phi^*$, or $F^*/\Delta U^*$ versus $\Phi^*$ for various values of the characteristic length of the weighting function $h_g$. Note that in our problem, we have $\mathbf{F}_D = n \langle f_1 \rangle_x^p$. The validity of Eq. (C2) would correspond to a perfect superimposition of the two curves. A better agreement is observed as $h_g$ is decreased; however, the fluctuations of the force also increase in the static region within the granular bed. Conversely, for large $h_g$, the slip velocity $\langle \mathbf{u} \rangle^f - \langle \mathbf{u} \rangle^p$ is artificially increased, in line with Fig. 20(a). This leads to a larger discrepancy between the curves in Fig. 23. Nevertheless, the present results suggest that Richardson and Zaki's (1954) drag law may be a good approximation for predicting the drag force on grains which are transported as bedload.
<figure><img src="image_24.png" /><figcaption>FIG. 22. $\mathcal{F} \equiv P^p/(\eta_f \dot{\gamma})$ as a function of $\phi$ for $h_g = d$ (black symbols), $h_g = d/2$ (blue symbols), $h_g = d/4$ (magenta symbols), $h_g = d/10$ (red symbols). As $d/4$ and $d/10$ are shown to be dispersed (see Fig. 21), $h_g = d$ and $h_g = d/2$ symbols are slightly larger for clarity. Here, $\text{Re}_p = 1$ and $\theta = 0.67$ (Ga = 1.22) and we use the setup of case A of Table III. Dash line and solid are reported from Fig. 15, i.e., from Ref. [12] and the extended model (57).</figcaption></figure>

Similarly, the influence of $h_g$ on the friction-to-viscosity parameter $\mathcal{F}$ is shown in Fig. 22. Clearly, $h_g$ affects the obtained results at small $\phi$, which indicates that the extended model (57) includes the influence of the transition from the granular bed to the pure fluid layer. However, this separation is mostly visible for $\phi < 0.1$, even lower. Then, it clearly shows that the model (57) is, at least, required to capture the rapid transition observed in $0.05 < \phi < 0.4$ whatever $h_g$.

## APPENDIX C: DRAG FORCE ON THE PARTICLES

In this Appendix, we present the results for the drag force acting on the particles within the granular bed. To this end, we assume that the term $n \langle f_1 \rangle_x^p$ present in Eqs. (44) and (48) contains the drag force only, that is we neglect the added-mass force. This is relevant in the present work since we verified that the local acceleration of the apparent fluid and granular phase, namely $D_f \langle \mathbf{u} \rangle^f / Dt$ and $D_p \langle \mathbf{u} \rangle^p / Dt$ were negligible.

Richardson and Zaki (1954) [58] proposed a drag law for a suspension of mono-dispersed spherical particles in an infinite medium which is valid for arbitrary particle Reynolds numbers, which reads

$$ \mathbf{F}_D = \frac{(\rho_p - \rho_f)\phi g}{V_t (1 - \phi)^{l-2}} (\langle \mathbf{u} \rangle^f - \langle \mathbf{u} \rangle^p), \quad (C1) $$

where $V_t$ is the terminal velocity of the particle and $l$ is an exponent which depends on the particle Reynolds number. In the present range of $\text{Re}_p$, one can assume that (i) the terminal velocity follows the Stokes law, i.e., $V_t = V_s = (\rho_p - \rho_f)gd^2/(18\eta_f)$; (ii) the value of $l$ is roughly constant and we set $l = 4.65$. The mean drag force, after Ref. [58], reads

$$ \underbrace{\frac{\mathbf{F}_D}{\eta_f \dot{\Gamma}/d}}_{(F^*)} = \underbrace{\frac{18\phi}{(1-\phi)^{2.65}}}_{(\Phi^*)} \underbrace{\frac{\langle \mathbf{u} \rangle^f - \langle \mathbf{u} \rangle^p}{\dot{\Gamma} d}}_{(\Delta U^*)}. \quad (C2) $$

Figure 23 assesses the validity of Eq. (C2) on one configuration, namely $\text{Re}_p = 1$, $\theta = 0.67$, $Ga = 1.22$, by showing either $\Delta U^*$ versus $F^*/\Phi^*$, or $F^*/\Delta U^*$ versus $\Phi^*$ for various values of the characteristic length of the weighting function $h_g$. Note that in our problem, we have $\mathbf{F}_D = n \langle f_1 \rangle_x^p$. The validity of Eq. (C2) would correspond to a perfect superimposition of the two curves. A better agreement is observed as $h_g$ is decreased; however, the fluctuations of the force also increase in the static region within the granular bed. Conversely, for large $h_g$, the slip velocity $\langle \mathbf{u} \rangle^f - \langle \mathbf{u} \rangle^p$ is artificially increased, in line with Fig. 20(a). This leads to a larger discrepancy between the curves in Fig. 23. Nevertheless, the present results suggest that Richardson and Zaki's (1954) drag law may be a good approximation for predicting the drag force on grains which are transported as bedload.
<figure><img src="image_25.png" /><figcaption>FIG. 23. Vertical profile of the drag force on the grains compared with Richardson and Zaki's drag law [58] given in Eq. (C2) ($Re_p = 1, \theta = 0.67, Ga = 1.22$). Left: $\Delta U^*$ vs $F^*/\Phi^*$. Right: $F^*/\Delta U^*$ vs $\Phi^*$ with $F^* = \frac{F_p}{\eta_f \Gamma / d}$, $\Delta U^* = \frac{(\mathbf{u})^f - (\mathbf{u})^p}{\Gamma/d}$ and $\Phi^* = \frac{18\phi}{(1-\phi)^{2.65}}$. Each row corresponds to one value of the characteristic length of the weighting function $h_g$, namely (from top to bottom) $h_g = \{d, d/2, d/4, d/10\}$.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>FIG. 23. Vertical profile of the drag force on the grains compared with Richardson and Zaki's drag law [58] given in Eq. (C2) ($Re_p = 1, \theta = 0.67, Ga = 1.22$). Left: $\Delta U^*$ vs $F^*/\Phi^*$. Right: $F^*/\Delta U^*$ vs $\Phi^*$ with $F^* = \frac{F_p}{\eta_f \Gamma / d}$, $\Delta U^* = \frac{(\mathbf{u})^f - (\mathbf{u})^p}{\Gamma / d}$ and $\Phi^* = \frac{18\phi}{(1-\phi)^{2.65}}$. Each row corresponds to one value of the characteristic length of the weighting function $h_g$, namely (from top to bottom) $h_g = \{d, d/2, d/4, d/10\}$.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>FIG. 24. Effect of including a lubrication force on the phases rheologies: (a) Apparent viscosity $\eta_{eq}^f/\eta_f$ of the equivalent fluid phase vs solid volume fraction $\phi$. (b) Apparent friction coefficient $\mu_{eq}^p$ of the equivalent granular phase vs $\mathcal{J}$. (c) Apparent viscosity $\eta_{eq}^m/\eta_f$ of the equivalent mixture phase vs $\mathcal{J}$: present IBM-DEM simulations for (◇) $\text{Re}_p = 1$, $\theta = 0.67$, $Ga = 1.2$ and (*) $\text{Re}_p = 1$, $\theta = 0.33$, $Ga = 1.7$, with a lubrication force (small symbols) and without lubrication force (large symbols).</figcaption></figure>

## APPENDIX D: EFFECT OF INCLUDING A LUBRICATION FORCE ON THE MESOSCALE RESULTS

In this Appendix, we assess the effect of including a normal lubrication force in the forces applied to the grains (4). Here, only the normal component of the lubrication force is considered between particle i and j of velocity $\mathbf{u}_{pi}$ and $\mathbf{u}_{pj}$ and radius $R_i$ and $R_j$, respectively, and is denoted $\mathbf{f}_{lub}$. This force here reads [59]

$$ \mathbf{f}_{lub} = -\frac{6\pi \mu (\mathbf{u}_{pi} \cdot \mathbf{n} - \mathbf{u}_{pj} \cdot \mathbf{n})}{\delta_n + \eta_e} \left(\frac{R_i R_j}{R_i + R_j}\right)^2 \mathbf{n}, \quad (D1) $$

where $\mathbf{n}$ is the unit normal vector parallel to the grains centers, $\delta_n$ is the minimum distance between the grain surfaces, $\eta_e$ is an effective roughness height which is purposely added in Eq. (D1) to mimic real particles and avoid the divergence of the force when contact occurs ($\delta_n = 0$). In the following, we set $\eta_e/R = 4 \times 10^{-4}$. The present lubrication force is switched on when the distance between particles is such as $0 \le \delta_n \le R/2$. More details about the present modeling of the lubrication force can be found in Ref. [41].

Figure 24 illustrates the effect of including a lubrication force on the rheology of the fluid phase, the granular phase and the mixture phase, respectively, for one value of the Reynolds number
<figure><img src="image_26.png" /><figcaption>FIG. 24. Effect of including a lubrication force on the phases rheologies: (a) Apparent viscosity $\eta_{eq}^f/\eta_f$ of the equivalent fluid phase vs solid volume fraction $\phi$. (b) Apparent friction coefficient $\mu_{eq}^p$ of the equivalent granular phase vs $\mathcal{J}$. (c) Apparent viscosity $\eta_{eq}^m/\eta_f$ of the equivalent mixture phase vs $\mathcal{J}$: present IBM-DEM simulations for (◇) $\text{Re}_p = 1$, $\theta = 0.67$, $Ga = 1.2$ and (*) $\text{Re}_p = 1$, $\theta = 0.33$, $Ga = 1.7$, with a lubrication force (small symbols) and without lubrication force (large symbols).</figcaption></figure>

## APPENDIX D: EFFECT OF INCLUDING A LUBRICATION FORCE ON THE MESOSCALE RESULTS

In this Appendix, we assess the effect of including a normal lubrication force in the forces applied to the grains (4). Here, only the normal component of the lubrication force is considered between particle i and j of velocity $\mathbf{u}_{pi}$ and $\mathbf{u}_{pj}$ and radius $R_i$ and $R_j$, respectively, and is denoted $\mathbf{f}_{lub}$. This force here reads [59]

$$ \mathbf{f}_{lub} = -\frac{6\pi \mu (\mathbf{u}_{pi} \cdot \mathbf{n} - \mathbf{u}_{pj} \cdot \mathbf{n})}{\delta_n + \eta_e} \left(\frac{R_i R_j}{R_i + R_j}\right)^2 \mathbf{n}, \quad (D1) $$

where $\mathbf{n}$ is the unit normal vector parallel to the grains centers, $\delta_n$ is the minimum distance between the grain surfaces, $\eta_e$ is an effective roughness height which is purposely added in Eq. (D1) to mimic real particles and avoid the divergence of the force when contact occurs ($\delta_n = 0$). In the following, we set $\eta_e/R = 4 \times 10^{-4}$. The present lubrication force is switched on when the distance between particles is such as $0 \le \delta_n \le R/2$. More details about the present modeling of the lubrication force can be found in Ref. [41].

Figure 24 illustrates the effect of including a lubrication force on the rheology of the fluid phase, the granular phase and the mixture phase, respectively, for one value of the Reynolds number
(Re<sub>p</sub> = 1) and two values of the Shields number (θ = 0.33 and 0.67). The friction coefficient of the granular phase [Fig. 24(b)] is observed to be roughly unchanged when adding a lubrication force. The apparent viscosity of the fluid phase and the mixture phase is slightly smaller when lubrication force is included [Figs. 24(a) and 24(c)], while the friction coefficient of the mixture phase seems to be slightly larger [Fig. 24(d), θ = 0.67]. However, the variation of these parameters due to the addition of a lubrication force is observed to be smaller than that due to the Shields number θ. It is therefore reasonable to conclude that in the present problem and the present range of parameters, the effect of adding a normal lubrication force is marginal.

[1] J. S. Bridge, *Rivers and Floodplains: Forms, Processes, and Sedimentary Record* (John Wiley & Sons, New York, NY, 2009).

[2] R. Soulsby, *Dynamics of Marine Sands* (Thomas Telford, London, UK, 1997).

[3] B. Andreotti, Y. Forterre, and O. Pouliquen, *Les milieux granulaires: Entre fluide et solide* (EDP Sciences, Les Ulis, France, 2012).

[4] É. Guazzelli and O. Pouliquen, *Rheology of dense granular suspensions*, J. Fluid Mech. **852**, P1 (2018).

[5] R. Jackson, *The Dynamics of Fluidized Particles* (Cambridge University Press, Cambridge, UK, 2000).

[6] R. Jackson, Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid, *Chem. Eng. Sci.* **52**, 2457 (1997).

[7] Y. Kim, Z. Cheng, T.-J. Hsu, and J. Chauchat, A numerical study of sheet flow under monochromatic nonbreaking waves using a free surface resolving Eulerian two-phase flow model, *J. Geophys. Res. [Oceans]* **123**, 4693 (2018).

[8] A. Einstein, Eine neue bestimmung der moleküldimensionen (in german), *Annalen der Physik* **324**, 289 (1906).

[9] A. Einstein, Correction of the paper 'eine neue bestimmung der moleküldimensionen (1906) (in german), *Annalen der Physik* **339**, 591 (1911).

[10] G. Batchelor and J. Green, The determination of the bulk stress in a suspension of spherical particles to order $c^2$, *J. Fluid Mech.* **56**, 401 (1972).

[11] L. Gibilaro, K. Gallucci, R. Di Felice, and P. Pagliai, On the apparent viscosity of a fluidized bed, *Chem. Eng. Sci.* **62**, 294 (2007).

[12] F. Boyer, E. Guazzelli, and O. Pouliquen, Unifying suspension and granular rheology, *Phys. Rev. Lett.* **107**, 188301 (2011).

[13] M. Trulsson, B. Andreotti, and P. Claudin, Transition from the viscous to inertial regime in dense suspensions, *Phys. Rev. Lett.* **109**, 118305 (2012).

[14] G. Ovarlez, F. Bertrand, and S. Rodts, Local determination of the constitutive law of a dense suspension of noncolloidal particles through magnetic resonance imaging, *J. Rheol.* **50**, 259 (2006).

[15] P. Aussillous, J. Chauchat, M. Pailha, M. Médale, and E. Guazzelli, Investigation of the mobile granular layer in bedload transport by laminar shearing flows, *J. Fluid Mech.* **736**, 594 (2013).

[16] M. Houssais, C. P. Ortiz, D. J. Durian, and D. J. Jerolmack, Rheology of sediment transported by a laminar flow, *Phys. Rev. E* **94**, 062609 (2016).

[17] M. Houssais and D. J. Jerolmack, Toward a unifying constitutive relation for sediment transport across environments, *Geomorphology* **277**, 251 (2017).

[18] R. Maurin, J. Chauchat, B. Chareyre, and P. Frey, A minimal coupled fluid-discrete element model for bedload transport, *Phys. Fluids* **27**, 113302 (2015).

[19] R. Maurin, J. Chauchat, and P. Frey, Dense granular flow rheology in turbulent bedload transport, *J. Fluid Mech.* **804**, 490 (2016).

[20] T. Pähtz and O. Durán, Universal friction law at granular solid-gas transition explains scaling of sediment transport load with excess fluid shear stress, *Phys. Rev. Fluids* **3**, 104302 (2018).
(Re$_p$ = 1) and two values of the Shields number ($\theta$ = 0.33 and 0.67). The friction coefficient of the granular phase [Fig. 24(b)] is observed to be roughly unchanged when adding a lubrication force. The apparent viscosity of the fluid phase and the mixture phase is slightly smaller when lubrication force is included [Figs. 24(a) and 24(c)], while the friction coefficient of the mixture phase seems to be slightly larger [Fig. 24(d), $\theta$ = 0.67]. However, the variation of these parameters due to the addition of a lubrication force is observed to be smaller than that due to the Shields number $\theta$. It is therefore reasonable to conclude that in the present problem and the present range of parameters, the effect of adding a normal lubrication force is marginal.

[1] J. S. Bridge, *Rivers and Floodplains: Forms, Processes, and Sedimentary Record* (John Wiley & Sons, New York, NY, 2009).

[2] R. Soulsby, *Dynamics of Marine Sands* (Thomas Telford, London, UK, 1997).

[3] B. Andreotti, Y. Forterre, and O. Pouliquen, *Les milieux granulaires: Entre fluide et solide* (EDP Sciences, Les Ulis, France, 2012).

[4] É. Guazzelli and O. Pouliquen, Rheology of dense granular suspensions, *J. Fluid Mech.* **852**, P1 (2018).

[5] R. Jackson, *The Dynamics of Fluidized Particles* (Cambridge University Press, Cambridge, UK, 2000).

[6] R. Jackson, Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid, *Chem. Eng. Sci.* **52**, 2457 (1997).

[7] Y. Kim, Z. Cheng, T.-J. Hsu, and J. Chauchat, A numerical study of sheet flow under monochromatic nonbreaking waves using a free surface resolving Eulerian two-phase flow model, *J. Geophys. Res. [Oceans]* **123**, 4693 (2018).

[8] A. Einstein, Eine neue bestimmung der moleküldimensionen (in german), *Annalen der Physik* **324**, 289 (1906).

[9] A. Einstein, Correction of the paper 'eine neue bestimmung der moleküldimensionen (1906) (in german), *Annalen der Physik* **339**, 591 (1911).

[10] G. Batchelor and J. Green, The determination of the bulk stress in a suspension of spherical particles to order $c^2$, *J. Fluid Mech.* **56**, 401 (1972).

[11] L. Gibilaro, K. Gallucci, R. Di Felice, and P. Pagliai, On the apparent viscosity of a fluidized bed, *Chem. Eng. Sci.* **62**, 294 (2007).

[12] F. Boyer, E. Guazzelli, and O. Pouliquen, Unifying suspension and granular rheology, *Phys. Rev. Lett.* **107**, 188301 (2011).

[13] M. Trulsson, B. Andreotti, and P. Claudin, Transition from the viscous to inertial regime in dense suspensions, *Phys. Rev. Lett.* **109**, 118305 (2012).

[14] G. Ovarlez, F. Bertrand, and S. Rodts, Local determination of the constitutive law of a dense suspension of noncolloidal particles through magnetic resonance imaging, *J. Rheol.* **50**, 259 (2006).

[15] P. Aussillous, J. Chauchat, M. Pailha, M. Médale, and E. Guazzelli, Investigation of the mobile granular layer in bedload transport by laminar shearing flows, *J. Fluid Mech.* **736**, 594 (2013).

[16] M. Houssais, C. P. Ortiz, D. J. Durian, and D. J. Jerolmack, Rheology of sediment transported by a laminar flow, *Phys. Rev. E* **94**, 062609 (2016).

[17] M. Houssais and D. J. Jerolmack, Toward a unifying constitutive relation for sediment transport across environments, *Geomorphology* **277**, 251 (2017).

[18] R. Maurin, J. Chauchat, B. Chareyre, and P. Frey, A minimal coupled fluid-discrete element model for bedload transport, *Phys. Fluids* **27**, 113302 (2015).

[19] R. Maurin, J. Chauchat, and P. Frey, Dense granular flow rheology in turbulent bedload transport, *J. Fluid Mech.* **804**, 490 (2016).

[20] T. Pähtz and O. Durán, Universal friction law at granular solid-gas transition explains scaling of sediment transport load with excess fluid shear stress, *Phys. Rev. Fluids* **3**, 104302 (2018).
[21] T. Pähtz and O. Durán, The cessation threshold of nonsuspended sediment transport across aeolian and fluvial environments, *J. Geophys. Res.: Earth Surf.* **123**, 1638 (2018).

[22] T. Pähtz, O. Durán, D. N. de Klerk, I. Govender, and M. Trulsson, Local rheology relation with variable yield stress ratio across dry, wet, dense, and dilute granular flows, *Phys. Rev. Lett.* **123**, 048001 (2019).

[23] A. G. Kidanemariam and M. Uhlmann, Direct numerical simulation of pattern formation in subaqueous sediment, *J. Fluid Mech.* **750**, R2 (2014).

[24] A. G. Kidanemariam and M. Uhlmann, Formation of sediment patterns in channel flow: Minimal unstable systems and their temporal evolution, *J. Fluid Mech.* **818**, 716 (2017).

[25] M. Scherer, A. G. Kidanemariam, and M. Uhlmann, On the scaling of the instability of a flat sediment bed with respect to ripple-like patterns, *J. Fluid Mech.* **900**, A1 (2020).

[26] B. Vowinckel, E. Biegert, E. Meiburg, P. Aussillous, and É. Guazzelli, Rheology of mobile sediment beds sheared by viscous, pressure-driven flows, *J. Fluid Mech.* **921**, A20 (2021).

[27] Q. Zhang, E. Deal, J. T. Perron, J. G. Venditti, S. J. Benavides, M. Rushlow, and K. Kamrin, Fluid-driven transport of round sediment particles: From discrete simulations to continuum modeling, *J. Geophys. Res.: Earth Surf.* **127**, e2021JF006504 (2022).

[28] B. Bigot, T. Bonometti, L. Lacaze, and O. Thual, A simple immersed-boundary method for solid-fluid interaction in constant- and stratified-density flows, *Comput. Fluids* **97**, 126 (2014).

[29] J.-L. Pierson and J. Magnaudet, Inertial settling of a sphere through an interface. Part 2. sphere and tail dynamics, *J. Fluid Mech.* **835**, 808 (2018).

[30] E. Izard, Modélisation numérique des écoulements granulaires denses immergés dans un fluide, Ph.D. thesis, Dynamique des fluides, Institut National Polytechnique de Toulouse - INPT, Toulouse INP, 2014, http://ethesis.inp-toulouse.fr/archive/00002887/.

[31] P. Cundall and O. Strack, A discrete numerical model for granular assemblies, *Géotechnique* **29**, 47 (1979).

[32] O. Durán, B. Andreotti, and P. Claudin, Numerical simulation of turbulent sediment transport, from bed load to saltation, *Phys. Fluids* **24**, 103306 (2012).

[33] J.-L. Pierson, Traversée d'une interface entre deux fluides par une sphère, Ph.D. thesis, Dynamique des fluides, Institut National Polytechnique de Toulouse - INPT, 2015, https://oatao.univ-toulouse.fr/15754/.

[34] Y. Yuki, S. Takeuchi, and T. Kajishima, Efficient immersed boundary method for strong interaction problem of arbitrary shape object with the self-induced flow, *J. Fluid Sci. Technol.* **2**, 1 (2007).

[35] I. Calmet, Analyse par simulation des grandes échelles des mouvements turbulents et du transfert de masse sous une interface plane, Ph.D. thesis, Dynamique des fluides, Institut National Polytechnique de Toulouse - INPT, Toulouse INP, 1995, https://www.theses.fr/1995INPT110H.

[36] T. Ikeno and T. Kajishima, Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations, *J. Comput. Phys.* **226**, 1485 (2007).

[37] J. Shäfer, S. Dippel, and D. Wolf, Force schemes in simulations of granular materials, *J. Phys. I (France)* **6**, 5 (1996).

[38] X. Fang, J. Tang, and H. Luo, Granular damping analysis using an improved discrete element approach, *J. Sound Vib.* **308**, 112 (2007).

[39] F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Rheophysics of dense granular materials: Discrete simulation of plane shear flows, *Phys. Rev. E* **72**, 021309 (2005).

[40] O. Baran, D. Ertaş, T. C. Halsey, G. S. Grest, and J. B. Lechman, Velocity correlations in dense gravity-driven granular chute flow, *Phys. Rev. E* **74**, 051302 (2006).

[41] E. Izard, T. Bonometti, and L. Lacaze, Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid, *J. Fluid Mech.* **747**, 422 (2014).

[42] L. Lacaze, J. Bouteloup, B. Fry, and E. Izard, Immersed granular collapse: From viscous to free-fall unsteady granular flows, *J. Fluid Mech.* **912**, A15 (2021).

[43] C. W. Gear, *Numerical Initial Value Problems in Ordinary Differential Equations* (Prentice Hall PTR, Hoboken, NJ, 1971).

[44] M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids* (Oxford University Press, Oxford, UK, 1987).
[21] T. Pähtz and O. Durán, The cessation threshold of nonsuspended sediment transport across aeolian and fluvial environments, *J. Geophys. Res.: Earth Surf.* **123**, 1638 (2018).

[22] T. Pähtz, O. Durán, D. N. de Klerk, I. Govender, and M. Trulsson, Local rheology relation with variable yield stress ratio across dry, wet, dense, and dilute granular flows, *Phys. Rev. Lett.* **123**, 048001 (2019).

[23] A. G. Kidanemariam and M. Uhlmann, Direct numerical simulation of pattern formation in subaqueous sediment, *J. Fluid Mech.* **750**, R2 (2014).

[24] A. G. Kidanemariam and M. Uhlmann, Formation of sediment patterns in channel flow: Minimal unstable systems and their temporal evolution, *J. Fluid Mech.* **818**, 716 (2017).

[25] M. Scherer, A. G. Kidanemariam, and M. Uhlmann, On the scaling of the instability of a flat sediment bed with respect to ripple-like patterns, *J. Fluid Mech.* **900**, A1 (2020).

[26] B. Vowinckel, E. Biegert, E. Meiburg, P. Aussillous, and É. Guazzelli, Rheology of mobile sediment beds sheared by viscous, pressure-driven flows, *J. Fluid Mech.* **921**, A20 (2021).

[27] Q. Zhang, E. Deal, J. T. Perron, J. G. Venditti, S. J. Benavides, M. Rushlow, and K. Kamrin, Fluid-driven transport of round sediment particles: From discrete simulations to continuum modeling, *J. Geophys. Res.: Earth Surf.* **127**, e2021JF006504 (2022).

[28] B. Bigot, T. Bonometti, L. Lacaze, and O. Thual, A simple immersed-boundary method for solid-fluid interaction in constant- and stratified-density flows, *Comput. Fluids* **97**, 126 (2014).

[29] J.-L. Pierson and J. Magnaudet, Inertial settling of a sphere through an interface. Part 2. sphere and tail dynamics, *J. Fluid Mech.* **835**, 808 (2018).

[30] E. Izard, Modélisation numérique des écoulements granulaires denses immergés dans un fluide, Ph.D. thesis, Dynamique des fluides, Institut National Polytechnique de Toulouse - INPT, Toulouse INP, 2014, http://ethesis.inp-toulouse.fr/archive/00002887/.

[31] P. Cundall and O. Strack, A discrete numerical model for granular assemblies, *Géotechnique* **29**, 47 (1979).

[32] O. Durán, B. Andreotti, and P. Claudin, Numerical simulation of turbulent sediment transport, from bed load to saltation, *Phys. Fluids* **24**, 103306 (2012).

[33] J.-L. Pierson, Traversée d'une interface entre deux fluides par une sphère, Ph.D. thesis, Dynamique des fluides, Institut National Polytechnique de Toulouse - INPT, 2015, https://oatao.univ-toulouse.fr/15754/.

[34] Y. Yuki, S. Takeuchi, and T. Kajishima, Efficient immersed boundary method for strong interaction problem of arbitrary shape object with the self-induced flow, *J. Fluid Sci. Technol.* **2**, 1 (2007).

[35] I. Calmet, Analyse par simulation des grandes échelles des mouvements turbulents et du transfert de masse sous une interface plane, Ph.D. thesis, Dynamique des fluides, Institut National Polytechnique de Toulouse - INPT, Toulouse INP, 1995, https://www.theses.fr/1995INPT110H.

[36] T. Ikeno and T. Kajishima, Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations, *J. Comput. Phys.* **226**, 1485 (2007).

[37] J. Shäfer, S. Dippel, and D. Wolf, Force schemes in simulations of granular materials, *J. Phys. I (France)* **6**, 5 (1996).

[38] X. Fang, J. Tang, and H. Luo, Granular damping analysis using an improved discrete element approach, *J. Sound Vib.* **308**, 112 (2007).

[39] F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Rheophysics of dense granular materials: Discrete simulation of plane shear flows, *Phys. Rev. E* **72**, 021309 (2005).

[40] O. Baran, D. Ertaş, T. C. Halsey, G. S. Grest, and J. B. Lechman, Velocity correlations in dense gravity-driven granular chute flow, *Phys. Rev. E* **74**, 051302 (2006).

[41] E. Izard, T. Bonometti, and L. Lacaze, Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid, *J. Fluid Mech.* **747**, 422 (2014).

[42] L. Lacaze, J. Bouteloup, B. Fry, and E. Izard, Immersed granular collapse: From viscous to free-fall unsteady granular flows, *J. Fluid Mech.* **912**, A15 (2021).

[43] C. W. Gear, *Numerical Initial Value Problems in Ordinary Differential Equations* (Prentice Hall PTR, Hoboken, NJ, 1971).

[44] M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids* (Oxford University Press, Oxford, UK, 1987).
[45] T. Pöschel and T. Schwager, *Computational Granular Dynamics: Models and Algorithms* (Springer Science & Business Media, Cham, 2005).

[46] E. Izard, T. Bonometti, and L. Lacaze, Simulation of an avalanche in a fluid with a soft-sphere/immersed boundary method including a lubrication force, *J. Comput. Multiphase Flows* **6**, 391 (2014).

[47] S. B. Savage, M. H. Babaei, and T. Dabros, Modeling gravitational collapse of rectangular granular piles in air and water, *Mech. Res. Commun.* **56**, 1 (2014).

[48] C.-J. Lin, J. H. Peery, and W. R. Schowalter, Simple shear flow round a rigid sphere: Inertial effects and suspension rheology, *J. Fluid Mech.* **44**, 1 (1970).

[49] I. M. Krieger and T. J. Dougherty, A mechanism for non-Newtonian flow in suspensions of rigid spheres, *Trans. Soc. Rheol.* **3**, 137 (1959).

[50] P. Jop, Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, *Nature (London)* **441**, 727 (2006).

[51] GDR MiDi, On dense granular flows, Eur. Phys. J. E **14**, 341 (2004).

[52] C. Cassar, M. Nicolas, and O. Pouliquen, Submarine granular flows down inclined planes, *Phys. Fluids* **17**, 103301 (2005).

[53] H. Brinkman, A calculation of the viscosity and the sedimentation constant for solutions of large chain molecules taking into account the hampered flow of the solvent through these molecules, *Physica* **13**, 447 (1947).

[54] M. Ouriemi, P. Aussillous, and E. Guazzelli, Sediment dynamics. Part 1. Bed-load transport by laminar shearing flows, *J. Fluid Mech.* **636**, 295 (2009).

[55] F. Charru, J. Bouteloup, T. Bonometti, and L. Lacaze, Sediment transport and bedforms: A numerical study of two-phase viscous shear flow, *Meccanica* **51**, 3055 (2016).

[56] F. Tapia, M. Ichihara, O. Pouliquen, and É. Guazzelli, Viscous to inertial transition in dense granular suspension, *Phys. Rev. Lett.* **129**, 078001 (2022).

[57] A. S. Baumgarten and K. Kamrin, A general fluid-sediment mixture model and constitutive theory validated in many flow regimes, *J. Fluid Mech.* **861**, 721 (2019).

[58] J. Richardson and W. Zaki, The sedimentation of a suspension of uniform spheres under conditions of viscous flow, *Chem. Eng. Sci.* **3**, 65 (1954).

[59] H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, *Chem. Eng. Sci.* **16**, 242 (1961).
[45] T. Pöschel and T. Schwager, *Computational Granular Dynamics: Models and Algorithms* (Springer Science & Business Media, Cham, 2005).

[46] E. Izard, T. Bonometti, and L. Lacaze, Simulation of an avalanche in a fluid with a soft-sphere/immersed boundary method including a lubrication force, *J. Comput. Multiphase Flows* **6**, 391 (2014).

[47] S. B. Savage, M. H. Babaei, and T. Dabros, Modeling gravitational collapse of rectangular granular piles in air and water, *Mech. Res. Commun.* **56**, 1 (2014).

[48] C.-J. Lin, J. H. Peery, and W. R. Schowalter, Simple shear flow round a rigid sphere: Inertial effects and suspension rheology, *J. Fluid Mech.* **44**, 1 (1970).

[49] I. M. Krieger and T. J. Dougherty, A mechanism for non-Newtonian flow in suspensions of rigid spheres, *Trans. Soc. Rheol.* **3**, 137 (1959).

[50] P. Jop, Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, *Nature (London)* **441**, 727 (2006).

[51] GDR MiDi, On dense granular flows, Eur. Phys. J. E **14**, 341 (2004).

[52] C. Cassar, M. Nicolas, and O. Pouliquen, Submarine granular flows down inclined planes, *Phys. Fluids* **17**, 103301 (2005).

[53] H. Brinkman, A calculation of the viscosity and the sedimentation constant for solutions of large chain molecules taking into account the hampered flow of the solvent through these molecules, *Physica* **13**, 447 (1947).

[54] M. Ouriemi, P. Aussillous, and E. Guazzelli, Sediment dynamics. Part 1. Bed-load transport by laminar shearing flows, *J. Fluid Mech.* **636**, 295 (2009).

[55] F. Charru, J. Bouteloup, T. Bonometti, and L. Lacaze, Sediment transport and bedforms: A numerical study of two-phase viscous shear flow, *Meccanica* **51**, 3055 (2016).

[56] F. Tapia, M. Ichihara, O. Pouliquen, and É. Guazzelli, Viscous to inertial transition in dense granular suspension, *Phys. Rev. Lett.* **129**, 078001 (2022).

[57] A. S. Baumgarten and K. Kamrin, A general fluid-sediment mixture model and constitutive theory validated in many flow regimes, *J. Fluid Mech.* **861**, 721 (2019).

[58] J. Richardson and W. Zaki, The sedimentation of a suspension of uniform spheres under conditions of viscous flow, *Chem. Eng. Sci.* **3**, 65 (1954).

[59] H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, *Chem. Eng. Sci.* **16**, 242 (1961).