<figure><img src="image_1.png" /></figure>

HAL
open science

# From discrete to continuum description of weakly inertial bedload transport

Benjamin Fry, Laurent Lacaze, Thomas Bonometti, Pierre Elyakime, François Charru

► To cite this version:

Benjamin Fry, Laurent Lacaze, Thomas Bonometti, Pierre Elyakime, François Charru. From discrete to continuum description of weakly inertial bedload transport. Physical Review Fluids, 2024, 9 (2), pp.024304. 10.1103/physrevfluids.9.024304. hal-04532432

HAL Id: hal-04532432

https://ut3-toulouseinp.hal.science/hal-04532432

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# From discrete to continuum description of weakly inertial bedload transport

Benjamin Fry, Laurent Lacaze, Thomas Bonometti, Pierre Elyakime, François Charru

► To cite this version:

Benjamin Fry, Laurent Lacaze, Thomas Bonometti, Pierre Elyakime, François Charru. From discrete to continuum description of weakly inertial bedload transport. Physical Review Fluids, 2024, 9 (2), pp.024304. 10.1103/physrevfluids.9.024304. hal-04532432

HAL Id: hal-04532432

https://ut3-toulouseinp.hal.science/hal-04532432v1

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# From discrete to continuum description of weakly inertial bedload transport

Benjamin Fry, Laurent Lacaze,\*, Thomas Bonometti, Pierre Elyakime, and François Charru

Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse,
CNRS, Toulouse INP, UPS, Toulouse, France

<figure><img src="image_2.png" /></figure>

(Received 21 November 2019; accepted 8 January 2024; published 20 February 2024)

Granular bed motion induced by a liquid shear flow is considered as a model of bedload. This flow is characterized by a localization of the granular flow close to the upper surface of the bed. In such a framework, the paper aims at discussing and proposing continuum effective models describing the physical properties, as the rheology, of the two-phase granular suspension at a so-called mesoscale larger than the grain size. The main questions addressed here concern the relevance and the extension of effective rheological models for granular suspensions often extracted from simple-shear flow configurations in a viscous limit, to situations of both granular shear localization as in bedload and weakly inertial flow regime. For this purpose, we consider the steady transport of a granular bed by a laminar Couette fluid flow above it and close to the onset of motion, for a grain-to-fluid density ratio of 2.5, as silica in water, and a range of particle Reynolds numbers $Re_p \in [0.1, 10]$ and Shields numbers $\theta \in [0.1, 0.7]$. To provide accurate continuum models, the dynamics into the granular shear layer has to be first known down to a scale smaller than each grain, a so-called microscale. Numerical simulations are thus performed at the microscale at which individual grain dynamics is resolved, using an immersed boundary method (IBM) coupled to a discrete element method (DEM) granular solver. Upscaling is then performed to obtain the equivalent momentum balance at the mesoscale, characterized by continuum phases, using a spatial averaging method on volume element larger than the grain diameter. This approach allows us obtaining stresses, strains, and their relationships for the fluid phase, the granular phase, and the equivalent mixture, independently. The main contribution of this work is threefold: (i) we highlight the relevance of mesoscopic rheological continuum law for localized granular shear flow; (ii) we extract rheological models from direct numerical simulations (IBM/DEM) in a weakly inertial regime, going beyond purely viscous situations; and (iii) we extend Coulomb-like model $\mu(I)$ of a granular suspension to incorporate fluid/particle inertial effects showing a different dependence of fluid phase and granular phase contributions with dimensionless numbers.

DOI: 10.1103/PhysRevFluids.9.024304

## I. INTRODUCTION

Immersed granular transport occurs in many environmental situations and industrial processes, such as sediment transport in rivers and oceans [1,2], oil extraction, pharmaceutic or food processes [3], to mention but a few. The modeling of granular flows at these scales imposes a minimum scale of resolution which is greater than the grain size. For this reason, the granular phase is usually modeled as a continuum phase.

Yet, modeling immersed granular flows at a scale large enough so the granular medium can be viewed as a continuum (denoted mesoscale hereafter) remains a challenging task [4]. On the one

<footnote>\*laurent.lacaze@imft.fr</footnote>
# From discrete to continuum description of weakly inertial bedload transport

Benjamin Fry, Laurent Lacaze,\*, Thomas Bonometti, Pierre Elyakime, and François Charru

Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse,
CNRS, Toulouse INP, UPS, Toulouse, France

<figure><img src="image_2.png" /></figure>

(Received 21 November 2019; accepted 8 January 2024; published 20 February 2024)

Granular bed motion induced by a liquid shear flow is considered as a model of bedload. This flow is characterized by a localization of the granular flow close to the upper surface of the bed. In such a framework, the paper aims at discussing and proposing continuum effective models describing the physical properties, as the rheology, of the two-phase granular suspension at a so-called mesoscale larger than the grain size. The main questions addressed here concern the relevance and the extension of effective rheological models for granular suspensions often extracted from simple-shear flow configurations in a viscous limit, to situations of both granular shear localization as in bedload and weakly inertial flow regime. For this purpose, we consider the steady transport of a granular bed by a laminar Couette fluid flow above it and close to the onset of motion, for a grain-to-fluid density ratio of 2.5, as silica in water, and a range of particle Reynolds numbers $Re_p \in [0.1, 10]$ and Shields numbers $\theta \in [0.1, 0.7]$. To provide accurate continuum models, the dynamics into the granular shear layer has to be first known down to a scale smaller than each grain, a so-called microscale. Numerical simulations are thus performed at the microscale at which individual grain dynamics is resolved, using an immersed boundary method (IBM) coupled to a discrete element method (DEM) granular solver. Upscaling is then performed to obtain the equivalent momentum balance at the mesoscale, characterized by continuum phases, using a spatial averaging method on volume element larger than the grain diameter. This approach allows us obtaining stresses, strains, and their relationships for the fluid phase, the granular phase, and the equivalent mixture, independently. The main contribution of this work is threefold: (i) we highlight the relevance of mesoscopic rheological continuum law for localized granular shear flow; (ii) we extract rheological models from direct numerical simulations (IBM/DEM) in a weakly inertial regime, going beyond purely viscous situations; and (iii) we extend Coulomb-like model $\mu(I)$ of a granular suspension to incorporate fluid/particle inertial effects showing a different dependence of fluid phase and granular phase contributions with dimensionless numbers.

DOI: 10.1103/PhysRevFluids.9.024304

## I. INTRODUCTION

Immersed granular transport occurs in many environmental situations and industrial processes, such as sediment transport in rivers and oceans [1,2], oil extraction, pharmaceutical or food processes [3], to mention but a few. The modeling of granular flows at these scales imposes a minimum scale of resolution which is greater than the grain size. For this reason, the granular phase is usually modeled as a continuum phase.

Yet, modeling immersed granular flows at a scale large enough so the granular medium can be viewed as a continuum (denoted mesoscale hereafter) remains a challenging task [4]. On the one

<footnote>\*laurent.lacaze@imft.fr</footnote>
hand, two-phase Eulerian models such as those developed by R. Jackson [5,6] suffer from the lack of accurate and robust closure models [see, e.g., Ref. [7]]. In particular, the apparent rheological behavior of the various phases (fluid, granular, mixture) still raises some open questions. Moreover, the validity of this mesoscale description usually requires a clear separation between the grain size and the macroscopic scale of deformation of the granular material. This assumption can be drastic for localized shear and sharp interface as for instance when a bed of grains is sheared and transported by a fluid flow as encountered in sediment transport. The strategy to adopt to model the granular phase then remains uncertain.

In the viscous limit, the effective fluid phase at a scale larger than the grains is usually described as an equivalent fluid whose effective viscosity is a function of the particle volume fraction $\phi$. This allows to recover well-known results of the effective viscosity at low particle concentration [8–10]. The apparent viscosity of the fluid phase at larger particle concentration is approximated by semiempirical relations [11]. The influence of inertia on these models remains an open-question. As for the granular phase, several rheological models have been proposed. Recent advances on the understanding of dense granular flows have led to some rheological models for suspensions [12,13]. Finally, the rheology of the mixture is generally modeled by the Krieger and Dougherty (1959)'s correlation [14]. More recently, a mixture model has been proposed following recent developments on the rheology of dense granular medium [12]. All these models need to be assessed for a wider range of granular flows, including bedload sediment transport, for which the mesoscale description can suffer grain size effects as mentioned previously.

Yet, recent experimental works have shown the pertinence of large-scale mixture rheological model even for this drastic configuration [15–17]. At the same time, Euler-Lagrange approaches have been used to assess these rheological models for the granular phase [18,19] and the mixture phase [20]. These approaches are useful to cover a wide range of granular transport regimes (bedload, saltation, suspension). For instance, Pähtz and coworkers used such an approach to propose a model predicting the bulk transport cessation threshold for both aeolian and fluvial systems [21]. They were also able to assess the local rheological models in a large range of granular flows including wet, dry, dense, and dilute granular flows [22]. Note, however, that these Euler-Lagrange approaches make use of some approximations as the flow is solved at a larger scale than the particle scale. For this purpose, numerical simulations resolved at the grain scales are needed to complement these approaches. Bedload transport has been investigated by Uhlmann and coworkers using fully resolved simulations [23–25], however, to our knowledge, the question of the relevance of the mesoscale models was not addressed. Recently, Vowinckel *et al.* (2021) [26] used such an approach to investigate the rheology of mobile sediment beds in viscous pressure driven flows, for which the layer of carried grains remain thick. The latter study was focused on a specific flow regime associated with small particle Reynolds number for comparison with experiments presented in Ref. [15]. On the other hand, a similar approach, but based on lattice Boltzmann method for the fluid phase resolution, has been used to study the case of sediment transport in an inertial regime [27].

Here, we consider the generic problem of bedload sediment transport of a bed of grains, induced by a linear shear-driven fluid flow. The problem is sketched in Fig. 1(a). A bed of grains of mean diameter $d$ and density $\rho_p$ is subjected to a Couette flow with shear rate $\dot{\Gamma}$, in a fluid of density $\rho_f$ and viscosity $\eta_f$. For such typical flow, the subsequent granular transport is mostly confined to a relatively thin layer at the top of the bed. The present problem depends on three dimensionless parameters, namely the density ratio $r = \rho_p/\rho_f$, the particle Reynolds number $Re_p = \rho_f \dot{\Gamma}^2 d^2 / \eta_f$ and the Shields number $\theta = \eta_f \dot{\Gamma} / (\Delta\rho g d)$, where $\Delta\rho = \rho_p - \rho_f$ and $g$ is the gravitational acceleration. These numbers characterize the local dynamics at the grain scale, they are thus based on the particle diameter $d$ as the characteristic length scale and $\dot{\Gamma}$ as the inverse time scale. In the following, $\dot{\Gamma}$ is defined far from the granular bed where the fluid shear rate is constant. The two independent parameters used here, namely $Re_p$ and $\theta$, both depend on the shear rate. Alternatively, one may replace $Re_p$ or $\theta$ by the Galileo number $Ga = \sqrt{gd^3 \rho_f \Delta\rho / \eta_f} = \sqrt{Re_p / \theta}$ which only depends on
hand, two-phase Eulerian models such as those developed by R. Jackson [5,6] suffer from the lack of accurate and robust closure models [see, e.g., Ref. [7]]. In particular, the apparent rheological behavior of the various phases (fluid, granular, mixture) still raises some open questions. Moreover, the validity of this mesoscale description usually requires a clear separation between the grain size and the macroscopic scale of deformation of the granular material. This assumption can be drastic for localized shear and sharp interface as for instance when a bed of grains is sheared and transported by a fluid flow as encountered in sediment transport. The strategy to adopt to model the granular phase then remains uncertain.

In the viscous limit, the effective fluid phase at a scale larger than the grains is usually described as an equivalent fluid whose effective viscosity is a function of the particle volume fraction $\phi$. This allows to recover well-known results of the effective viscosity at low particle concentration [8–10]. The apparent viscosity of the fluid phase at larger particle concentration is approximated by semiempirical relations [11]. The influence of inertia on these models remains an open-question. As for the granular phase, several rheological models have been proposed. Recent advances on the understanding of dense granular flows have led to some rheological models for suspensions [12,13]. Finally, the rheology of the mixture is generally modeled by the Krieger and Dougherty (1959)'s correlation [14]. More recently, a mixture model has been proposed following recent developments on the rheology of dense granular medium [12]. All these models need to be assessed for a wider range of granular flows, including bedload sediment transport, for which the mesoscale description can suffer grain size effects as mentioned previously.

Yet, recent experimental works have shown the pertinence of large-scale mixture rheological model even for this drastic configuration [15–17]. At the same time, Euler-Lagrange approaches have been used to assess these rheological models for the granular phase [18,19] and the mixture phase [20]. These approaches are useful to cover a wide range of granular transport regimes (bedload, saltation, suspension). For instance, Pähtz and coworkers used such an approach to propose a model predicting the bulk transport cessation threshold for both aeolian and fluvial systems [21]. They were also able to assess the local rheological models in a large range of granular flows including wet, dry, dense, and dilute granular flows [22]. Note, however, that these Euler-Lagrange approaches make use of some approximations as the flow is solved at a larger scale than the particle scale. For this purpose, numerical simulations resolved at the grain scales are needed to complement these approaches. Bedload transport has been investigated by Uhlmann and coworkers using fully resolved simulations [23–25], however, to our knowledge, the question of the relevance of the mesoscale models was not addressed. Recently, Vowinckel *et al.* (2021) [26] used such an approach to investigate the rheology of mobile sediment beds in viscous pressure driven flows, for which the layer of carried grains remain thick. The latter study was focused on a specific flow regime associated with small particle Reynolds number for comparison with experiments presented in Ref. [15]. On the other hand, a similar approach, but based on lattice Boltzmann method for the fluid phase resolution, has been used to study the case of sediment transport in an inertial regime [27].

Here, we consider the generic problem of bedload sediment transport of a bed of grains, induced by a linear shear-driven fluid flow. The problem is sketched in Fig. 1(a). A bed of grains of mean diameter $d$ and density $\rho_p$ is subjected to a Couette flow with shear rate $\dot{\Gamma}$, in a fluid of density $\rho_f$ and viscosity $\eta_f$. For such typical flow, the subsequent granular transport is mostly confined to a relatively thin layer at the top of the bed. The present problem depends on three dimensionless parameters, namely the density ratio $r = \rho_p/\rho_f$, the particle Reynolds number $Re_p = \rho_f \dot{\Gamma}^2 d^2 / \eta_f$ and the Shields number $\theta = \eta_f \dot{\Gamma} / (\Delta\rho g d)$, where $\Delta\rho = \rho_p - \rho_f$ and $g$ is the gravitational acceleration. These numbers characterize the local dynamics at the grain scale, they are thus based on the particle diameter $d$ as the characteristic length scale and $\dot{\Gamma}$ as the inverse time scale. In the following, $\dot{\Gamma}$ is defined far from the granular bed where the fluid shear rate is constant. The two independent parameters used here, namely $Re_p$ and $\theta$, both depend on the shear rate. Alternatively, one may replace $Re_p$ or $\theta$ by the Galileo number $Ga = \sqrt{gd^3 \rho_f \Delta\rho / \eta_f} = \sqrt{Re_p / \theta}$ which only depends on
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Scheme of the considered problem at the microscale. (b) Domain of simulation. (c) Scheme of the considered problem at the mesoscale. (d) Vertical profile of the streamwise velocity in the fluid phase (blue, $\langle u \rangle^f$), the granular phase (red, $\langle u \rangle^p$), and solid volume fraction (black, $\phi$) for $Re_p = 1$ and $\theta = 0.67$. Note that the velocity profile of the granular phase is represented with a solid line for $\phi \ge 0.025$ and with a dashed line otherwise. Here, we use a characteristic length of the weighting function $h_g = d$. $h_{\text{bed}}$ is the mean height of the granular bed defined as the barycentre of the local granular flow rate $h_{\text{bed}} = \int y \, dq / \int dq = \int y \, \phi \, u_p \, dy / \int \phi \, u_p \, dy$ as suggested by Durán et al. (2012) [32].</figcaption></figure>

the fluid and particle properties. To finish with, note that Stokes numbers can also be found in the literature, characterizing the ratio between particle inertia and fluid dissipation, and can be defined here as combinations of the independent dimensionless numbers defined previously. Depending on the considered inertial term of the particle, vertical one due to the weight of the particle or horizontal one due to the shear-driven fluid forcing, the Stokes number can take different form. The former is $\sqrt{r}Ga$, the latter being $rRe_p$, respectively.

In the present work, we perform three-dimensional “fully resolved” simulations of this problem using an immersed boundary method (IBM) [28,29] coupled to a granular solver [30] based on a discrete element method (DEM) [31]. We refer to “fully resolved” simulations because almost all the scales of the flow are resolved (except some part of the lubrication effects below the fluid mesh resolution). By temporally and spatially averaging these results, we can compute the flow properties at the mesoscale, which only depend on the vertical direction.

The motivations of the present work are twofold: (1) to compute all the terms of the Eulerian two-phase equations and propose some parametrization of the closure terms; (2) to assess the existing rheological models of the various phases in the generic case of bedload transport localised to a thin (few grains) shear layer.

The paper first presents the IBM-DEM numerical approach and setup used for computing the flow at the microscale. In Sec. III, we describe some theoretical background, namely the upscaling procedure giving the Eulerian equations at the mesoscale and rheological models for the various phases. In Sec. IV, spatial and temporal averaging is performed and the flow properties and rheology at the mesoscale are obtained for each phase. Finally, the influence of the particle Reynolds number and Shields number on the mesoscale behavior is discussed.

## II. GOVERNING EQUATIONS AND NUMERICAL SETUP

### A. Fluid solver

The fluid solver is the in-house code JADIM based on an immersed boundary method (IBM) for computing the fluid flow around moving particles on a fixed Cartesian grid [28–30,33]. Note that the present IBM method is of the body-force type similar to that of Ref. [34]. Assuming
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Scheme of the considered problem at the microscale. (b) Domain of simulation. (c) Scheme of the considered problem at the mesoscale. (d) Vertical profile of the streamwise velocity in the fluid phase (blue, $\langle u \rangle^f$), the granular phase (red, $\langle u \rangle^p$), and solid volume fraction (black, $\phi$) for $Re_p = 1$ and $\theta = 0.67$. Note that the velocity profile of the granular phase is represented with a solid line for $\phi \ge 0.025$ and with a dashed line otherwise. Here, we use a characteristic length of the weighting function $h_g = d$. $h_{\text{bed}}$ is the mean height of the granular bed defined as the barycentre of the local granular flow rate $h_{\text{bed}} = \int y \, dq / \int dq = \int y \, \phi \, u_p \, dy / \int \phi \, u_p \, dy$ as suggested by Durán et al. (2012) [32].</figcaption></figure>

the fluid and particle properties. To finish with, note that Stokes numbers can also be found in the literature, characterizing the ratio between particle inertia and fluid dissipation, and can be defined here as combinations of the independent dimensionless numbers defined previously. Depending on the considered inertial term of the particle, vertical one due to the weight of the particle or horizontal one due to the shear-driven fluid forcing, the Stokes number can take different form. The former is $\sqrt{r}Ga$, the latter being $rRe_p$, respectively.

In the present work, we perform three-dimensional “fully resolved” simulations of this problem using an immersed boundary method (IBM) [28,29] coupled to a granular solver [30] based on a discrete element method (DEM) [31]. We refer to “fully resolved” simulations because almost all the scales of the flow are resolved (except some part of the lubrication effects below the fluid mesh resolution). By temporally and spatially averaging these results, we can compute the flow properties at the mesoscale, which only depend on the vertical direction.

The motivations of the present work are twofold: (1) to compute all the terms of the Eulerian two-phase equations and propose some parametrization of the closure terms; (2) to assess the existing rheological models of the various phases in the generic case of bedload transport localised to a thin (few grains) shear layer.

The paper first presents the IBM-DEM numerical approach and setup used for computing the flow at the microscale. In Sec. III, we describe some theoretical background, namely the upscaling procedure giving the Eulerian equations at the mesoscale and rheological models for the various phases. In Sec. IV, spatial and temporal averaging is performed and the flow properties and rheology at the mesoscale are obtained for each phase. Finally, the influence of the particle Reynolds number and Shields number on the mesoscale behavior is discussed.

## II. GOVERNING EQUATIONS AND NUMERICAL SETUP

### A. Fluid solver

The fluid solver is the in-house code JADIM based on an immersed boundary method (IBM) for computing the fluid flow around moving particles on a fixed Cartesian grid [28–30,33]. Note that the present IBM method is of the body-force type similar to that of Ref. [34]. Assuming
an incompressible Newtonian fluid of density $\rho_f$ and dynamic viscosity $\eta_f$, the flow evolution is described by

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{u}) = \frac{1}{\rho_f} \nabla \cdot \boldsymbol{\sigma} + \mathbf{f}_{\text{IBM}} + \mathbf{g}, \qquad (2) $$

where $\boldsymbol{\sigma} = -P\mathbf{I} + \eta_f(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ is the local hydrodynamic stress tensor. In the previous expressions, $\mathbf{u}$ and $P$ denote the local velocity and pressure, respectively, $\mathbf{g}$ is the gravitational acceleration and $\mathbf{f}_{\text{IBM}}$ is a volume force term used to account for the presence of solid particles in the domain, as will be explained in the following.

The present approach is based on a finite-volume method with a three steps mixed third-order Runge-Kutta/Crank-Nicolson scheme. Pressure is solved using a projection method [35]. Domain decomposition and message passing interface parallelization are performed.

Equations (1) and (2) are resolved in the whole domain, in fluid cells but also in cells containing the solid particles. As the boundary location of the moving spheres almost never coincide with the mesh, so-called local solid volume fraction $\alpha_{\text{IBM}}$ is introduced to materialize the particles. $\alpha_{\text{IBM}}$ is equal to one in cells filled with solid, zero in cells filled with fluid, and in the boundary region $\alpha_{\text{IBM}}$ is varying with a sinusoidal shape from one to zero. The characteristic length $\Delta_{\alpha}$ of the transition is typically of the order of one to three mesh cells of size $\Delta x$. In the present work, $\Delta_{\alpha} = 1.3\Delta x$ [28,34].

The body force term $\mathbf{f}_{\text{IBM}}$ is then defined as

$$ \mathbf{f}_{\text{IBM}} = \alpha_{\text{IBM}} \frac{\mathbf{U}_d - \tilde{\mathbf{u}}}{\Delta t}, \qquad (3) $$

where $\Delta t$ is the numerical time step, $\tilde{\mathbf{u}}$ is a predictor fluid velocity, and $\mathbf{U}_d$ is a prescribed particle velocity. $\mathbf{f}_{\text{IBM}}$ basically ensures the velocity to be the actual particle velocity within solid particles. Ikeno and Kajishima (2007) [36] showed that using the local particle velocity as the prescribed velocity, namely $\mathbf{U}_d = \mathbf{u}^p + \omega_p \times \mathbf{r}$ may tend to overestimate boundary layer effects (hence the drag on the body) in the fluid region near the particle surface. To reduce this effect, we impose $\mathbf{U}_d = \alpha_{\text{IBM}}(\mathbf{u}^p + \omega_p \times \mathbf{r}) + (1 - \alpha_{\text{IBM}})\mathbf{U}_I$, as done by these authors and Pierson (2015) [33], where $\mathbf{U}_I$ is a local fluid velocity resulting from a multidirectional linear interpolation of $\tilde{\mathbf{u}}$ over the neighboring cells.

A subloop on the entire forcing/Crank-Nicolson step is performed to improve the enforcement of the no-slip condition at the particle boundary [see Refs. [29,33], for more details]. More specifically, $\tilde{\mathbf{u}}$ is first computed without considering the particles, then the IBM force is computed, thirdly an unsteady Stokes equation is solved and finally $\tilde{\mathbf{u}}$ is updated. This forcing subloop is computationally expensive and is therefore used only twice per Runge-Kutta loop (tests showed that using more subloops did not improve the results significantly).

## B. Granular solver

The granular solver is the in-house code GraDyM [30] based on a soft-sphere discrete element method (DEM) which allows to describe multicontact grains interaction [31,37,38]. It solves the Newton's equations for linear and angular momentum, with velocity $\mathbf{u}_p$ and angular velocity $\omega_p$, which read

$$ \rho_p V_p \frac{d\mathbf{u}_p}{dt} = \rho_p V_p \mathbf{g} + \mathbf{f}_h + \sum_{q \neq p} \mathbf{f}_c^{pq}, \qquad (4) $$

$$ \mathbf{I}_p \frac{d\omega_p}{dt} = \mathcal{T}_h + \sum_{q \neq p} \mathcal{T}_c^{pq}, \qquad (5) $$
an incompressible Newtonian fluid of density $\rho_f$ and dynamic viscosity $\eta_f$, the flow evolution is described by

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{u}) = \frac{1}{\rho_f} \nabla \cdot \boldsymbol{\sigma} + \mathbf{f}_{\text{IBM}} + \mathbf{g}, \qquad (2) $$

where $\boldsymbol{\sigma} = -P\mathbf{I} + \eta_f(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ is the local hydrodynamic stress tensor. In the previous expressions, $\mathbf{u}$ and $P$ denote the local velocity and pressure, respectively, $\mathbf{g}$ is the gravitational acceleration and $\mathbf{f}_{\text{IBM}}$ is a volume force term used to account for the presence of solid particles in the domain, as will be explained in the following.

The present approach is based on a finite-volume method with a three steps mixed third-order Runge-Kutta/Crank-Nicolson scheme. Pressure is solved using a projection method [35]. Domain decomposition and message passing interface parallelization are performed.

Equations (1) and (2) are resolved in the whole domain, in fluid cells but also in cells containing the solid particles. As the boundary location of the moving spheres almost never coincide with the mesh, so-called local solid volume fraction $\alpha_{\text{IBM}}$ is introduced to materialize the particles. $\alpha_{\text{IBM}}$ is equal to one in cells filled with solid, zero in cells filled with fluid, and in the boundary region $\alpha_{\text{IBM}}$ is varying with a sinusoidal shape from one to zero. The characteristic length $\Delta_{\alpha}$ of the transition is typically of the order of one to three mesh cells of size $\Delta x$. In the present work, $\Delta_{\alpha} = 1.3\Delta x$ [28,34].

The body force term $\mathbf{f}_{\text{IBM}}$ is then defined as

$$ \mathbf{f}_{\text{IBM}} = \alpha_{\text{IBM}} \frac{\mathbf{U}_d - \tilde{\mathbf{u}}}{\Delta t}, \qquad (3) $$

where $\Delta t$ is the numerical time step, $\tilde{\mathbf{u}}$ is a predictor fluid velocity, and $\mathbf{U}_d$ is a prescribed particle velocity. $\mathbf{f}_{\text{IBM}}$ basically ensures the velocity to be the actual particle velocity within solid particles. Ikeno and Kajishima (2007) [36] showed that using the local particle velocity as the prescribed velocity, namely $\mathbf{U}_d = \mathbf{u}^p + \omega_p \times \mathbf{r}$ may tend to overestimate boundary layer effects (hence the drag on the body) in the fluid region near the particle surface. To reduce this effect, we impose $\mathbf{U}_d = \alpha_{\text{IBM}}(\mathbf{u}^p + \omega_p \times \mathbf{r}) + (1 - \alpha_{\text{IBM}})\mathbf{U}_I$, as done by these authors and Pierson (2015) [33], where $\mathbf{U}_I$ is a local fluid velocity resulting from a multidirectional linear interpolation of $\tilde{\mathbf{u}}$ over the neighboring cells.

A subloop on the entire forcing/Crank-Nicolson step is performed to improve the enforcement of the no-slip condition at the particle boundary [see Refs. [29,33], for more details]. More specifically, $\tilde{\mathbf{u}}$ is first computed without considering the particles, then the IBM force is computed, thirdly an unsteady Stokes equation is solved and finally $\tilde{\mathbf{u}}$ is updated. This forcing subloop is computationally expensive and is therefore used only twice per Runge-Kutta loop (tests showed that using more subloops did not improve the results significantly).

## B. Granular solver

The granular solver is the in-house code GraDyM [30] based on a soft-sphere discrete element method (DEM) which allows to describe multicontact grains interaction [31,37,38]. It solves the Newton's equations for linear and angular momentum, with velocity $\mathbf{u}_p$ and angular velocity $\omega_p$, which read

$$ \rho_p V_p \frac{d\mathbf{u}_p}{dt} = \rho_p V_p \mathbf{g} + \mathbf{f}_h + \sum_{q \neq p} \mathbf{f}_c^{pq}, \qquad (4) $$

$$ \mathbf{I}_p \frac{d\omega_p}{dt} = \mathcal{T}_h + \sum_{q \neq p} \mathcal{T}_c^{pq}, \qquad (5) $$
<table><caption>TABLE I. Parameters used in the simulations. Here, Re<sub>p</sub> = ρ<sub>f</sub> ḡd<sup>2</sup>/η<sub>f</sub> is particle Reynolds number, θ = η<sub>f</sub> ḡ/(Δρgd) is the Shields number, St = ρ<sub>p</sub> ḡd<sup>2</sup>/η<sub>f</sub> = rRe<sub>p</sub> is the Stokes number, and Ga = √gd<sup>3</sup>ρ<sub>f</sub>Δρ/η<sub>f</sub> = √Re<sub>p</sub>/θ is the Galileo number. In all the cases, the density ratio is r = ρ<sub>p</sub>/ρ<sub>f</sub> = 2.5.</caption><thead><tr><th></th><th>Re<sub>p</sub> = 0.1 (St = 0.25)</th><th>Re<sub>p</sub> = 1 (St = 2.5)</th><th>Re<sub>p</sub> = 10 (St = 25)</th></tr></thead><tbody><tr><td>θ = 0.11</td><td></td><td></td><td>+(Ga = 9.5)</td></tr><tr><td>θ = 0.22</td><td>× (Ga = 0.67)</td><td>× (Ga = 2.1)</td><td>× (Ga = 6.7)</td></tr><tr><td>θ = 0.33</td><td>*</td><td>*</td><td>*</td></tr><tr><td>θ = 0.4</td><td>○ (Ga = 0.50)</td><td>○ (Ga = 1.6)</td><td>○ (Ga = 5.0)</td></tr><tr><td>θ = 0.5</td><td>□ (Ga = 0.45)</td><td>□ (Ga = 1.4)</td><td>□ (Ga = 4.5)</td></tr><tr><td>θ = 0.67</td><td>◇ (Ga = 0.39)</td><td>◇ (Ga = 1.2)</td><td>◇ (Ga = 3.9)</td></tr></tbody></table>

where $V_p = \pi d^3/6$ is the volume of particle $p$, $I_p$ is the moment of inertia tensor, $\sum_{q \neq p} \mathbf{f}_c^{pq}$ ($\sum_{q \neq p} \mathcal{T}_c^{pq}$) is the total contact force (torque) induced by particles $q$ in contact with the particle $p$, and $\mathbf{f}_h$ ($\mathcal{T}_h$) is the total hydrodynamic force (torque) acting on the particle.

The contact force and torque $\mathbf{f}_c$ and $\mathcal{T}_c$ are computed using a soft-sphere model where the particles are allowed to slightly overlap, with a linear mass-spring model with dissipation in the normal direction and a mass-spring model coupled to a Coulomb type behavior above the slip threshold.

This model introduces four new dimensionless parameters, namely the dry restitution coefficient $e$, the characteristic contact time $\dot{\Gamma}t_{col}$ which are both related to dissipation during bouncing and contact stiffness; the ratio between the normal and tangential stiffnesses $k_t/k_n$ and the friction coefficient $\mu$ used for the Coulomb threshold. We choose a dry restitution coefficient of $e = 0.7$, a collision time small enough to ensure rigid behavior of the particles and loose dependency with this parameter [39,40], typically $t_{col} \le 10^{-3}\sqrt{d/g}$; here $k_t/k_n = 0.2857$ as suggested by Schafer et al. (1996) [37] and $\mu = 0.25$.

The hydrodynamic force and torque $\mathbf{f}_h$ and $\mathcal{T}_h$ acting on each particle are computed using $\mathbf{f}_{IBM}$ as

$$ \mathbf{f}_h = \int_{S_p} \sigma \cdot \mathbf{n} \, dS = \rho_f \left[ - \int_{V_p} \mathbf{f}_{IBM} \, d\nu + \frac{d}{dt} \left( \int_{V_p} \mathbf{u} \, d\nu \right) - V_p \mathbf{g} \right], \quad (6) $$

$$ \mathcal{T}_h = \rho_f \left[ - \int_{V_p} \mathbf{r} \times \mathbf{f}_{IBM} d\nu + \frac{d}{dt} \int_{V_p} \mathbf{r} \times \mathbf{u} d\nu \right], \quad (7) $$

where $S_p$ ($V_p$) is the surface (volume) of particle $p$, $\mathbf{f}_{IBM}$ is the body-force source term used in the fluid solver to impose the local particle velocity at the particle location. As explained in previous studies [41], an extra lubrication force close to contact between particles can be required to capture unavailable details of the fluid flow in the gap between two bouncing particles. However, as done in Ref. [42], no lubrication force is used here as it is shown not to affect significantly the present results, in particular those related to the rheological laws (see Appendix D).

Equations (4) and (5) are solved using the third-order Gear's predictor-corrector scheme [39,43–45]. The time step for the time advancement of the procedure is set to $\Delta t_{DEM} = t_{col}/50$.

The reader is referred to Refs. [28,29] for a validation of the present immersed boundary method without any contact and Refs. [41,46] for validations in the case where contacts are present.

### C. Numerical setup

In the present work, we set the particle-to-fluid density ratio to $r = 2.5$, we vary the Shields number in the range $0.1 \le \theta \le 0.7$ and the particle Reynolds number in the range $0.1 \le Re_p \le 10$ ($0.39 \le Ga \le 10$). Most of the explored values of $\theta$ are above the transport threshold and below the threshold of suspension initiation; hence the granular bed is moving as bedload. The various simulations are summarized in Table I. As the particle Reynolds numbers is lower than 10, this study
<table><caption>TABLE I. Parameters used in the simulations. Here, $Re_p = \rho_f \dot{\Gamma} d^2 / \eta_f$ is particle Reynolds number, $\theta = \eta_f \dot{\Gamma} / (\Delta \rho g d)$ is the Shields number, $St = \rho_p \dot{\Gamma} d^2 / \eta_f = r Re_p$ is the Stokes number, and $Ga = \sqrt{gd^3 \rho_f \Delta \rho} / \eta_f = \sqrt{Re_p / \theta}$ is the Galileo number. In all the cases, the density ratio is $r = \rho_p / \rho_f = 2.5$.</caption><thead><tr><th></th><th>Re<sub>p</sub> = 0.1 (St = 0.25)</th><th>Re<sub>p</sub> = 1 (St = 2.5)</th><th>Re<sub>p</sub> = 10 (St = 25)</th></tr></thead><tbody><tr><td>θ = 0.11</td><td></td><td></td><td>+(Ga = 9.5)</td></tr><tr><td>θ = 0.22</td><td>× (Ga = 0.67)</td><td>× (Ga = 2.1)</td><td>× (Ga = 6.7)</td></tr><tr><td>θ = 0.33</td><td>*</td><td>*</td><td>*</td></tr><tr><td>θ = 0.4</td><td>○ (Ga = 0.50)</td><td>○ (Ga = 1.6)</td><td>○ (Ga = 5.0)</td></tr><tr><td>θ = 0.5</td><td>□ (Ga = 0.45)</td><td>□ (Ga = 1.4)</td><td>□ (Ga = 4.5)</td></tr><tr><td>θ = 0.67</td><td>◇ (Ga = 0.39)</td><td>◇ (Ga = 1.2)</td><td>◇ (Ga = 3.9)</td></tr></tbody></table>

where $V_p = \pi d^3/6$ is the volume of particle $p$, $I_p$ is the moment of inertia tensor, $\sum_{q \neq p} \mathbf{f}_c^{pq}$ ($\sum_{q \neq p} \mathcal{T}_c^{pq}$) is the total contact force (torque) induced by particles $q$ in contact with the particle $p$, and $\mathbf{f}_h$ ($\mathcal{T}_h$) is the total hydrodynamic force (torque) acting on the particle.

The contact force and torque $\mathbf{f}_c$ and $\mathcal{T}_c$ are computed using a soft-sphere model where the particles are allowed to slightly overlap, with a linear mass-spring model with dissipation in the normal direction and a mass-spring model coupled to a Coulomb type behavior above the slip threshold.

This model introduces four new dimensionless parameters, namely the dry restitution coefficient $e$, the characteristic contact time $\dot{\Gamma} t_{col}$ which are both related to dissipation during bouncing and contact stiffness; the ratio between the normal and tangential stiffnesses $k_t/k_n$ and the friction coefficient $\mu$ used for the Coulomb threshold. We choose a dry restitution coefficient of $e = 0.7$, a collision time small enough to ensure rigid behavior of the particles and loose dependency with this parameter [39,40], typically $t_{col} \le 10^{-3} \sqrt{d/g}$; here $k_t/k_n = 0.2857$ as suggested by Schafer et al. (1996) [37] and $\mu = 0.25$.

The hydrodynamic force and torque $\mathbf{f}_h$ and $\mathcal{T}_h$ acting on each particle are computed using $\mathbf{f}_{IBM}$ as

$$ \mathbf{f}_h = \int_{S_p} \sigma \cdot \mathbf{n} \, dS = \rho_f \left[ - \int_{V_p} \mathbf{f}_{IBM} \, d\nu + \frac{d}{dt} \left( \int_{V_p} \mathbf{u} \, d\nu \right) - V_p \mathbf{g} \right], \quad (6) $$

$$ \mathcal{T}_h = \rho_f \left[ - \int_{V_p} \mathbf{r} \times \mathbf{f}_{IBM} d\nu + \frac{d}{dt} \int_{V_p} \mathbf{r} \times \mathbf{u} d\nu \right], \quad (7) $$

where $S_p$ ($V_p$) is the surface (volume) of particle $p$, $\mathbf{f}_{IBM}$ is the body-force source term used in the fluid solver to impose the local particle velocity at the particle location. As explained in previous studies [41], an extra lubrication force close to contact between particles can be required to capture unavailable details of the fluid flow in the gap between two bouncing particles. However, as done in Ref. [42], no lubrication force is used here as it is shown not to affect significantly the present results, in particular those related to the rheological laws (see Appendix D).

Equations (4) and (5) are solved using the third-order Gear's predictor-corrector scheme [39,43–45]. The time step for the time advancement of the procedure is set to $\Delta t_{DEM} = t_{col}/50$.

The reader is referred to Refs. [28,29] for a validation of the present immersed boundary method without any contact and Refs. [41,46] for validations in the case where contacts are present.

### C. Numerical setup

In the present work, we set the particle-to-fluid density ratio to $r = 2.5$, we vary the Shields number in the range $0.1 \le \theta \le 0.7$ and the particle Reynolds number in the range $0.1 \le Re_p \le 10$ ($0.39 \le Ga \le 10$). Most of the explored values of $\theta$ are above the transport threshold and below the threshold of suspension initiation; hence the granular bed is moving as bedload. The various simulations are summarized in Table I. As the particle Reynolds numbers is lower than 10, this study
<figure><img src="image_4.png" /><figcaption>FIG. 2. Time evolution of the instantaneous granular flow rate per unit depth scaled by $V_s$ with $V_s = (\rho_p - \rho_f)gd^2/(18\eta_f)$ the Stokes velocity ($\text{Re}_p = 1, \theta = 0.67, \text{Ga} = 1.22$). Black line: IBM-DEM results. Red line: mean value at steady state. Blue line: $q_{\text{conv}}$ as defined using Eq. (8) with $t_{\text{max}} = 1300\dot{\Gamma}^{-1}$).</figcaption></figure>

may be considered as relevant for all the situations where the bedload layer lies within the viscous sublayer of turbulent flows. Note that in Ref. [26], the simulations were for a fixed $\text{Ga} = 0.85$ and a range of Reynolds number $0.1 \lesssim \text{Re}_p \lesssim 1$ (using our scaling) and hence, our simulations covered a wider range of Galileo numbers and Reynolds numbers.

To vary the Shields number $\theta$, we varied the intensity of the gravitational acceleration $g$ in the range $[1, 6]d\dot{\Gamma}^2$, which does not modify the particle Reynolds number $\text{Re}_p$. We varied $\dot{\Gamma}$ and $\eta_f$ simultaneously to vary $\text{Re}_p$ without modifying the Shields parameter $\theta$. For this, we kept the product $\dot{\Gamma}\eta_f$ constant while modifying the ratio $\dot{\Gamma}/\eta_f$.

The flow configuration is sketched in Fig. 1(a). We define the Cartesian frame axes $(x, y, z)$ such as the $x$ and $z$ axes lie in the horizontal plane (perpendicular to gravity), $x$ being aligned with the imposed fluid flow velocity, and $y$ is the vertical direction. The simulations are performed in a $V_{\text{tot}} = L_x \times L_y \times L_z = 4d \times 20d \times 4d$ computational domain. In addition, we set the spatial resolution of the fluid solver to $\Delta x = \Delta y = \Delta z = 0.1d$. Periodic boundary conditions are imposed along the horizontal $x$ and $z$ directions. Note that we verified that varying the size of the domain and spatial resolution did not change the results significantly (see Appendix A for more details). The bottom of the domain, of size $L_x \times L_z = 4d \times 4d$, consists in a rough wall composed of a network of fixed half-grains. Moreover, we benefit from the steadiness of the flow to extract converged average quantities of the granular phase for such a small extension of the spatial domain. Half of the domain is initially filled with $N_p = 170$ grains. The particles have a mean diameter $d$ with a dispersion of 5% around the mean value to avoid crystallization and segregation phenomena. The speed at the top of the domain is imposed as $\dot{\Gamma}L_y/2$ (where $L_y/2$ is the fluid thickness), or equivalently $L_y/(2d) = 10$ in its dimensionless form. Then the bulk Reynolds number of the fluid phase $\text{Re}_B = \rho_f\dot{\Gamma}L_y^2/8\eta_f \equiv 50\text{Re}_p \in [5, 500]$, confirming the laminar state of the fluid phase flow.

The initial state of the simulation is obtained as follows. First we let the particles settle down without any interstitial fluid (i.e., using the GraDyM code only in the dry configuration). Second, we fix the particles (the granular bed thus becomes a porous medium) and simulate the Couette flow until a steady state is reached (i.e., using the JADIM code only). Finally, at $t = 0$, the grains are released and become free to move and rotate.

The physical duration of a run is then about $1300\dot{\Gamma}^{-1}$ including a transient regime of about $200\dot{\Gamma}^{-1}$ prior a quasi-steady state of the fluid-particle flow is reached. As an example, Fig. 2 presents the instantaneous dimensionless granular flow rate $q/V_s d$, with $V_s$ the Stokes velocity of a settling particle as defined in the caption of Fig. 2, for $\text{Re}_p = 1$ and $\theta = 0.67$ ($\text{Ga} = 1.22$). To ensure statistical convergence, we compute the cumulative time averaging of the considered quantity by
<figure><img src="image_4.png" /><figcaption>FIG. 2. Time evolution of the instantaneous granular flow rate per unit depth scaled by $V_s$ with $V_s = (\rho_p - \rho_f)gd^2/(18\eta_f)$ the Stokes velocity ($\text{Re}_p = 1, \theta = 0.67, \text{Ga} = 1.22$). Black line: IBM-DEM results. Red line: mean value at steady state. Blue line: $q_{\text{conv}}$ as defined using Eq. (8) with $t_{\text{max}} = 1300\dot{\Gamma}^{-1}$).</figcaption></figure>

may be considered as relevant for all the situations where the bedload layer lies within the viscous sublayer of turbulent flows. Note that in Ref. [26], the simulations were for a fixed $\text{Ga} = 0.85$ and a range of Reynolds number $0.1 \lesssim \text{Re}_p \lesssim 1$ (using our scaling) and hence, our simulations covered a wider range of Galileo numbers and Reynolds numbers.

To vary the Shields number $\theta$, we varied the intensity of the gravitational acceleration $g$ in the range $[1, 6]d\dot{\Gamma}^2$, which does not modify the particle Reynolds number $\text{Re}_p$. We varied $\dot{\Gamma}$ and $\eta_f$ simultaneously to vary $\text{Re}_p$ without modifying the Shields parameter $\theta$. For this, we kept the product $\dot{\Gamma}\eta_f$ constant while modifying the ratio $\dot{\Gamma}/\eta_f$.

The flow configuration is sketched in Fig. 1(a). We define the Cartesian frame axes $(x, y, z)$ such as the $x$ and $z$ axes lie in the horizontal plane (perpendicular to gravity), $x$ being aligned with the imposed fluid flow velocity, and $y$ is the vertical direction. The simulations are performed in a $V_{\text{tot}} = L_x \times L_y \times L_z = 4d \times 20d \times 4d$ computational domain. In addition, we set the spatial resolution of the fluid solver to $\Delta x = \Delta y = \Delta z = 0.1d$. Periodic boundary conditions are imposed along the horizontal $x$ and $z$ directions. Note that we verified that varying the size of the domain and spatial resolution did not change the results significantly (see Appendix A for more details). The bottom of the domain, of size $L_x \times L_z = 4d \times 4d$, consists in a rough wall composed of a network of fixed half-grains. Moreover, we benefit from the steadiness of the flow to extract converged average quantities of the granular phase for such a small extension of the spatial domain. Half of the domain is initially filled with $N_p = 170$ grains. The particles have a mean diameter $d$ with a dispersion of 5% around the mean value to avoid crystallization and segregation phenomena. The speed at the top of the domain is imposed as $\dot{\Gamma}L_y/2$ (where $L_y/2$ is the fluid thickness), or equivalently $L_y/(2d) = 10$ in its dimensionless form. Then the bulk Reynolds number of the fluid phase $\text{Re}_B = \rho_f\dot{\Gamma}L_y^2/8\eta_f \equiv 50\text{Re}_p \in [5, 500]$, confirming the laminar state of the fluid phase flow.

The initial state of the simulation is obtained as follows. First we let the particles settle down without any interstitial fluid (i.e., using the GraDyM code only in the dry configuration). Second, we fix the particles (the granular bed thus becomes a porous medium) and simulate the Couette flow until a steady state is reached (i.e., using the JADIM code only). Finally, at $t = 0$, the grains are released and become free to move and rotate.

The physical duration of a run is then about $1300\dot{\Gamma}^{-1}$ including a transient regime of about $200\dot{\Gamma}^{-1}$ prior a quasi-steady state of the fluid-particle flow is reached. As an example, Fig. 2 presents the instantaneous dimensionless granular flow rate $q/V_s d$, with $V_s$ the Stokes velocity of a settling particle as defined in the caption of Fig. 2, for $\text{Re}_p = 1$ and $\theta = 0.67$ ($\text{Ga} = 1.22$). To ensure statistical convergence, we compute the cumulative time averaging of the considered quantity by
starting from the end of the simulation (corresponding to a time denoted $t_{max}$). In particular, we define $f_{conv}$ for an arbitrary time-dependent quantity $f$, as

$$f_{conv}(t) = \frac{1}{t_{max} - t} \int_{t}^{t_{max}} f(\tau) d\tau. \qquad (8)$$

For instance, $q_{conv}/V_s d$ is shown in Fig. 2 (solid blue line). Starting from $t_{max}$, the function converge to a mean value as $t$ decreases. Note that the convergence should disappear when $t \to 0$ as the initial transient regime is then included in the integral of Eq. (8). Therefore, we can consider that the statistical convergence is reached when a plateau is observed on a significant interval of time $t$. In Fig. 2, the plateau of statical convergence of $q/V_s d$ is obtained between $200\bar{\Gamma}^{-1}$ and $600\bar{\Gamma}^{-1}$. In our case, we therefore perform temporal averaging over the time interval $[550\bar{\Gamma}^{-1}, t_{max}]$ using 150 samples separated by a time laps of $5\bar{\Gamma}t$. The obtained mean value for the granular flow rate $q/V_s d$ is shown in Fig. 2 (solid red line).

### III. MESOSCALE STRESSES: THEORETICAL BACKGROUND

#### A. Volume average approach

To obtain the mesoscopic mass and momentum equations from the microscopic (particle size) ones, Jackson proposed a spatial averaging procedure which is summarized below, for a steady flow, and keeping his notation. The reader is referred to Jackson (1997, 2000) [5,6] for more details. This spatial averaging allows to obtain an equivalent fluid phase and an equivalent solid phase which are both defined at any position $\mathbf{x}$ in the domain $V_{tot}$. These two phases are therefore weighted by mean of a local volume fraction for both the fluid phase $\epsilon(\mathbf{x})$ and the solid phase $\phi(\mathbf{x}) = 1 - \epsilon(\mathbf{x})$ (as both fluid and solid particle are considered incompressible). The latter quantities represent the portion of each phases on a length scale $l$ around $\mathbf{x}$ prescribed through a spatial averaging—or weighting—function $\mathcal{G}$, such as

$$\epsilon(\mathbf{x}) = \int_{V_f} \hat{\mathcal{G}}(|\mathbf{x} - \mathbf{x}'|) dv_{\mathbf{x}'}, \qquad (9)$$

$$\phi(\mathbf{x}) = \sum_{p} \int_{V_p} \hat{\mathcal{G}}(|\mathbf{x} - \mathbf{x}'|) dv_{\mathbf{x}'}, \qquad (10)$$

where $\mathbf{x} = (x, y, z)^T$, $\mathbf{x}' = (x', y', z')^T$, $V_p$ is the volume of the particle $p$, and $V_f$ the part of the volume $V_{tot}$ filled with fluid. To be consistent, $\mathcal{G}$ has to be defined as a decreasing function from the considered position $\mathbf{x}$ with a characteristic spatial extension $l$, normalized ($\int_{V_{tot}} \hat{\mathcal{G}} dv = 1$) and differentiable. Note that due to the symmetries of the considered system, averaged quantities shall only be a function of the vertical direction $y$, while remaining invariant of the horizontal directions $(x, z)$. The averaging length scale $l$ has thus only to be applied in the vertical direction, $l = h_g$. In other words, we choose a weighting function invariant in the horizontal directions, as the flow is, and Gaussian in the vertical direction, as

$$\hat{\mathcal{G}}(|\mathbf{x} - \mathbf{x}'|) \equiv \mathcal{G}(|y - y'|) = \frac{1}{L_x L_z} \frac{e^{-(\frac{y-y'}{h_g})^2}}{\int_0^{L_y} e^{-(\frac{y-y''}{h_g})^2} dy''}. \qquad (11)$$

Using Eqs. (11), (9), and (10) then read

$$\epsilon(y) = \int_{V_{tot}} \mathcal{G}(|y - y'|)[1 - \alpha_{IBM}(\mathbf{x}')] dv_{\mathbf{x}'}, \qquad (12)$$

$$\phi(y) = \sum_{p} \int_{V_{tot}} \mathcal{G}(|y - y'|) \alpha_{IBM}(\mathbf{x}') dv_{\mathbf{x}'}, \qquad (13)$$

recalling $\alpha_{IBM} = 0$ in the fluid phase and $\alpha_{IBM} = 1$ in the particles.
starting from the end of the simulation (corresponding to a time denoted $t_{max}$). In particular, we define $f_{conv}$ for an arbitrary time-dependent quantity $f$, as

$$f_{conv}(t) = \frac{1}{t_{max} - t} \int_{t}^{t_{max}} f(\tau) d\tau. \qquad (8)$$

For instance, $q_{conv}/V_s d$ is shown in Fig. 2 (solid blue line). Starting from $t_{max}$, the function converge to a mean value as $t$ decreases. Note that the convergence should disappear when $t \to 0$ as the initial transient regime is then included in the integral of Eq. (8). Therefore, we can consider that the statistical convergence is reached when a plateau is observed on a significant interval of time $t$. In Fig. 2, the plateau of statical convergence of $q/V_s d$ is obtained between $200\bar{\Gamma}^{-1}$ and $600\bar{\Gamma}^{-1}$. In our case, we therefore perform temporal averaging over the time interval $[550\bar{\Gamma}^{-1}, t_{max}]$ using 150 samples separated by a time laps of $5\bar{\Gamma}t$. The obtained mean value for the granular flow rate $q/V_s d$ is shown in Fig. 2 (solid red line).

### III. MESOSCALE STRESSES: THEORETICAL BACKGROUND

#### A. Volume average approach

To obtain the mesoscopic mass and momentum equations from the microscopic (particle size) ones, Jackson proposed a spatial averaging procedure which is summarized below, for a steady flow, and keeping his notation. The reader is referred to Jackson (1997, 2000) [5,6] for more details. This spatial averaging allows to obtain an equivalent fluid phase and an equivalent solid phase which are both defined at any position $\mathbf{x}$ in the domain $V_{tot}$. These two phases are therefore weighted by mean of a local volume fraction for both the fluid phase $\epsilon(\mathbf{x})$ and the solid phase $\phi(\mathbf{x}) = 1 - \epsilon(\mathbf{x})$ (as both fluid and solid particle are considered incompressible). The latter quantities represent the portion of each phases on a length scale $l$ around $\mathbf{x}$ prescribed through a spatial averaging—or weighting—function $\mathcal{G}$, such as

$$\epsilon(\mathbf{x}) = \int_{V_f} \hat{\mathcal{G}}(|\mathbf{x} - \mathbf{x}'|) dv_{\mathbf{x}'}, \qquad (9)$$

$$\phi(\mathbf{x}) = \sum_{p} \int_{V_p} \hat{\mathcal{G}}(|\mathbf{x} - \mathbf{x}'|) dv_{\mathbf{x}'}, \qquad (10)$$

where $\mathbf{x} = (x, y, z)^T$, $\mathbf{x}' = (x', y', z')^T$, $V_p$ is the volume of the particle $p$, and $V_f$ the part of the volume $V_{tot}$ filled with fluid. To be consistent, $\mathcal{G}$ has to be defined as a decreasing function from the considered position $\mathbf{x}$ with a characteristic spatial extension $l$, normalized ($\int_{V_{tot}} \hat{\mathcal{G}} dv = 1$) and differentiable. Note that due to the symmetries of the considered system, averaged quantities shall only be a function of the vertical direction $y$, while remaining invariant of the horizontal directions $(x, z)$. The averaging length scale $l$ has thus only to be applied in the vertical direction, $l = h_g$. In other words, we choose a weighting function invariant in the horizontal directions, as the flow is, and Gaussian in the vertical direction, as

$$\hat{\mathcal{G}}(|\mathbf{x} - \mathbf{x}'|) \equiv \mathcal{G}(|y - y'|) = \frac{1}{L_x L_z} \frac{e^{-(\frac{y-y'}{h_g})^2}}{\int_0^{L_y} e^{-(\frac{y-y''}{h_g})^2} dy''}. \qquad (11)$$

Using Eqs. (11), (9), and (10) then read

$$\epsilon(y) = \int_{V_{tot}} \mathcal{G}(|y - y'|)[1 - \alpha_{IBM}(\mathbf{x}')] dv_{\mathbf{x}'}, \qquad (12)$$

$$\phi(y) = \sum_{p} \int_{V_{tot}} \mathcal{G}(|y - y'|) \alpha_{IBM}(\mathbf{x}') dv_{\mathbf{x}'}, \qquad (13)$$

recalling $\alpha_{IBM} = 0$ in the fluid phase and $\alpha_{IBM} = 1$ in the particles.
An example of the vertical profile of the steady solid volume fraction $\phi(y)$ obtained in the present configuration of a sheared granular bed for $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22) is shown in Fig. 1(d) (black line). Here $h_g = d$ has been chosen as the characteristic length scale of the weight function (11). As expected, one observes that $\phi$ is close to its maximum packing fraction $\approx 0.6$ in a static region close to the bottom and rapidly decreases close to its upper surface to reach zero in top fluid layer. The vertical variation of $\phi$ depends on the dynamics of the moving layer which will be discussed later on in this section.

As the grains are rigid, we can also compute the number of particles per unit volume $n(y)$ at the altitude $y$ as

$$n(y) = \sum_{p} \mathcal{G}(|y - y^p|), \qquad (14)$$

where $y^p$ is the vertical coordinate of the position of particle $p$ center. Note that assuming $y - y^p \ll h_g$ over a particle diameter $d$, i.e., $\mathcal{G}$ varies slowly over a particle diameter and thus $d/2 \ll h_g$, one can easily link $\phi$ and $n$ as [5]

$$\phi(y) = V_p n(y) + \mathcal{O}\left[\left(\frac{d}{2h_g}\right)^2\right]. \qquad (15)$$

It is worth mentioning here that $d/2 \ll h_g$ would be the first condition of scale separation necessary for the validity of the upscaling approach. As a second condition, the characteristic spatial extension $h_g$ of the weighting function $\mathcal{G}$ must be much smaller than the macroscopic scale $L$, namely $h_g \ll L$. The latter condition can be seen as preventing unexpected smoothing of flow properties gradients at the mesoscale. These conditions are not straightforwardly fulfilled in the present configuration as $L$ is typically associated with the vertical gradient length scale, which could be typically of order $d$.

## B. Eulerian two-phase equations

Following the procedure described in Sec. III A, the fluid phase average $\langle f \rangle^f$ of an arbitrary quantity $f$ obtained from the IBM solver is then defined as

$$\langle f \rangle^f(y) = \frac{1}{\epsilon(y)} \int_{V_{\text{tot}}} [1 - \alpha_{\text{IBM}}(\mathbf{x}')] f(\mathbf{x}') \mathcal{G}(y - y') d\mathbf{v}'_{\mathbf{x}}. \qquad (16)$$

As explained in Ref. [5], one may define two types of averages for the granular phase, namely the *particle phase average* (which uses $n$) and the *solid phase average* (which uses $\phi$). For simplicity, we only present the particle phase average. In the following, we refer indifferently the “solid” or “particulate” quantities to as “granular.” Then, the particle phase average $\langle f \rangle^p$ of an arbitrary quantity $f^p$ of particle $p$ is defined by

$$\langle f \rangle^p(y) = \frac{1}{n(y)} \sum_p f^p \mathcal{G}(|y - y^p|). \qquad (17)$$

For example, Fig. 1(d) shows the vertical profiles of the streamwise steady fluid velocity $\langle u_x \rangle^f(y)$ (blue line) and the streamwise steady granular velocity $\langle u_{p|x} \rangle^p(y)$ (red line) obtained for $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22). Again, $h_g = d$ has been chosen as the characteristic length scale of the weight function (11). The fluid phase is shown to be everywhere moving faster than the granular, as it is expected for such a fluid-driven transport configuration. Moreover, velocity profiles rapidly decreases in the moving layer to reach zero simultaneously with the volume fraction reaching its maximum packing fraction in the static layer.

Applying the averaging procedure to the equations solved at the microscale, namely Eqs. (2) and (4), leads to the mesoscopic system of Eulerian two-phase equations. The complete set of equations are found in Ref. [5]. For sake of clarity, we simplify the momentum equation for the
An example of the vertical profile of the steady solid volume fraction $\phi(y)$ obtained in the present configuration of a sheared granular bed for $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22) is shown in Fig. 1(d) (black line). Here $h_g = d$ has been chosen as the characteristic length scale of the weight function (11). As expected, one observes that $\phi$ is close to its maximum packing fraction $\approx 0.6$ in a static region close to the bottom and rapidly decreases close to its upper surface to reach zero in top fluid layer. The vertical variation of $\phi$ depends on the dynamics of the moving layer which will be discussed later on in this section.

As the grains are rigid, we can also compute the number of particles per unit volume $n(y)$ at the altitude $y$ as

$$n(y) = \sum_{p} \mathcal{G}(|y - y^p|), \qquad (14)$$

where $y^p$ is the vertical coordinate of the position of particle $p$ center. Note that assuming $y - y^p \ll h_g$ over a particle diameter $d$, i.e., $\mathcal{G}$ varies slowly over a particle diameter and thus $d/2 \ll h_g$, one can easily link $\phi$ and $n$ as [5]

$$\phi(y) = V_p n(y) + \mathcal{O}\left[\left(\frac{d}{2h_g}\right)^2\right]. \qquad (15)$$

It is worth mentioning here that $d/2 \ll h_g$ would be the first condition of scale separation necessary for the validity of the upscaling approach. As a second condition, the characteristic spatial extension $h_g$ of the weighting function $\mathcal{G}$ must be much smaller than the macroscopic scale $L$, namely $h_g \ll L$. The latter condition can be seen as preventing unexpected smoothing of flow properties gradients at the mesoscale. These conditions are not straightforwardly fulfilled in the present configuration as $L$ is typically associated with the vertical gradient length scale, which could be typically of order $d$.

## B. Eulerian two-phase equations

Following the procedure described in Sec. III A, the fluid phase average $\langle f \rangle^f$ of an arbitrary quantity $f$ obtained from the IBM solver is then defined as

$$\langle f \rangle^f(y) = \frac{1}{\epsilon(y)} \int_{V_{\text{tot}}} [1 - \alpha_{\text{IBM}}(\mathbf{x}')] f(\mathbf{x}') \mathcal{G}(y - y') d\mathbf{v}'_{\mathbf{x}}. \qquad (16)$$

As explained in Ref. [5], one may define two types of averages for the granular phase, namely the *particle phase average* (which uses $n$) and the *solid phase average* (which uses $\phi$). For simplicity, we only present the particle phase average. In the following, we refer indifferently the “solid” or “particulate” quantities to as “granular.” Then, the particle phase average $\langle f \rangle^p$ of an arbitrary quantity $f^p$ of particle $p$ is defined by

$$\langle f \rangle^p(y) = \frac{1}{n(y)} \sum_{p} f^p \mathcal{G}(|y - y^p|). \qquad (17)$$

For example, Fig. 1(d) shows the vertical profiles of the streamwise steady fluid velocity $\langle u_x \rangle^f(y)$ (blue line) and the streamwise steady granular velocity $\langle u_{p|x} \rangle^p(y)$ (red line) obtained for $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22). Again, $h_g = d$ has been chosen as the characteristic length scale of the weight function (11). The fluid phase is shown to be everywhere moving faster than the granular, as it is expected for such a fluid-driven transport configuration. Moreover, velocity profiles rapidly decreases in the moving layer to reach zero simultaneously with the volume fraction reaching its maximum packing fraction in the static layer.

Applying the averaging procedure to the equations solved at the microscale, namely Eqs. (2) and (4), leads to the mesoscopic system of Eulerian two-phase equations. The complete set of equations are found in Ref. [5]. For sake of clarity, we simplify the momentum equation for the
fluid phase and the granular phase, obtained at the order $O((d/L)^2)$, for the specific case of the steady flow considered here. These read, respectively,

$$ \rho_f \in \langle \mathbf{u} \rangle^f \cdot \nabla \langle \mathbf{u} \rangle^f = \nabla \cdot \mathbf{\Sigma}^f - n \langle \mathbf{f}_h \rangle^p + \rho_f \in \mathbf{g}, \quad (18) $$

$$ \rho_p \phi \langle \mathbf{u}_p \rangle^p \cdot \nabla \langle \mathbf{u}_p \rangle^p = \nabla \cdot \mathbf{\Sigma}^p + n \langle \mathbf{f}_h \rangle^p + \rho_p \phi \mathbf{g}, \quad (19) $$

which are 2D and for which variables only depend on y. Hence $\nabla = (0, d/dy)$, $\langle \mathbf{u} \rangle^f = (\langle u_x \rangle^f, \langle u_y \rangle^f)$, $\langle \mathbf{u}_p \rangle^p = (\langle u_{p|x} \rangle^p, \langle u_{p|y} \rangle^p)$ and so on. The fluid-particle interaction force denoted $\langle \mathbf{f}_h \rangle^p$ reads

$$ n(y) \langle \mathbf{f}_h \rangle^p(y) = \sum_p \mathcal{G}(|y - y^p|) \int_{S_p} (\boldsymbol{\sigma} \cdot \mathbf{n}) \, dS = \sum_p \mathcal{G}(|y - y^p|) \mathbf{f}_h(\mathbf{x}^p), \quad (20) $$

which can be explicitly calculated from the microscale using Eq. (6). Note that this fluid-particle force is usually split as $n \langle \mathbf{f}_h \rangle^p = n \langle \mathbf{f}_1 \rangle^p + \phi \nabla \cdot \mathbf{\Sigma}^f$, with $\phi \nabla \cdot \mathbf{\Sigma}^f$ the so-called generalized buoyancy. Then $\langle \mathbf{f}_1 \rangle^p$ corresponds to the rest of local hydrodynamics interaction including drag among others.

$\mathbf{\Sigma}^f$ and $\mathbf{\Sigma}^p$ are the 2D apparent stress tensors of the equivalent fluid phase and granular phase, respectively. They can be split into different contributions, defined from microscale, as

$$ \mathbf{\Sigma}^f = \epsilon \langle \boldsymbol{\sigma} \rangle^f + n \langle \mathbf{s}^f \rangle^p - \rho_f \epsilon \langle \mathbf{u}' \otimes \mathbf{u}' \rangle^f, \quad (21) $$

$$ \mathbf{\Sigma}^p = n \langle \mathbf{s}^s \rangle^p - \rho_p \phi \langle \mathbf{u}'_p \otimes \mathbf{u}'_p \rangle^p, \quad (22) $$

where $\mathbf{u}' = (1 - \alpha_{\text{IBM}}(\mathbf{y}))\mathbf{u} - \langle \mathbf{u} \rangle^f$ (respectively, $\mathbf{u}'_p = \mathbf{u}_p - \langle \mathbf{u}_p \rangle^p$) is the spatially fluctuating part of the local fluid (respectively, granular) velocity and the tensors $\langle \mathbf{s}^f \rangle^p$ and $\langle \mathbf{s}^s \rangle^p$ read

$$ n \langle \mathbf{s}^f \rangle^p = \frac{d}{2} \sum_p \mathcal{G}(|y - y^p|) \int_{S_p} (\boldsymbol{\sigma} \cdot \mathbf{n}) \otimes \mathbf{n} \, dS = \frac{d}{2} \sum_p \mathcal{G}(|y - y^p|) \\ \times \int_{V_{\text{tot}}} \left[ \alpha_{\text{IBM}} \rho_f \left( \frac{d\mathbf{u}}{dt} - \mathbf{f}_{\text{IBM}} \right) \otimes \mathbf{r} \right] dv, \quad (23) $$

$$ n \langle \mathbf{s}^s \rangle^p = \frac{d}{2} \sum_p \mathcal{G}(|y - y^p|) \sum_{q \neq p} \mathbf{f}_c^{pq} \otimes \mathbf{n}^{pq}, \quad (24) $$

where $\mathbf{n}^{pq} = (\mathbf{x}^q - \mathbf{x}^p)/(|\mathbf{x}^q - \mathbf{x}^p|)$ is the unit outward vector normal to the surface of the particle $p$ at the contact point.

It shall be noted that the effective fluid stress tensor (21) contains a contribution $n \langle \mathbf{s}^f \rangle^p$ coming from the presence of the grains, leading to a differential stress contribution arising from the vertical variation of the local fluid stress at the particle scale, and probably leading to a $\phi$-dependent effective fluid viscosity as discussed in the following.

Moreover, mass conservation for the fluid phase and the granular phase shall impose, for the shear flow configuration considered here, $\epsilon \langle u_y \rangle^f = 0$ and $\phi \langle u_{p|y} \rangle^p = 0$, respectively. Then, the left-hand sides of Eqs. (18) and (19) have to be zero. We suggest here to keep these terms for now, to verify in the following the validity of averaged momentum equations obtained from the microscale.

### C. Eulerian equations for the mixture

In the previous section, upscaling has been performed in a way that both the fluid phase and the granular phase remain separate. One may also represent the system by a single equivalent phase that contains both the fluid and grains. This phase is here called the mixture phase and its velocity is defined as

$$ \bar{\rho} \langle \mathbf{u} \rangle^m = \rho_p \phi \langle \mathbf{u}_p \rangle^p + \rho_f \epsilon \langle \mathbf{u} \rangle^f, \quad (25) $$
fluid phase and the granular phase, obtained at the order $O((d/L)^2)$, for the specific case of the steady flow considered here. These read, respectively,

$$ \rho_f \in \langle \mathbf{u} \rangle^f \cdot \nabla \langle \mathbf{u} \rangle^f = \nabla \cdot \mathbf{\Sigma}^f - n \langle \mathbf{f}_h \rangle^p + \rho_f \in \mathbf{g}, \quad (18) $$

$$ \rho_p \phi \langle \mathbf{u}_p \rangle^p \cdot \nabla \langle \mathbf{u}_p \rangle^p = \nabla \cdot \mathbf{\Sigma}^p + n \langle \mathbf{f}_h \rangle^p + \rho_p \phi \mathbf{g}, \quad (19) $$

which are 2D and for which variables only depend on y. Hence $\nabla = (0, d/dy)$, $\langle \mathbf{u} \rangle^f = (\langle u_x \rangle^f, \langle u_y \rangle^f)$, $\langle \mathbf{u}_p \rangle^p = (\langle u_{p|x} \rangle^p, \langle u_{p|y} \rangle^p)$ and so on. The fluid-particle interaction force denoted $\langle \mathbf{f}_h \rangle^p$ reads

$$ n(y) \langle \mathbf{f}_h \rangle^p(y) = \sum_p \mathcal{G}(|y - y^p|) \int_{S_p} (\boldsymbol{\sigma} \cdot \mathbf{n}) \, dS \equiv \sum_p \mathcal{G}(|y - y^p|) \mathbf{f}_h(\mathbf{x}^p), \quad (20) $$

which can be explicitly calculated from the microscale using Eq. (6). Note that this fluid-particle force is usually split as $n \langle \mathbf{f}_h \rangle^p = n \langle \mathbf{f}_1 \rangle^p + \phi \nabla \cdot \mathbf{\Sigma}^f$, with $\phi \nabla \cdot \mathbf{\Sigma}^f$ the so-called generalized buoyancy. Then $\langle \mathbf{f}_1 \rangle^p$ corresponds to the rest of local hydrodynamics interaction including drag among others.

$\mathbf{\Sigma}^f$ and $\mathbf{\Sigma}^p$ are the 2D apparent stress tensors of the equivalent fluid phase and granular phase, respectively. They can be split into different contributions, defined from microscale, as

$$ \mathbf{\Sigma}^f = \epsilon \langle \boldsymbol{\sigma} \rangle^f + n \langle \mathbf{s}^f \rangle^p - \rho_f \epsilon \langle \mathbf{u}' \otimes \mathbf{u}' \rangle^f, \quad (21) $$

$$ \mathbf{\Sigma}^p = n \langle \mathbf{s}^s \rangle^p - \rho_p \phi \langle \mathbf{u}'_p \otimes \mathbf{u}'_p \rangle^p, \quad (22) $$

where $\mathbf{u}' = (1 - \alpha_{\text{IBM}}(\mathbf{y}))\mathbf{u} - \langle \mathbf{u} \rangle^f$ (respectively, $\mathbf{u}'_p = \mathbf{u}_p - \langle \mathbf{u}_p \rangle^p$) is the spatially fluctuating part of the local fluid (respectively, granular) velocity and the tensors $\langle \mathbf{s}^f \rangle^p$ and $\langle \mathbf{s}^s \rangle^p$ read

$$ n \langle \mathbf{s}^f \rangle^p = \frac{d}{2} \sum_p \mathcal{G}(|y - y^p|) \int_{S_p} (\boldsymbol{\sigma} \cdot \mathbf{n}) \otimes \mathbf{n} \, dS = \frac{d}{2} \sum_p \mathcal{G}(|y - y^p|) \\ \times \int_{V_{\text{tot}}} \left[ \alpha_{\text{IBM}} \rho_f \left( \frac{d\mathbf{u}}{dt} - \mathbf{f}_{\text{IBM}} \right) \otimes \mathbf{r} \right] dv, \quad (23) $$

$$ n \langle \mathbf{s}^s \rangle^p = \frac{d}{2} \sum_p \mathcal{G}(|y - y^p|) \sum_{q \neq p} \mathbf{f}_c^{pq} \otimes \mathbf{n}^{pq}, \quad (24) $$

where $\mathbf{n}^{pq} = (\mathbf{x}^q - \mathbf{x}^p)/(|\mathbf{x}^q - \mathbf{x}^p|)$ is the unit outward vector normal to the surface of the particle $p$ at the contact point.

It shall be noted that the effective fluid stress tensor (21) contains a contribution $n \langle \mathbf{s}^f \rangle^p$ coming from the presence of the grains, leading to a differential stress contribution arising from the vertical variation of the local fluid stress at the particle scale, and probably leading to a $\phi$-dependent effective fluid viscosity as discussed in the following.

Moreover, mass conservation for the fluid phase and the granular phase shall impose, for the shear flow configuration considered here, $\epsilon \langle u_y \rangle^f = 0$ and $\phi \langle u_{p|y} \rangle^p = 0$, respectively. Then, the left-hand sides of Eqs. (18) and (19) have to be zero. We suggest here to keep these terms for now, to verify in the following the validity of averaged momentum equations obtained from the microscale.

### C. Eulerian equations for the mixture

In the previous section, upscaling has been performed in a way that both the fluid phase and the granular phase remain separate. One may also represent the system by a single equivalent phase that contains both the fluid and grains. This phase is here called the mixture phase and its velocity is defined as

$$ \bar{\rho} \langle \mathbf{u} \rangle^m = \rho_p \phi \langle \mathbf{u}_p \rangle^p + \rho_f \epsilon \langle \mathbf{u} \rangle^f, \quad (25) $$
where $\bar{\rho} = \phi \rho_p + \epsilon \rho_f$ is the local density of the mixture. Using Eqs. (18), (19), (25) and mass continuity, one can show that the momentum equation for the mixture phase reads [5]

$$ \bar{\rho} \langle \mathbf{u} \rangle^m \cdot \nabla \langle \mathbf{u} \rangle^m = \nabla \cdot \mathbf{\Sigma}^{\text{tot}} + \bar{\rho} \mathbf{g}, \quad (26) $$

where $\mathbf{\Sigma}^{\text{tot}}$ is the total apparent stress. Again, due to the properties of the configuration studied here, all variables are only y-dependent, and the spatial gradient can be written as $\nabla = (0, d/dy)$.

Mass conservation for the mixture phase implies $d\langle u_y \rangle^m/dy = 0$ leading to $\langle u_y \rangle^m = 0$ over the domain to satisfy upper and bottom boundary conditions, as for the fluid and granular phases (see previous section). This means that Eq. (26) can be simplified as

$$ 0 = \nabla \cdot \mathbf{\Sigma}^{\text{tot}} + \bar{\rho} \mathbf{g}. \quad (27) $$

Then assuming left-hand sides of Eqs. (18) and (19) to be nought, as it shall be, it thus follows that $\mathbf{\Sigma}^{\text{tot}} = \mathbf{\Sigma}^f + \mathbf{\Sigma}^p$ with $\mathbf{\Sigma}^f$ and $\mathbf{\Sigma}^p$ defined as in Eqs. (21) and (22), respectively. The latter equality is not necessarily straightforward as the contribution of the effective stress emanating from the inertial nonlinear term in the mixture phase is not the sum of the contribution of each phase in a general case [5,47]. Note that here, this is not a viscous approximation disregarding velocity differences between granular phase and fluid phase, but that any contribution due to this velocity differences are included in $\mathbf{\Sigma}^f + \mathbf{\Sigma}^p$.

### D. Constitutive laws for the various phases (fluid, granular, mixture)

In the paper, we assume all effective stresses to be symmetric Cauchy-type tensors, as done for generalized Newtonian fluid models, whatever the phase considered (fluid, granular, or mixture). This means that stresses would be the sum of an isotropic pressure and a deviatoric viscous part, characterized by an effective (or equivalent) viscosity. The difference between phases (fluid, granular, mixture) being the structure of this equivalent viscosity and its dependence with state variable, as the solid fraction $\phi$ for instance. This assumption, which is quite strong, allows to test and possibly to extend some of the rheological models obtained in the literature, which are described in the following.

We present here the rheological models that we will use to characterize the different phases. They are not exhaustive but share many similarities with most of the available models of the literature. They are considered here as the base state of our current knowledge on suspension and granular rheology, as simple models for a description as a generalized fluid. Let us first consider the mixture phase. At low particle concentration (and zero particle Reynolds number), Einstein [8,9] showed that the apparent viscosity of the mixture in the presence of a dilute suspension is

$$ \eta_{\text{eq}}^m / \eta_f = 1 + \frac{5}{2} \phi. \quad (28) $$

Reference [48] included small inertial expansion to this solution, highlighting an extra $Re_p^{3/2}$ term into the $\phi$ contribution of Eq. (28). Batchelor and Green (1972) [10] made the second-order development, in the viscous limit, which extends the range of applicability of the model to larger particle concentrations (typically $\phi \le 0.1$). For a shear flow, their model reads

$$ \eta_{\text{eq}}^m / \eta_f = 1 + \frac{5}{2} \phi + 7.6 \phi^2. \quad (29) $$

Models (28) and (29) are exact solutions based on theoretical development, which obviously implies some limitations, as pairwise long-distance interaction and no solid contact. The consequence is that they can only be applied to situations of small solid fraction. For denser situation, the previous assumptions are no longer expected to be valid leading to a significant increase of the dissipation with solid fraction which is not captured by the previous models. In this case, theoretical development of the effective viscosity is not conceivable. Then empirical correlations based on experimental or numerical results have been proposed and can be found in the literature. One of the popular ones for the case of dense and neutrally buoyant suspension is the empirical correlation of
where $\bar{\rho} = \phi \rho_p + \epsilon \rho_f$ is the local density of the mixture. Using Eqs. (18), (19), (25) and mass continuity, one can show that the momentum equation for the mixture phase reads [5]

$$ \bar{\rho} \langle \mathbf{u} \rangle^m \cdot \nabla \langle \mathbf{u} \rangle^m = \nabla \cdot \mathbf{\Sigma}^{\text{tot}} + \bar{\rho} \mathbf{g}, \quad (26) $$

where $\mathbf{\Sigma}^{\text{tot}}$ is the total apparent stress. Again, due to the properties of the configuration studied here, all variables are only y-dependent, and the spatial gradient can be written as $\nabla = (0, d/dy)$.

Mass conservation for the mixture phase implies $d\langle u_y \rangle^m/dy = 0$ leading to $\langle u_y \rangle^m = 0$ over the domain to satisfy upper and bottom boundary conditions, as for the fluid and granular phases (see previous section). This means that Eq. (26) can be simplified as

$$ 0 = \nabla \cdot \mathbf{\Sigma}^{\text{tot}} + \bar{\rho} \mathbf{g}. \quad (27) $$

Then assuming left-hand sides of Eqs. (18) and (19) to be nought, as it shall be, it thus follows that $\mathbf{\Sigma}^{\text{tot}} = \mathbf{\Sigma}^f + \mathbf{\Sigma}^p$ with $\mathbf{\Sigma}^f$ and $\mathbf{\Sigma}^p$ defined as in Eqs. (21) and (22), respectively. The latter equality is not necessarily straightforward as the contribution of the effective stress emanating from the inertial nonlinear term in the mixture phase is not the sum of the contribution of each phase in a general case [5,47]. Note that here, this is not a viscous approximation disregarding velocity differences between granular phase and fluid phase, but that any contribution due to this velocity differences are included in $\mathbf{\Sigma}^f + \mathbf{\Sigma}^p$.

### D. Constitutive laws for the various phases (fluid, granular, mixture)

In the paper, we assume all effective stresses to be symmetric Cauchy-type tensors, as done for generalized Newtonian fluid models, whatever the phase considered (fluid, granular, or mixture). This means that stresses would be the sum of an isotropic pressure and a deviatoric viscous part, characterized by an effective (or equivalent) viscosity. The difference between phases (fluid, granular, mixture) being the structure of this equivalent viscosity and its dependence with state variable, as the solid fraction $\phi$ for instance. This assumption, which is quite strong, allows to test and possibly to extend some of the rheological models obtained in the literature, which are described in the following.

We present here the rheological models that we will use to characterize the different phases. They are not exhaustive but share many similarities with most of the available models of the literature. They are considered here as the base state of our current knowledge on suspension and granular rheology, as simple models for a description as a generalized fluid. Let us first consider the mixture phase. At low particle concentration (and zero particle Reynolds number), Einstein [8,9] showed that the apparent viscosity of the mixture in the presence of a dilute suspension is

$$ \eta_{\text{eq}}^m / \eta_f = 1 + \frac{5}{2} \phi. \quad (28) $$

Reference [48] included small inertial expansion to this solution, highlighting an extra $Re_p^{3/2}$ term into the $\phi$ contribution of Eq. (28). Batchelor and Green (1972) [10] made the second-order development, in the viscous limit, which extends the range of applicability of the model to larger particle concentrations (typically $\phi \le 0.1$). For a shear flow, their model reads

$$ \eta_{\text{eq}}^m / \eta_f = 1 + \frac{5}{2} \phi + 7.6 \phi^2. \quad (29) $$

Models (28) and (29) are exact solutions based on theoretical development, which obviously implies some limitations, as pairwise long-distance interaction and no solid contact. The consequence is that they can only be applied to situations of small solid fraction. For denser situation, the previous assumptions are no longer expected to be valid leading to a significant increase of the dissipation with solid fraction which is not captured by the previous models. In this case, theoretical development of the effective viscosity is not conceivable. Then empirical correlations based on experimental or numerical results have been proposed and can be found in the literature. One of the popular ones for the case of dense and neutrally buoyant suspension is the empirical correlation of
Krieger and Dougherty (1959) [49]

$$ \frac{\eta_{eq}^m}{\eta_f} = \left(1 - \frac{\phi}{\phi_m}\right)^{-\frac{5}{2}\phi_m}, \qquad (30) $$

with $\phi_m = 0.62$. This correlation provides a good approximation but it does not allow to separate the different contributions to the dissipation in the system (viscosity of the fluid, solid frictional contacts).

To understand and theorize the rheology of dense suspensions, the contributions of fluid and particles to momentum transfer have to be separated. For this purpose, Boyer, Guazzelli, and Pouliquen (2011) [12] focused on the rheology of neutrally buoyant suspension with an imposed pressure on the particle phase to control the frictional viscosity usually observed for granular material. Based on independent viscosities for frictional dry granular medium [50] and expected-contactless correlations for the fluid phase, they proposed constitutive laws which unify suspension and granular rheology, modeled as the sum of two contributions, coming from hydrodynamic stresses and contact

$$ \eta_{eq}^m = \eta_{eq}^f + \eta_{eq}^p, \qquad (31) $$

where $\eta_{eq}^f$ and $\eta_{eq}^p$ are, respectively, the fluid and granular contributions. Summation of effective viscosity as in Eq. (31) is not necessarily justified. However, it is indeed reasonable if one assumes the shear rate of each phases to be roughly similar and therefore both equal to the shear rate of the mixture phase. As the total stress is split between each phases, then Eq. (31) holds.

Even if the approach of splitting the two contributions as in Eq. (31) could allow to improve our understanding of physical ingredients emanating from each phase, it still requires to model the effective viscosity of each contribution.

Let us first consider the apparent rheology of the fluid phase. At low concentration, the contribution of the contacts vanish and then, the apparent rheology of the fluid phase equals the one of the mixture. At larger particle concentrations and nonzero particle Reynolds numbers, Gibilaro et al. (2007) [11] propose an empirical expression of the apparent viscosity of the fluid phase, namely,

$$ \eta_{eq}^f / \eta_f = \epsilon^{-\beta} = (1 - \phi)^{-\beta}, \qquad (32) $$

with $\beta = 2.8$. This model can be seen as an empirical extension of models (28) and (29) at larger concentration, still disregarding solid contact between particles. Among other models, it has the advantage of being quite simple, even if it does not satisfy theoretical models when developed at small $\phi$.

Let us finally consider the apparent rheology of the granular phase. According to the present configuration, we focus here on rheological models dealing with dense granular flows of frictional particles, for which the stress is usually assumed to be Coulomb-type

$$ \Sigma_{xy}^p = \mu_{eq}^p P^p, \qquad (33) $$

where $\mu_{eq}^p$ is an effective friction to be modelled and $P^p$ is the granular pressure defined as $P^p = -(\Sigma_{xx}^p + \Sigma_{yy}^p + \Sigma_{zz}^p)/3$. In particular, we make use of the popular $\mu - \mathcal{I}$ rheology which has been extensively developed in the last decade for dry granular flows [39,50,51]. Trulsson et al. (2012) [13] extended this rheology by using both the inertial number

$$ \mathcal{I} = \frac{|\dot{\gamma}^p|d}{\sqrt{P^p / \rho_p}}, \qquad (34) $$

and the viscous number [12,52]

$$ \mathcal{J} = \frac{|\dot{\gamma}^p|\eta_f}{P^p}, \qquad (35) $$
Krieger and Dougherty (1959) [49]

$$ \frac{\eta_{eq}^m}{\eta_f} = \left(1 - \frac{\phi}{\phi_m}\right)^{-\frac{5}{2}\phi_m}, \qquad (30) $$

with $\phi_m = 0.62$. This correlation provides a good approximation but it does not allow to separate the different contributions to the dissipation in the system (viscosity of the fluid, solid frictional contacts).

To understand and theorize the rheology of dense suspensions, the contributions of fluid and particles to momentum transfer have to be separated. For this purpose, Boyer, Guazzelli, and Pouliquen (2011) [12] focused on the rheology of neutrally buoyant suspension with an imposed pressure on the particle phase to control the frictional viscosity usually observed for granular material. Based on independent viscosities for frictional dry granular medium [50] and expected-contactless correlations for the fluid phase, they proposed constitutive laws which unify suspension and granular rheology, modeled as the sum of two contributions, coming from hydrodynamic stresses and contact

$$ \eta_{eq}^m = \eta_{eq}^f + \eta_{eq}^p, \qquad (31) $$

where $\eta_{eq}^f$ and $\eta_{eq}^p$ are, respectively, the fluid and granular contributions. Summation of effective viscosity as in Eq. (31) is not necessarily justified. However, it is indeed reasonable if one assumes the shear rate of each phases to be roughly similar and therefore both equal to the shear rate of the mixture phase. As the total stress is split between each phases, then Eq. (31) holds.

Even if the approach of splitting the two contributions as in Eq. (31) could allow to improve our understanding of physical ingredients emanating from each phase, it still requires to model the effective viscosity of each contribution.

Let us first consider the apparent rheology of the fluid phase. At low concentration, the contribution of the contacts vanish and then, the apparent rheology of the fluid phase equals the one of the mixture. At larger particle concentrations and nonzero particle Reynolds numbers, Gibilaro et al. (2007) [11] propose an empirical expression of the apparent viscosity of the fluid phase, namely,

$$ \eta_{eq}^f / \eta_f = \epsilon^{-\beta} = (1 - \phi)^{-\beta}, \qquad (32) $$

with $\beta = 2.8$. This model can be seen as an empirical extension of models (28) and (29) at larger concentration, still disregarding solid contact between particles. Among other models, it has the advantage of being quite simple, even if it does not satisfy theoretical models when developed at small $\phi$.

Let us finally consider the apparent rheology of the granular phase. According to the present configuration, we focus here on rheological models dealing with dense granular flows of frictional particles, for which the stress is usually assumed to be Coulomb-type

$$ \Sigma_{xy}^p = \mu_{eq}^p P^p, \qquad (33) $$

where $\mu_{eq}^p$ is an effective friction to be modelled and $P^p$ is the granular pressure defined as $P^p = -(\Sigma_{xx}^p + \Sigma_{yy}^p + \Sigma_{zz}^p)/3$. In particular, we make use of the popular $\mu - \mathcal{I}$ rheology which has been extensively developed in the last decade for dry granular flows [39,50,51]. Trulsson et al. (2012) [13] extended this rheology by using both the inertial number

$$ \mathcal{I} = \frac{|\dot{\gamma}^p|d}{\sqrt{P^p / \rho_p}}, \qquad (34) $$

and the viscous number [12,52]

$$ \mathcal{J} = \frac{|\dot{\gamma}^p|\eta_f}{P^p}, \qquad (35) $$
in their model, where $|\dot{\gamma}^p|$ is the norm of the local shear of the granular phase, y-dependent in our situation. This leads to a new dimensionless number

$$ \mathcal{K} = \mathcal{J} + \alpha \mathcal{I}^2, \qquad (36) $$

where $\alpha$ is some constant of order unity. This theory attempts to unify rheological models obtained in viscous ($\mathcal{J}$) and inertial ($\mathcal{I}$) configurations, including the transition from one regime to another. Trulsson et al. (2012) [13] obtained a value of $\alpha = 0.635$ from two-dimensional Euler-Lagrange simulations. Note that Lacaze et al. [42] recently performed fully resolved three-dimensional simulations of the collapse of a granular column in a viscous fluid which seems to support the present value of $\alpha$.

It is worth noting that $\mathcal{K}$ can be expressed as a function of the local Reynold number $Re_l^p = \rho_f |\dot{\gamma}^p| d^2 / \eta_f$ as $\mathcal{K} = \mathcal{J}(1 + \alpha r Re_l^p)$. $\mathcal{K}$ can thus be interpreted as a viscous number $\mathcal{J}$ including some inertial effects via the extra term $\alpha r Re_l^p$. In this study, this extra term was observed to be of order one or more. For instance, in the case $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22), $\alpha r Re_l^p$ ranges from 0 in the static bed of grain up to 1.6 at the top of the moving layer. This means that both regimes are considered here, namely the viscous and the inertial regime, as well as the transition between them.

Following previous results from dry configurations, the authors of Ref. [13] propose models for the effective friction coefficient $\mu_{eq}^p$ and the volume fraction $\phi$ depending on $\mathcal{K}$. As their model for $\phi(\mathcal{K})$ was obtained in the small $\mathcal{K}$ limit, their result can be combined to the viscous model proposed by Ref. [12] to obtain

$$ \phi(\mathcal{K}) = \frac{\phi_c}{1 + b\sqrt{\mathcal{K}}}, \qquad (37) $$

where $\phi_c$ is the jamming volume fraction, which shall be expected to be equivalent to $\phi_m$ in Eq. (30), and the value of $b$ then varies from authors. In particular, $b = 1$ to recover the model of Ref. [12] and $b = 0.42/\phi_c$ for Ref. [13]. The model usually proposes that the friction can be written as

$$ \mu_{eq}^p(\mathcal{K}) = \mu_1 + \frac{\mu_2 - \mu_1}{1 + \sqrt{\mathcal{K}_0/\mathcal{K}}}, \qquad (38) $$

where $\mu_1$, $\mu_2$, and $\mathcal{K}_0$ are fitting parameters. Note that, for instance, Trulsson et al. found $\mu_1 \approx 0.277$, $\mu_2 \approx 0.85$, and $\sqrt{\mathcal{K}_0} \approx 0.29$ from two-dimensional Euler-Lagrange simulations.

Suppressing $\mathcal{K}$ from Eqs. (37) and (38) allows us to express uniquely the friction as a function of the state $\phi$, only valid for steady state,

$$ \mu_{eq}^p(\phi) = \mu_1 + \frac{(\mu_2 - \mu_1)(\phi_c - \phi)}{\phi_c - (1 - b\sqrt{\mathcal{K}_0})\phi}. \qquad (39) $$

Based on this frictional granular model, an effective viscosity for the granular phase can be built as [50]

$$ \eta_{eq}^p(\phi) = \frac{\mu_{eq}^p(\phi)P^p}{|\dot{\gamma}^p|}. \qquad (40) $$

Expressing the pressure term as a function of $\phi$ using Eq. (37), one obtains

$$ \eta_{eq}^p / \eta_f = (1 + \alpha r \mathrm{Re}_l^p) \mu_{eq}^p(\phi) \left( \frac{b \phi}{\phi_c - \phi} \right)^2. \qquad (41) $$

Equation (41) is an extension of the model in Ref. [12] which was derived for the viscous case.

Finally, all effective models described in this section were obtained for viscous flow, even if some validation has been proposed for inertial situations, and for mostly homogenous suspensions. Their validity and extension to inertial and localised shear flows remain sparse and are not yet conclusive.
in their model, where $|\dot{\gamma}^p|$ is the norm of the local shear of the granular phase, y-dependent in our situation. This leads to a new dimensionless number

$$ \mathcal{K} = \mathcal{J} + \alpha \mathcal{I}^2, \qquad (36) $$

where $\alpha$ is some constant of order unity. This theory attempts to unify rheological models obtained in viscous ($\mathcal{J}$) and inertial ($\mathcal{I}$) configurations, including the transition from one regime to another. Trulsson et al. (2012) [13] obtained a value of $\alpha = 0.635$ from two-dimensional Euler-Lagrange simulations. Note that Lacaze et al. [42] recently performed fully resolved three-dimensional simulations of the collapse of a granular column in a viscous fluid which seems to support the present value of $\alpha$.

It is worth noting that $\mathcal{K}$ can be expressed as a function of the local Reynold number $Re_l^p = \rho_f |\dot{\gamma}^p| d^2 / \eta_f$ as $\mathcal{K} = \mathcal{J}(1 + \alpha r Re_l^p)$. $\mathcal{K}$ can thus be interpreted as a viscous number $\mathcal{J}$ including some inertial effects via the extra term $\alpha r Re_l^p$. In this study, this extra term was observed to be of order one or more. For instance, in the case $Re_p = 1$ and $\theta = 0.67$ (Ga = 1.22), $\alpha r Re_l^p$ ranges from 0 in the static bed of grain up to 1.6 at the top of the moving layer. This means that both regimes are considered here, namely the viscous and the inertial regime, as well as the transition between them.

Following previous results from dry configurations, the authors of Ref. [13] propose models for the effective friction coefficient $\mu_{eq}^p$ and the volume fraction $\phi$ depending on $\mathcal{K}$. As their model for $\phi(\mathcal{K})$ was obtained in the small $\mathcal{K}$ limit, their result can be combined to the viscous model proposed by Ref. [12] to obtain

$$ \phi(\mathcal{K}) = \frac{\phi_c}{1 + b\sqrt{\mathcal{K}}}, \qquad (37) $$

where $\phi_c$ is the jamming volume fraction, which shall be expected to be equivalent to $\phi_m$ in Eq. (30), and the value of $b$ then varies from authors. In particular, $b = 1$ to recover the model of Ref. [12] and $b = 0.42/\phi_c$ for Ref. [13]. The model usually proposes that the friction can be written as

$$ \mu_{eq}^p(\mathcal{K}) = \mu_1 + \frac{\mu_2 - \mu_1}{1 + \sqrt{\mathcal{K}_0/\mathcal{K}}}, \qquad (38) $$

where $\mu_1$, $\mu_2$, and $\mathcal{K}_0$ are fitting parameters. Note that, for instance, Trulsson et al. found $\mu_1 \approx 0.277$, $\mu_2 \approx 0.85$, and $\sqrt{\mathcal{K}_0} \approx 0.29$ from two-dimensional Euler-Lagrange simulations.

Suppressing $\mathcal{K}$ from Eqs. (37) and (38) allows us to express uniquely the friction as a function of the state $\phi$, only valid for steady state,

$$ \mu_{eq}^p(\phi) = \mu_1 + \frac{(\mu_2 - \mu_1)(\phi_c - \phi)}{\phi_c - (1 - b\sqrt{\mathcal{K}_0})\phi}. \qquad (39) $$

Based on this frictional granular model, an effective viscosity for the granular phase can be built as [50]

$$ \eta_{eq}^p(\phi) = \frac{\mu_{eq}^p(\phi)P^p}{|\dot{\gamma}^p|}. \qquad (40) $$

Expressing the pressure term as a function of $\phi$ using Eq. (37), one obtains

$$ \eta_{eq}^p / \eta_f = (1 + \alpha r \mathrm{Re}_l^p) \mu_{eq}^p(\phi) \left( \frac{b \phi}{\phi_c - \phi} \right)^2. \qquad (41) $$

Equation (41) is an extension of the model in Ref. [12] which was derived for the viscous case.

Finally, all effective models described in this section were obtained for viscous flow, even if some validation has been proposed for inertial situations, and for mostly homogenous suspensions. Their validity and extension to inertial and localised shear flows remain sparse and are not yet conclusive.
## IV. RESULTS

The results presented in this section were obtained by applying to the DNS the upscaling procedure based on the volume average approach as described in Sec. III. In particular, we use as characteristic length of the weighting function $h_g = d$ as was used for volume fraction and velocities profile in Fig. 1(d). Other values were tested (especially smaller ones) and some results with $h_g \neq d$ are presented in Appendix B. The value of $h_g$ can be crucial in the present configuration. In particular, larger values of $h_g$ could be chosen to satisfy the condition $h_g \gg d/2$. Recall that this condition is required to fulfill the scale separation necessary for the validity of the upscaling approach. However, Fig. 20 shows that increasing $h_g$ tends to underestimate the mesoscale granular velocity close to the granular bed surface. This can lead to an overestimation of the velocity difference, and then the drag, between the two phases (a presentation of the drag force in the granular bed is given in Appendix C). This highlights the main difficulty in performing spatial averaging in such shear-localised configuration. However, decreasing $h_g$, even below $d/2$ allows these profiles to converge. Nevertheless, strong oscillations of the volume fraction profiles, and actually of other mesoscopic quantities, are observed in the bed with an increase of their amplitude when reaching the static bed. Decreasing these oscillations would require to significantly increase the horizontal length scales of the domain. However, it will be shown in the following that this does not modify the physical properties of the sheared mixture. We thus choose in the following to use and discuss main results for $h_g = d$ for sake of clarity, but we will pay attention to discuss its influence according to the results shown in Appendix B.

### A. Momentum balance at the mesoscale

The relevance of the spatial averaging procedure for the configuration considered in this study is first discussed. This is done by verifying the validity of momentum balances at the mesoscale (18) and (19) which shall be satisfied by the equivalent phases, when each terms contributing to these equations are calculated from the resolved IBM simulation, as explained previously. This validity is quantified by two estimated inaccuracies, say a “balance error” $\zeta_b$ and an “integral cumulative error” $\zeta_{ic}$. The former corresponds to the error in the balance of the mesoscale momentum equations. The latter error characterizes the difference between the stress tensor calculated from the average procedure and the estimated stress tensor obtained by integrating momentum equation at the mesoscale.

All the results presented in this section are for $\theta = 0.67$ and $Re_p = 1$ (Ga = 1.22), but can be extended to all cases.

#### 1. The equivalent fluid phase

Let us consider the momentum balance of the equivalent fluid phase [see Eq. (18)] in the vertical direction

$$ \underbrace{\rho_f \epsilon \langle u_y \rangle^f}_{(1)} \frac{\partial \langle u_y \rangle^f}{\partial y} = \epsilon \underbrace{\frac{\partial \Sigma_{yy}^f}{\partial y}}_{(2)} - \underbrace{n \langle f_1 \rangle_y^p}_{(3)} + \underbrace{\rho_f \epsilon g}_{(4)} \quad (42) $$

with

$$ \underbrace{\Sigma_{yy}^f}_{(a)} = \epsilon \underbrace{\langle \sigma_{yy} \rangle^f}_{(b)} + \underbrace{n \langle s_{yy}^f \rangle^p}_{(c)} - \underbrace{\rho_f \epsilon \langle u_y^2 \rangle^f}_{(d)}. \quad (43) $$

Each individual term in Eqs. (42) and (43) are calculated from IBM-simulation results following the volume average procedure. In particular, the obtained vertical profiles of the different terms in Eq. (42) are shown in Fig. 3(a). It appears first that the inertial term (42.1) is negligible. As already mentioned, this proves that mass conservation is actually satisfied, as $\langle u_y \rangle^f$ shall be nought in the present configuration. Moreover, the magnitude of the lift force exerted on the granular phase
## IV. RESULTS

The results presented in this section were obtained by applying to the DNS the upscaling procedure based on the volume average approach as described in Sec. III. In particular, we use as characteristic length of the weighting function $h_g = d$ as was used for volume fraction and velocities profile in Fig. 1(d). Other values were tested (especially smaller ones) and some results with $h_g \neq d$ are presented in Appendix B. The value of $h_g$ can be crucial in the present configuration. In particular, larger values of $h_g$ could be chosen to satisfy the condition $h_g \gg d/2$. Recall that this condition is required to fulfill the scale separation necessary for the validity of the upscaling approach. However, Fig. 20 shows that increasing $h_g$ tends to underestimate the mesoscale granular velocity close to the granular bed surface. This can lead to an overestimation of the velocity difference, and then the drag, between the two phases (a presentation of the drag force in the granular bed is given in Appendix C). This highlights the main difficulty in performing spatial averaging in such shear-localised configuration. However, decreasing $h_g$, even below $d/2$ allows these profiles to converge. Nevertheless, strong oscillations of the volume fraction profiles, and actually of other mesoscopic quantities, are observed in the bed with an increase of their amplitude when reaching the static bed. Decreasing these oscillations would require to significantly increase the horizontal length scales of the domain. However, it will be shown in the following that this does not modify the physical properties of the sheared mixture. We thus choose in the following to use and discuss main results for $h_g = d$ for sake of clarity, but we will pay attention to discuss its influence according to the results shown in Appendix B.

### A. Momentum balance at the mesoscale

The relevance of the spatial averaging procedure for the configuration considered in this study is first discussed. This is done by verifying the validity of momentum balances at the mesoscale (18) and (19) which shall be satisfied by the equivalent phases, when each terms contributing to these equations are calculated from the resolved IBM simulation, as explained previously. This validity is quantified by two estimated inaccuracies, say a “balance error” $\zeta_b$ and an “integral cumulative error” $\zeta_{ic}$. The former corresponds to the error in the balance of the mesoscale momentum equations. The latter error characterizes the difference between the stress tensor calculated from the average procedure and the estimated stress tensor obtained by integrating momentum equation at the mesoscale.

All the results presented in this section are for $\theta = 0.67$ and $Re_p = 1$ ($Ga = 1.22$), but can be extended to all cases.

#### 1. The equivalent fluid phase

Let us consider the momentum balance of the equivalent fluid phase [see Eq. (18)] in the vertical direction

$$ \underbrace{\rho_f \epsilon \langle u_y \rangle^f}_{(1)} \frac{\partial \langle u_y \rangle^f}{\partial y} = \epsilon \underbrace{\frac{\partial \Sigma_{yy}^f}{\partial y}}_{(2)} - \underbrace{n \langle f_1 \rangle_y^p}_{(3)} + \underbrace{\rho_f \epsilon g}_{(4)} \quad (42) $$

with

$$ \underbrace{\Sigma_{yy}^f}_{(a)} = \epsilon \underbrace{\langle \sigma_{yy} \rangle^f}_{(b)} + \underbrace{n \langle s_{yy}^f \rangle^p}_{(c)} - \underbrace{\rho_f \epsilon \langle u_y^2 \rangle^f}_{(d)}. \quad (43) $$

Each individual term in Eqs. (42) and (43) are calculated from IBM-simulation results following the volume average procedure. In particular, the obtained vertical profiles of the different terms in Eq. (42) are shown in Fig. 3(a). It appears first that the inertial term (42.1) is negligible. As already mentioned, this proves that mass conservation is actually satisfied, as $\langle u_y \rangle^f$ shall be nought in the present configuration. Moreover, the magnitude of the lift force exerted on the granular phase
<figure><img src="image_5.png" /></figure>

FIG. 3. Vertical profiles of the vertical projection of the equivalent-fluid momentum balance at the mesoscale (42) and (43), for Re$_{ρ}$ = 1, θ = 0.67, Ga = 1.22. (a) Vertical profiles of all the terms in Eq. (42), labeled from (1) to (4). The black dotted line corresponds to the “balance error” profile $\zeta_b^{f,y}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (43), namely $\Sigma_{yy}^f$ (blue), $\epsilon \langle \sigma_{yy} \rangle^f$ (cyan), $n \langle s_{yy}^f \rangle^p$ (magenta), and $-\rho_f \epsilon \langle u_y^2 \rangle^f$ (yellow). The black dotted line is $\Sigma_{yy}^f$ obtained by vertically integrating Eq. (42), indicating the “integral cumulative” error $\zeta_{ic}^{f,y}$.

[Eq. (42.3); green curve in Fig. 3(a)] is also negligible. Therefore, the momentum balance in the vertical direction for the equivalent fluid phase is ensured by a balance between the buoyancy force (42.2) and the weight of the fluid (42.4), as actually expected.

The “balance error” $\zeta_b^{f,y}$ caused by the upscaling procedure on the vertical momentum equation at the mesoscale is obtained by summing the terms on the right hand side of Eq. (42) [black dotted line in Fig. 3(a)]. $\zeta_b^{f,y}$ is found to be reasonably small, i.e., with a maximum of 7.3% of $\rho_f g$. Not that $\zeta_b^{f,y}$ becomes large close to the bottom rough wall due to the averaging procedure limitation close to the boundaries of the domain. Nevertheless, this part of the domain is not accounted for in the following analysis as it does not influence the dynamics close to the bed surface.

Focusing on the stress contribution $\Sigma_{yy}^f$, Fig. 3(b) presents the obtained vertical profiles of the normal stress $\Sigma_{yy}^f$ and its different contributions obtained as in Eq. (43), namely the fluid phase normal stress $\epsilon \langle \sigma_{yy} \rangle^f$ (43.b), the contribution of the grains to the fluid normal stress $n \langle s_{yy}^f \rangle^p$ (43.c) and the Reynolds stress term $-\rho_f \epsilon \langle v^2 \rangle^f$ (43.d). All contributions are discussed according to their amplitude $|\cdot|$ which all increase from top to bottom due to the weight of the fluid. One first obtains that the normal stress satisfies a hydrostatic balance as expected from the previous discussion on Eq. (42), i.e., $|\Sigma_{yy}^f|$ linearly increases with decreasing y with a slope $\rho_f g$ [blue line (a) in Fig. 3(b)]. Contributions to this hydrostatic balance are as follows. Inertia $-\rho_f \epsilon \langle u_y^2 \rangle^f$ is almost negligible. $|n \langle s_{yy}^f \rangle^p|$ is zero above the bed of grains and increases as one goes deeper in the bed. Alternatively, $|\epsilon \langle \sigma_{yy} \rangle^f|$ linearly increases in the pure fluid column from the top wall down to the granular bed surface. Below the bed surface, $|\epsilon \langle \sigma_{yy} \rangle^f|$ suddenly decreases with decreasing height y in the moving layer while increasing again when reaching the static bed.

A vertical integration of Eq. (42) also leads to $\Sigma_{yy}^f$ estimated from integral values of the other force contributions. Its vertical profile is reported in Fig. 3(b) (black dotted line). The relative error between these two estimations of the normal stress [blue line and black dotted line in Fig. 3(b)] is referred to as its “cumulative integral error” $\zeta_{ic}^{f,y}$. Here, the maximum of $\zeta_{ic}^{f,y}$ is about 3.6% close to the fixed rough wall.

Let us now consider the momentum balance of the equivalent fluid phase in the streamwise direction. For the present flow configuration, the projection of Eq. (18) along the x axis reads

$$ \underbrace{\rho_f \epsilon \langle u_y \rangle^f \frac{\partial \langle u_x \rangle^f}{\partial y}}_{(1)} = \epsilon \underbrace{\frac{\partial \Sigma_{xy}^f}{\partial y}}_{(2)} - \underbrace{n \langle f_1 \rangle^p_x}_{(3)} \quad (44) $$
<figure><img src="image_5.png" /></figure>

FIG. 3. Vertical profiles of the vertical projection of the equivalent-fluid momentum balance at the mesoscale (42) and (43), for Re$_{ρ}$ = 1, θ = 0.67, Ga = 1.22. (a) Vertical profiles of all the terms in Eq. (42), labeled from (1) to (4). The black dotted line corresponds to the “balance error” profile $\zeta_b^{f,y}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (43), namely $\Sigma_{yy}^f$ (blue), $\epsilon(\sigma_{yy})^f$ (cyan), $n(s_{yy}^f)^p$ (magenta), and $-\rho_f \epsilon(u_y^2)^f$ (yellow). The black dotted line is $\Sigma_{yy}^f$ obtained by vertically integrating Eq. (42), indicating the “integral cumulative” error $\zeta_{ic}^{f,y}$.

[Eq. (42.3); green curve in Fig. 3(a)] is also negligible. Therefore, the momentum balance in the vertical direction for the equivalent fluid phase is ensured by a balance between the buoyancy force (42.2) and the weight of the fluid (42.4), as actually expected.

The “balance error” $\zeta_b^{f,y}$ caused by the upscaling procedure on the vertical momentum equation at the mesoscale is obtained by summing the terms on the right hand side of Eq. (42) [black dotted line in Fig. 3(a)]. $\zeta_b^{f,y}$ is found to be reasonably small, i.e., with a maximum of 7.3% of $\rho_f g$. Not that $\zeta_b^{f,y}$ becomes large close to the bottom rough wall due to the averaging procedure limitation close to the boundaries of the domain. Nevertheless, this part of the domain is not accounted for in the following analysis as it does not influence the dynamics close to the bed surface.

Focusing on the stress contribution $\Sigma_{yy}^f$, Fig. 3(b) presents the obtained vertical profiles of the normal stress $\Sigma_{yy}^f$ and its different contributions obtained as in Eq. (43), namely the fluid phase normal stress $\epsilon(\sigma_{yy})^f$ (43.b), the contribution of the grains to the fluid normal stress $n(s_{yy}^f)^p$ (43.c) and the Reynolds stress term $-\rho_f \epsilon(u_y^2)^f$ (43.d). All contributions are discussed according to their amplitude $|\cdot|$ which all increase from top to bottom due to the weight of the fluid. One first obtains that the normal stress satisfies a hydrostatic balance as expected from the previous discussion on Eq. (42), i.e., $|\Sigma_{yy}^f|$ linearly increases with decreasing y with a slope $\rho_f g$ [blue line (a) in Fig. 3(b)]. Contributions to this hydrostatic balance are as follows. Inertia $-\rho_f \epsilon(u_y^2)^f$ is almost negligible. $|n(s_{yy}^f)^p|$ is zero above the bed of grains and increases as one goes deeper in the bed. Alternatively, $|\epsilon(\sigma_{yy})^f|$ linearly increases in the pure fluid column from the top wall down to the granular bed surface. Below the bed surface, $|\epsilon(\sigma_{yy})^f|$ suddenly decreases with decreasing height y in the moving layer while increasing again when reaching the static bed.

A vertical integration of Eq. (42) also leads to $\Sigma_{yy}^f$ estimated from integral values of the other force contributions. Its vertical profile is reported in Fig. 3(b) (black dotted line). The relative error between these two estimations of the normal stress [blue line and black dotted line in Fig. 3(b)] is referred to as its “cumulative integral error” $\zeta_{ic}^{f,y}$. Here, the maximum of $\zeta_{ic}^{f,y}$ is about 3.6% close to the fixed rough wall.

Let us now consider the momentum balance of the equivalent fluid phase in the streamwise direction. For the present flow configuration, the projection of Eq. (18) along the x axis reads

$$ \underbrace{\rho_f \epsilon(u_y)^f \frac{\partial \langle u_x \rangle^f}{\partial y}}_{(1)} = \epsilon \underbrace{\frac{\partial \Sigma_{xy}^f}{\partial y}}_{(2)} - \underbrace{n \langle f_1 \rangle_x^p}_{(3)} \quad (44) $$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Vertical profiles of the horizontal projection of the equivalent-fluid momentum balance at the mesoscale (44) and (45), for Re$_p = 1$, $\theta = 0.67$, $Ga = 1.22$. (a) Vertical profiles of all the terms in Eq. (44), labeled from (1) to (3). The black dotted line corresponds to the “balance error” profile $\zeta_b^{f,x}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (45), namely $\Sigma_{xy}^f$ (blue), $\epsilon\langle\sigma_{xy}\rangle^f$ (cyan), $n\langle s_{xy}^f \rangle^p$ (magenta), and $-\rho_f\epsilon\langle u_x' u_y' \rangle^f$ (yellow). The black dotted line is $\Sigma_{xy}^f$ obtained by vertically integrating Eq. (44), indicating the “integral cumulative” error $\zeta_{ic}^{f,x}$.</figcaption></figure>

with

$$ \underbrace{\Sigma_{xy}^f}_{(a)} = \underbrace{\epsilon\langle\sigma_{xy}\rangle^f}_{(b)} + \underbrace{n\langle s_{xy}^f \rangle^p}_{(c)} - \underbrace{\rho_f\epsilon\langle u_x' u_y' \rangle^f}_{(d)}. \quad (45) $$

Following a similar approach as for the previous discussion on vertical momentum balance, results for the streamwise momentum balance of the fluid phase are presented in Fig. 4. Note that we scaled the forces in the streamwise direction by $\eta_f \Gamma/d$ while in the vertical direction $\rho_f g$ was used, the ratio between these scaling being $(r-1)\theta$. The vertical profiles of the terms of Eq. (44) are shown in Fig. 4(a). As in Fig. 3(a), we plot as black dotted line $\zeta_b^{f,x}$ as an estimation of the local balance error. Figure 4(b) presents the vertical profile of the shear stress $\Sigma_{xy}^f$ (blue) and its three contributions given in Eq. (45). The black dotted line corresponds to the “integral cumulative error” in the streamwise direction $\zeta_{ic}^{f,x}$.

Again, the inertial term in the streamwise direction (44.1) is negligible, as mass balance imposes zero vertical velocity [blue dashed line in Fig. 4(a)]. Unlike vertical balance, the fluid-particle forces only occur due to motion of grains in the streamwise direction, and are therefore nonzero in the bedload region of the moving layers only. In this region, the drag force acting on the grains (44.3) is balanced by the divergence of the shear stress in the fluid (44.2) [see Fig. 4(a)]. As in the vertical direction, the error balance $\zeta_b^{f,x}$ is small but, here, the maximum error is located at the top of the bedload layer. The shear stress $\Sigma_{xy}^f$ profile (45) is shown to decrease from a constant value in the pure fluid phase scaling with $\eta_f \Gamma$ towards zero in the static granular bed, i.e., for decreasing y [blue line (a) in Fig. 4(b)]. As will be discussed later, the equivalent fluid shear stress $\Sigma_{xy}^f$ is actually transferred to the granular one $\Sigma_{xy}^p$ in the granular bed, leading to a constant total shear stress (mixture) in the whole domain, as expected for such configuration (27). The different contributions of this equivalent fluid shear stress profile can also be seen in Fig. 4(b). First, the streamwise contribution of the Reynolds stress is again negligible for the present range of parameters. The main contribution to the shear stress is shown to be the viscous shear stress of the interstitial fluid $\epsilon\langle\sigma_{xy}\rangle^f$ (45.b), at least in the pure fluid layer where it is equal to the equivalent total fluid stress. This contribution then decreases to zero in the static bed, as does the equivalent fluid stress, but it is compensated by an extra contribution $n\langle s_{xy}^f \rangle^p$ (45.c) in the moving granular layer [magenta line in Fig. 4(b)], which can actually becomes larger. This granular contribution shifts, in a way, the influence of the total fluid shear stress $\Sigma_{xy}^f$ deeper in the granular bed, which therefore decreases slower than $\epsilon\langle\sigma_{xy}\rangle^f$. Note that
<figure><img src="image_6.png" /><figcaption>FIG. 4. Vertical profiles of the horizontal projection of the equivalent-fluid momentum balance at the mesoscale (44) and (45), for Re<sub>p</sub> = 1, θ = 0.67, Ga = 1.22. (a) Vertical profiles of all the terms in Eq. (44), labeled from (1) to (3). The black dotted line corresponds to the “balance error” profile ζ<sub>b</sub><sup>f,x</sup> (see text for details). (b) Vertical profiles of all the terms in Eq. (45), namely Σ<sub>xy</sub><sup>f</sup> (blue), ε⟨σ<sub>xy</sub>⟩<sup>f</sup> (cyan), n⟨s<sub>xy</sub><sup>f</sup>⟩<sup>p</sup> (magenta), and −ρ<sub>f</sub>ε⟨u'<sub>x</sub>u'<sub>y</sub>⟩<sup>f</sup> (yellow). The black dotted line is Σ<sub>xy</sub><sup>f</sup> obtained by vertically integrating Eq. (44), indicating the “integral cumulative” error ζ<sub>ic</sub><sup>f,x</sup>.</figcaption></figure>

with

$$ \underbrace{\Sigma_{xy}^f}_{(a)} = \underbrace{\epsilon\langle\sigma_{xy}\rangle^f}_{(b)} + \underbrace{n\langle s_{xy}^f \rangle^p}_{(c)} - \underbrace{\rho_f\epsilon\langle u_x' u_y' \rangle^f}_{(d)}. \quad (45) $$

Following a similar approach as for the previous discussion on vertical momentum balance, results for the streamwise momentum balance of the fluid phase are presented in Fig. 4. Note that we scaled the forces in the streamwise direction by η<sub>f</sub>Γ/d while in the vertical direction ρ<sub>f</sub>g was used, the ratio between these scaling being (r − 1)θ. The vertical profiles of the terms of Eq. (44) are shown in Fig. 4(a). As in Fig. 3(a), we plot as black dotted line ζ<sub>b</sub><sup>f,x</sup> as an estimation of the local balance error. Figure 4(b) presents the vertical profile of the shear stress Σ<sub>xy</sub><sup>f</sup> (blue) and its three contributions given in Eq. (45). The black dotted line corresponds to the “integral cumulative error” in the streamwise direction ζ<sub>ic</sub><sup>f,x</sup>.

Again, the inertial term in the streamwise direction (44.1) is negligible, as mass balance imposes zero vertical velocity [blue dashed line in Fig. 4(a)]. Unlike vertical balance, the fluid-particle forces only occur due to motion of grains in the streamwise direction, and are therefore nonzero in the bedload region of the moving layers only. In this region, the drag force acting on the grains (44.3) is balanced by the divergence of the shear stress in the fluid (44.2) [see Fig. 4(a)]. As in the vertical direction, the error balance ζ<sub>b</sub><sup>f,x</sup> is small but, here, the maximum error is located at the top of the bedload layer. The shear stress Σ<sub>xy</sub><sup>f</sup> profile (45) is shown to decrease from a constant value in the pure fluid phase scaling with η<sub>f</sub>Γ towards zero in the static granular bed, i.e., for decreasing y [blue line (a) in Fig. 4(b)]. As will be discussed later, the equivalent fluid shear stress Σ<sub>xy</sub><sup>f</sup> is actually transferred to the granular one Σ<sub>xy</sub><sup>p</sup> in the granular bed, leading to a constant total shear stress (mixture) in the whole domain, as expected for such configuration (27). The different contributions of this equivalent fluid shear stress profile can also be seen in Fig. 4(b). First, the streamwise contribution of the Reynolds stress is again negligible for the present range of parameters. The main contribution to the shear stress is shown to be the viscous shear stress of the interstitial fluid ε⟨σ<sub>xy</sub>⟩<sup>f</sup> (45.b), at least in the pure fluid layer where it is equal to the equivalent total fluid stress. This contribution then decreases to zero in the static bed, as does the equivalent fluid stress, but it is compensated by an extra contribution n⟨s<sub>xy</sub><sup>f</sup>⟩<sup>p</sup> (45.c) in the moving granular layer [magenta line in Fig. 4(b)], which can actually becomes larger. This granular contribution shifts, in a way, the influence of the total fluid shear stress Σ<sub>xy</sub><sup>f</sup> deeper in the granular bed, which therefore decreases slower than ε⟨σ<sub>xy</sub>⟩<sup>f</sup>. Note that
<figure><img src="image_7.png" /><figcaption>FIG. 5. Vertical profiles of the vertical projection of the equivalent-granular momentum balance at the mesoscale (46) and (47), for Re$_p = 1$, $\theta = 0.67$, $Ga = 1.22$. (a) Vertical profiles of all the terms in Eq. (46), labeled from (1) to (5). The black dotted line corresponds to the “balance error” profile $\zeta_b^{p,y}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (47), namely $\Sigma_{yy}^p$ (red), $n\langle s_{yy}^p \rangle^p$ (magenta), and $-\rho_p\phi\langle u_{p|y}^2 \rangle^f$ (yellow). Note that here $\Sigma_{yy}^p \approx n\langle s_{yy}^p \rangle^p$. The black dotted line is $\Sigma_{yy}^p$ obtained by vertically integrating Eq. (46), indicating the “integral cumulative” error $\zeta_{ic}^{p,y}$.</figcaption></figure>

here, the maximum of $\zeta_{ic}^{f,x}$ is around 3.4% of the total shear stress (black dotted line). This is again due to small discrepancies close to the bed surface which is cumulated when integrating down to the static bed.

## 2. The equivalent granular phase

The procedure to analyze the granular phase is exactly the same as the one described in the previous section and is not recalled here for the sake of conciseness. The momentum balance for the granular phase in the vertical direction (19) reads

$$ \underbrace{\rho_p \phi \langle u_{p|y} \rangle^p}_{(1)} - \underbrace{\frac{\partial \langle u_{p|y} \rangle^p}{\partial y}}_{(2)} + \underbrace{\phi \frac{\partial \Sigma_{yy}^f}{\partial y}}_{(3)} + \underbrace{n \langle f_1 \rangle_y^p}_{(4)} + \underbrace{\rho_p \phi g_y}_{(5)}, \quad (46) $$

with

$$ \Sigma_{yy}^p = n\langle s_{yy}^p \rangle^p - \rho_p \epsilon \langle u_{p|y}^2 \rangle^p. \quad (47) $$

The vertical profiles of the terms (46.1–5) and the relative contribution of the particle stress (47) are plotted in Figs. 5(a) and 5(b), respectively. Again, the inertial term (46.1) vanishes due to mass conservation, i.e., $\langle u_{p|y} \rangle^p = 0$. The apparent weight of the granular phase [grain weight plus generalized buoyancy, i.e., (46.3)+(46.5)] is shown to be balanced by the divergence of the granular normal stress $\Sigma_{yy}^p$ (46.2), often referred to as granular pressure. This implies a granulostatic equilibrium in the granular bed, i.e., a linear profile of $\Sigma_{yy}^p$ in the bed which becomes zero in the pure fluid upper layer [see Fig. 5(b)]. Note that only particle-particle contacts contribute to this granular stress, as the inertial term $\rho_p \epsilon \langle u_{p|y}^2 \rangle^p$ in the particle stress term (47) is found negligible here. Again, the vertical contribution of the fluid-particle force $n \langle f_1 \rangle_y^p$ (46.4) is found to be negligible, even in the moving layer. Errors $\zeta_b^{p,y}$ and $\zeta_{ic}^{p,y}$ estimated as for the fluid phase remain small in the moving layer of grains. However, they can significantly increase in the static bed [see black dotted lines in Figs. 5(a) and 5(b)]. The reason for that inaccuracy has already been mentioned previously as a limitation in upscaling approach using the present procedure in the static bed, which is not a limitation according to the purpose of the paper.
<figure><img src="image_7.png" /><figcaption>FIG. 5. Vertical profiles of the vertical projection of the equivalent-granular momentum balance at the mesoscale (46) and (47), for Re$_p = 1$, $\theta = 0.67$, $Ga = 1.22$. (a) Vertical profiles of all the terms in Eq. (46), labeled from (1) to (5). The black dotted line corresponds to the “balance error” profile $\zeta_b^{p,y}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (47), namely $\Sigma_{yy}^p$ (red), $n\langle s_{yy}^p \rangle^p$ (magenta), and $-\rho_p\phi\langle u_{p|y}^2 \rangle^f$ (yellow). Note that here $\Sigma_{yy}^p \approx n\langle s_{yy}^p \rangle^p$. The black dotted line is $\Sigma_{yy}^p$ obtained by vertically integrating Eq. (46), indicating the “integral cumulative” error $\zeta_{ic}^{p,y}$.</figcaption></figure>

here, the maximum of $\zeta_{ic}^{f,x}$ is around 3.4% of the total shear stress (black dotted line). This is again due to small discrepancies close to the bed surface which is cumulated when integrating down to the static bed.

## 2. The equivalent granular phase

The procedure to analyze the granular phase is exactly the same as the one described in the previous section and is not recalled here for the sake of conciseness. The momentum balance for the granular phase in the vertical direction (19) reads

$$ \underbrace{\rho_p \phi \langle u_{p|y} \rangle^p}_{(1)} - \underbrace{\frac{\partial \langle u_{p|y} \rangle^p}{\partial y}}_{(2)} + \underbrace{\phi \frac{\partial \Sigma_{yy}^f}{\partial y}}_{(3)} + \underbrace{n \langle f_1 \rangle_y^p}_{(4)} + \underbrace{\rho_p \phi g_y}_{(5)}, \quad (46) $$

with

$$ \Sigma_{yy}^p = n\langle s_{yy}^p \rangle^p - \rho_p \epsilon \langle u_{p|y}^2 \rangle^p. \quad (47) $$

The vertical profiles of the terms (46.1–5) and the relative contribution of the particle stress (47) are plotted in Figs. 5(a) and 5(b), respectively. Again, the inertial term (46.1) vanishes due to mass conservation, i.e., $\langle u_{p|y} \rangle^p = 0$. The apparent weight of the granular phase [grain weight plus generalized buoyancy, i.e., (46.3)+(46.5)] is shown to be balanced by the divergence of the granular normal stress $\Sigma_{yy}^p$ (46.2), often referred to as granular pressure. This implies a granulostatic equilibrium in the granular bed, i.e., a linear profile of $\Sigma_{yy}^p$ in the bed which becomes zero in the pure fluid upper layer [see Fig. 5(b)]. Note that only particle-particle contacts contribute to this granular stress, as the inertial term $\rho_p \epsilon \langle u_{p|y}^2 \rangle^p$ in the particle stress term (47) is found negligible here. Again, the vertical contribution of the fluid-particle force $n \langle f_1 \rangle_y^p$ (46.4) is found to be negligible, even in the moving layer. Errors $\zeta_b^{p,y}$ and $\zeta_{ic}^{p,y}$ estimated as for the fluid phase remain small in the moving layer of grains. However, they can significantly increase in the static bed [see black dotted lines in Figs. 5(a) and 5(b)]. The reason for that inaccuracy has already been mentioned previously as a limitation in upscaling approach using the present procedure in the static bed, which is not a limitation according to the purpose of the paper.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Vertical profiles of the horizontal projection of the equivalent-granular momentum balance at the mesoscale (48) and (49), for Re_p = 1, θ = 0.67, Ga = 1.22. (a) Vertical profiles of all the terms in Eq. (48), labeled from (1) to (4). The black dotted line corresponds to the “balance error” profile ξ_b^{p,x} (see text for details). (b) Vertical profiles of all the terms in Eq. (49), namely Σ_{xy}^p (red), n⟨s_{yy}^p⟩^p (magenta), and −ρ_pφ⟨u_{p|y}^2⟩^f (yellow). Note that here Σ_{xy}^p ≈ n⟨s_{xy}^p⟩^p. The black dotted line is Σ_{xy}^p obtained by vertically integrating (48), indicating the “integral cumulative” error ξ_{ic}^{p,x}.</figcaption></figure>

The momentum balance in the streamwise direction (19) of the equivalent granular phase reads

$$ \underbrace{\rho_p \langle u_{p|y} \rangle^p}_{(1)} \frac{\partial \langle u_{p|x} \rangle^p}{\partial y} = \underbrace{\frac{\partial \Sigma_{xy}^p}{\partial y}}_{(2)} + \underbrace{\phi \frac{\partial \Sigma_{xy}^f}{\partial y}}_{(3)} + \underbrace{n \langle f_1 \rangle_x^p}_{(4)} \quad (48) $$

with

$$ \Sigma_{xy}^{p} = n \langle s_{xy}^{p} \rangle^{p} - \rho_{p} \epsilon \langle u'_{p|x} u'_{p|y} \rangle^{p}. \quad (49) $$

The vertical profiles of Eqs. (48) and (49) are shown in Figs. 6(a) and 6(b), respectively. Similarly to the vertical balance, Eq. (48.1) vanishes. Then, dissipation term associated with the vertical gradient of the shear stress (48.2) is balanced by the sum of the streamwise fluid-particle drag force (48.4) and the generalized buoyancy gradient (48.3), the latter playing also as a driving term emanating from the shear at the particle scale. Note that the drag contribution is dominant in this bedload moving layer and contributes to about 60% of the driving force in the streamwise direction. However, the presence of the buoyancy gradient, in the sense of Ref. [5], on the force balance in the streamwise direction, has to be associated with localised shear contributions at the particle scale, leading to significant stress variation from the top to the bottom of each grain. To finish with, the shear stress Σ_{xy}^p, whose contribution only comes from particle-particle contact [first term of Eq. (49)] as shown in Fig. 6(b), is zero in the pure fluid top layer while increasing with decreasing y in the bedload moving layer and reaches a roughly constant value of ≈η_fΓ in the static bed region. Hence, the fluid shear stress in the pure fluid region has been fully transferred, as expected and mentioned previously, to the granular phase in the static bed via the bedload moving layer. Again ξ_b^{p,x} and ξ_{ic}^{p,x} remain small only in the moving layer of grains.

To conclude on the mesoscale momentum equations estimated here, the expected balance between the two phases are obtained in both the vertical direction and the horizontal direction. Moreover, such approach allows to extract the main contributions to the fluid stress and the particle stress emanating from the fluid-particle interaction. It is worth noting that, a significant discrepancy in momentum balance for the granular phase, ξ_b^p and ξ_{ic}^p, has been observed in the whole static bed region. The latter is due to the time averaging of correlated samples when particles remain static (in particular when we time average forces for particles at the exact same location). Those fluctuations, however, are limited in amplitude and localized in the static bed (a more detailed discussion is given in Appendix A). Note that rheological flow curves for an equivalent generalized Newtonian
<figure><img src="image_8.png" /><figcaption>FIG. 6. Vertical profiles of the horizontal projection of the equivalent-granular momentum balance at the mesoscale (48) and (49), for Re_p = 1, θ = 0.67, Ga = 1.22. (a) Vertical profiles of all the terms in Eq. (48), labeled from (1) to (4). The black dotted line corresponds to the “balance error” profile ξ_b^{p,x} (see text for details). (b) Vertical profiles of all the terms in Eq. (49), namely Σ_{xy}^p (red), n⟨s_{yy}^p⟩^p (magenta), and −ρ_pφ⟨u_{p|y}^2⟩^f (yellow). Note that here Σ_{xy}^p ≈ n⟨s_{xy}^p⟩^p. The black dotted line is Σ_{xy}^p obtained by vertically integrating (48), indicating the “integral cumulative” error ζ_{ic}^{p,x}.</figcaption></figure>

The momentum balance in the streamwise direction (19) of the equivalent granular phase reads

$$ \underbrace{\rho_p \langle u_{p|y} \rangle^p}_{(1)} \frac{\partial \langle u_{p|x} \rangle^p}{\partial y} = \underbrace{\frac{\partial \Sigma_{xy}^p}{\partial y}}_{(2)} + \underbrace{\phi \frac{\partial \Sigma_{xy}^f}{\partial y}}_{(3)} + \underbrace{n \langle f_1 \rangle_x^p}_{(4)} \quad (48) $$

with

$$ \Sigma_{xy}^{p} = n \langle s_{xy}^{p} \rangle^{p} - \rho_{p} \epsilon \langle u'_{p|x} u'_{p|y} \rangle^{p}. \quad (49) $$

The vertical profiles of Eqs. (48) and (49) are shown in Figs. 6(a) and 6(b), respectively. Similarly to the vertical balance, Eq. (48.1) vanishes. Then, dissipation term associated with the vertical gradient of the shear stress (48.2) is balanced by the sum of the streamwise fluid-particle drag force (48.4) and the generalized buoyancy gradient (48.3), the latter playing also as a driving term emanating from the shear at the particle scale. Note that the drag contribution is dominant in this bedload moving layer and contributes to about 60% of the driving force in the streamwise direction. However, the presence of the buoyancy gradient, in the sense of Ref. [5], on the force balance in the streamwise direction, has to be associated with localised shear contributions at the particle scale, leading to significant stress variation from the top to the bottom of each grain. To finish with, the shear stress Σ_{xy}^p, whose contribution only comes from particle-particle contact [first term of Eq. (49)] as shown in Fig. 6(b), is zero in the pure fluid top layer while increasing with decreasing y in the bedload moving layer and reaches a roughly constant value of ≈η_fΓ in the static bed region. Hence, the fluid shear stress in the pure fluid region has been fully transferred, as expected and mentioned previously, to the granular phase in the static bed via the bedload moving layer. Again ζ_b^{p,x} and ζ_{ic}^{p,x} remain small only in the moving layer of grains.

To conclude on the mesoscale momentum equations estimated here, the expected balance between the two phases are obtained in both the vertical direction and the horizontal direction. Moreover, such approach allows to extract the main contributions to the fluid stress and the particle stress emanating from the fluid-particle interaction. It is worth noting that, a significant discrepancy in momentum balance for the granular phase, ζ_b^p and ζ_{ic}^p, has been observed in the whole static bed region. The latter is due to the time averaging of correlated samples when particles remain static (in particular when we time average forces for particles at the exact same location). Those fluctuations, however, are limited in amplitude and localized in the static bed (a more detailed discussion is given in Appendix A). Note that rheological flow curves for an equivalent generalized Newtonian
fluid as discussed in Sec. III D are only valid above the threshold of motion, i.e., in the moving layer here. Fluctuations in the static bed then do not affect such rheological description, and will be disregarded. Beyond these observations, discrepancies in momentum balance $\zeta_b^f$, $\zeta_{ic}^f$, $\zeta_b^p$, $\zeta_{ic}^p$ pointed out in Figs. 3–6 remain small in the moving layer. An improvement could however be expected by considering larger-order terms in the asymptotic expansion used to derive the mesoscopic equations (see Sec. III B). According to the small influence of these terms, such analyses remain beyond the scope of the present paper.

## B. Rheological models at the mesoscale

In this section, we discuss the rheology of the different equivalent phases computed from the microscale simulations, IBM-DEM. In all cases, to prevent over-interpretation of our results, we exclude the data coming from the static zone ($y \lesssim h_{bed} - 4d$), which corresponds to $\phi \gtrsim 0.601$ and $K \lesssim 10^{-5}$, and where fluctuations of the mesoscale quantities were observed in the mean vertical profiles (see related discussion at the end of the previous section).

### 1. The equivalent fluid phase

First consider the apparent rheology of the equivalent fluid phase. Assuming that the equivalent fluid phase is a generalized Netwonian fluid, its apparent viscosity $\eta_{eq}^f$ is a scalar which can be determined as $\eta_{eq}^f = \Sigma_{xy}^f / (\partial \langle u_x \rangle^f / \partial y) = \Sigma_{xy}^f / \dot{\gamma}^f$ for the 1D shear flow considered here. For convenience, it is usually preferred to use $(\Sigma_{xy}^f + \Sigma_{yx}^f)/2$ instead of $\Sigma_{xy}^f$ for symmetry reasons (even if those two terms are almost equal, the maximum difference being of the order of $0.02 \eta_f \dot{\Gamma}$ in our case). According to the previous definition of the effective fluid viscosity, the equivalent fluid stress is then proportional to the effective fluid strain rate. However, this link is far from being obvious, and remains an open question. In particular, using this definition for the effective stress would require to define a bulk viscosity according to the compressibility of the equivalent fluid phase. Another approach, which has been commonly used in the literature according to solution derived for dilute suspension in a viscous Stokes regime, is to make use of the incompressibility of the mixture phase $\langle \mathbf{u} \rangle = \phi \langle \mathbf{u}_p \rangle^p + \epsilon \langle \mathbf{u} \rangle^f$ to define the stress as a function of its strain rate even for the equivalent fluid phase [6,53,54]. Note that such a definition is also convenient as we anticipate defining the mixture phase viscosity as the sum of the equivalent phases contribution (31). Then, we use the following definition for computing the apparent viscosity of the equivalent fluid phase

$$ \eta_{eq}^f = \frac{\frac{1}{2}(\Sigma_{xy}^f + \Sigma_{yx}^f)}{\frac{\partial \langle u_x \rangle}{\partial y}} \equiv \frac{\frac{1}{2}(\Sigma_{xy}^f + \Sigma_{yx}^f)}{\dot{\gamma}}. \qquad (50) $$

More specifically, from the vertical evolution of $\Sigma_{xy}^f$, $\Sigma_{yx}^f$ and $\langle u_x \rangle$, we compute $\eta_{eq}^f(y)/\eta_f$. As the solid fraction also depends on the vertical position, $\phi(y)$ [see, e.g., the black line in Fig. 1(d)], one obtains $\eta_{eq}^f/\eta_f$ as a single-valued function of $\phi$, as shown in Fig. 7 (dot symbols) for $Re_p = 1$, $\theta = 0.67$ (Ga = 1.22).

The apparent viscosity of the fluid phase obtained in this latter case is very well fitted by the power law (32) with $\beta = 1.71$, at least for $\phi < 0.55$. Note that this value is rather different from $\beta = 2.8$ as obtained by Ref. [11], and do not moreover lead to the Einstein viscosity (28) at small $\phi$.

Figure 8(a) is similar to Fig. 7 but displays the apparent viscosity of the fluid phase for all the cases considered here, i.e., $Re_p = \{0.1, 1, 10\}$ and $\theta \in [0.1, 0.7]$. One observes that the rheological curve strongly varies with these two parameters. In any case, $\eta_{eq}^f/\eta_f$ can be fitted with a law of the form $(1 - \phi)^{-\beta}$ as previously, but now with $\beta = \beta(\theta, Re_p)$. In particular, increasing $\theta$ or $Re_p$ increase the effective fluid viscosity. In Fig. 8(a), green solid lines show the obtained fit for $(Re_p, \theta, Ga) = (0.1, 0.22, 0.67)$ and $(Re_p, \theta, Ga) = (10, 0.67, 3.9)$ corresponding to $\beta = 0.83$ and $\beta = 3.24$, respectively.
fluid as discussed in Sec. III D are only valid above the threshold of motion, i.e., in the moving layer here. Fluctuations in the static bed then do not affect such rheological description, and will be disregarded. Beyond these observations, discrepancies in momentum balance $\zeta_b^f$, $\zeta_{ic}^f$, $\zeta_b^p$, $\zeta_{ic}^p$ pointed out in Figs. 3–6 remain small in the moving layer. An improvement could however be expected by considering larger-order terms in the asymptotic expansion used to derive the mesoscopic equations (see Sec. III B). According to the small influence of these terms, such analyses remain beyond the scope of the present paper.

## B. Rheological models at the mesoscale

In this section, we discuss the rheology of the different equivalent phases computed from the microscale simulations, IBM-DEM. In all cases, to prevent over-interpretation of our results, we exclude the data coming from the static zone ($y \lesssim h_{\text{bed}} - 4d$), which corresponds to $\phi \gtrsim 0.601$ and $K \lesssim 10^{-5}$, and where fluctuations of the mesoscale quantities were observed in the mean vertical profiles (see related discussion at the end of the previous section).

### 1. The equivalent fluid phase

First consider the apparent rheology of the equivalent fluid phase. Assuming that the equivalent fluid phase is a generalized Netwonian fluid, its apparent viscosity $\eta_{eq}^f$ is a scalar which can be determined as $\eta_{eq}^f = \Sigma_{xy}^f/(\partial\langle u_x \rangle^f/\partial y) = \Sigma_{xy}^f/\dot{\gamma}^f$ for the 1D shear flow considered here. For convenience, it is usually preferred to use $(\Sigma_{xy}^f + \Sigma_{yx}^f)/2$ instead of $\Sigma_{xy}^f$ for symmetry reasons (even if those two terms are almost equal, the maximum difference being of the order of $0.02 \eta_f \dot{\Gamma}$ in our case). According to the previous definition of the effective fluid viscosity, the equivalent fluid stress is then proportional to the effective fluid strain rate. However, this link is far from being obvious, and remains an open question. In particular, using this definition for the effective stress would require to define a bulk viscosity according to the compressibility of the equivalent fluid phase. Another approach, which has been commonly used in the literature according to solution derived for dilute suspension in a viscous Stokes regime, is to make use of the incompressibility of the mixture phase $\langle \mathbf{u} \rangle = \phi \langle \mathbf{u}_p \rangle^p + \epsilon \langle \mathbf{u} \rangle^f$ to define the stress as a function of its strain rate even for the equivalent fluid phase [6,53,54]. Note that such a definition is also convenient as we anticipate defining the mixture phase viscosity as the sum of the equivalent phases contribution (31). Then, we use the following definition for computing the apparent viscosity of the equivalent fluid phase

$$ \eta_{eq}^f = \frac{\frac{1}{2}(\Sigma_{xy}^f + \Sigma_{yx}^f)}{\frac{\partial \langle u_x \rangle}{\partial y}} \equiv \frac{\frac{1}{2}(\Sigma_{xy}^f + \Sigma_{yx}^f)}{\dot{\gamma}}. \quad (50) $$

More specifically, from the vertical evolution of $\Sigma_{xy}^f$, $\Sigma_{yx}^f$ and $\langle u_x \rangle$, we compute $\eta_{eq}^f(y)/\eta_f$. As the solid fraction also depends on the vertical position, $\phi(y)$ [see, e.g., the black line in Fig. 1(d)], one obtains $\eta_{eq}^f/\eta_f$ as a single-valued function of $\phi$, as shown in Fig. 7 (dot symbols) for $Re_p = 1$, $\theta = 0.67$ (Ga = 1.22).

The apparent viscosity of the fluid phase obtained in this latter case is very well fitted by the power law (32) with $\beta = 1.71$, at least for $\phi < 0.55$. Note that this value is rather different from $\beta = 2.8$ as obtained by Ref. [11], and do not moreover lead to the Einstein viscosity (28) at small $\phi$.

Figure 8(a) is similar to Fig. 7 but displays the apparent viscosity of the fluid phase for all the cases considered here, i.e., $Re_p = \{0.1, 1, 10\}$ and $\theta \in [0.1, 0.7]$. One observes that the rheological curve strongly varies with these two parameters. In any case, $\eta_{eq}^f/\eta_f$ can be fitted with a law of the form $(1 - \phi)^{-\beta}$ as previously, but now with $\beta = \beta(\theta, Re_p)$. In particular, increasing $\theta$ or $Re_p$ increase the effective fluid viscosity. In Fig. 8(a), green solid lines show the obtained fit for $(Re_p, \theta, \text{Ga}) = (0.1, 0.22, 0.67)$ and $(Re_p, \theta, \text{Ga}) = (10, 0.67, 3.9)$ corresponding to $\beta = 0.83$ and $\beta = 3.24$, respectively.
<figure><img src="image_9.png" /><figcaption>FIG. 7. Apparent viscosity $\eta_{eq}^f/\eta_f$ of the fluid phase vs solid volume fraction $\phi$: (●), present IBM-DEM simulation for $Re_p = 1$, $\theta = 0.67$, Ga = 1.22 [Eq. (50)]; (dashed black) Einstein's model [8,9] $\eta_{eq}^f/\eta_f = 1 + \frac{5}{3}\phi$; (solid black) Gibilaro et al. (2007)'s correlation [11] $\eta_{eq}^f/\eta_f = (1-\phi)^{-\beta}$ with $\beta = 2.8$; and (red) adjustment with $\eta_{eq}^f/\eta_f = (1-\phi)^{-1.71}$. Confidence intervals of 95% are computed using the geometric standard deviation of the apparent viscosity for each uncorrelated time sample.</figcaption></figure>

Figure 8(b) presents the corresponding value of the $\beta$ exponent used in the constitutive law $\eta_{eq}^f/\eta_f = (1-\phi)^{-\beta}$ as a function of $\theta$ and $Re_p$ and obtained from a best fit of the data of Fig. 8(a). As already mentioned, the exponent $\beta$ obtained here strongly depends on the particle Reynolds number and the Shields number. It is found that $\beta$ has a power-law dependence on both $Re_p$ and $\theta$. In particular, the $\theta$-dependency seems linear and the $Re_p$-dependency is well described by a square root function. Note that at small $Re_p$ (black symbols for $Re_p = 0.1$ here), $\beta$ becomes nearly independent of $\theta$. To summarize and to highlight these observations, Fig. 8(c) shows the evolution

<figure><img src="image_10.png" /><figcaption>FIG. 8. (a) Apparent viscosity $\eta_{eq}^f/\eta_f$ of the fluid phase vs solid volume fraction $\phi$ for all the cases ($Re_p, \theta$, Ga) considered in the present work: see Table I for detail of the various symbols. Green solid lines show the obtained fit of the constitutive law $\eta_{eq}^f/\eta_f = (1-\phi)^{-\beta}$ for ($Re_p, \theta$, Ga) = (0.1, 0.67, 0.22) (black crosses) and ($Re_p, \theta$, Ga) = (10, 0.67, 3.9) (red diamonds), respectively. Green dashed lines in the insert are the small-$\phi$ expansion (52). (b) Corresponding exponent $\beta$ as a function of $\theta$. (c) $\beta - \beta_0$ as a function of $r\theta\sqrt{Re_p}$ with $\beta_0 = 0.565$. Lines correspond to 2/3 power-law functions with scale factor 1 (dotted line) and 0.85 (dashed line).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Apparent viscosity $\eta_{eq}^f/\eta_f$ of the fluid phase vs solid volume fraction $\phi$: (●), present IBM-DEM simulation for $Re_p = 1$, $\theta = 0.67$, Ga = 1.22 [Eq. (50)]; (dashed black) Einstein's model [8,9] $\eta_{eq}^f/\eta_f = 1 + \frac{5}{3}\phi$; (solid black) Gibilaro et al. (2007)'s correlation [11] $\eta_{eq}^f/\eta_f = (1-\phi)^{-\beta}$ with $\beta = 2.8$; and (red) adjustment with $\eta_{eq}^f/\eta_f = (1-\phi)^{-1.71}$. Confidence intervals of 95% are computed using the geometric standard deviation of the apparent viscosity for each uncorrelated time sample.</figcaption></figure>

Figure 8(b) presents the corresponding value of the $\beta$ exponent used in the constitutive law $\eta_{eq}^f/\eta_f = (1-\phi)^{-\beta}$ as a function of $\theta$ and $Re_p$ and obtained from a best fit of the data of Fig. 8(a). As already mentioned, the exponent $\beta$ obtained here strongly depends on the particle Reynolds number and the Shields number. It is found that $\beta$ has a power-law dependence on both $Re_p$ and $\theta$. In particular, the $\theta$-dependency seems linear and the $Re_p$-dependency is well described by a square root function. Note that at small $Re_p$ (black symbols for $Re_p = 0.1$ here), $\beta$ becomes nearly independent of $\theta$. To summarize and to highlight these observations, Fig. 8(c) shows the evolution

<figure><img src="image_10.png" /><figcaption>FIG. 8. (a) Apparent viscosity $\eta_{eq}^f/\eta_f$ of the fluid phase vs solid volume fraction $\phi$ for all the cases ($Re_p, \theta, \text{Ga}$) considered in the present work: see Table I for detail of the various symbols. Green solid lines show the obtained fit of the constitutive law $\eta_{eq}^f/\eta_f = (1-\phi)^{-\beta}$ for ($Re_p, \theta, \text{Ga}$) = (0.1, 0.67, 0.22) (black crosses) and ($Re_p, \theta, \text{Ga}$) = (10, 0.67, 3.9) (red diamonds), respectively. Green dashed lines in the insert are the small-$\phi$ expansion (52). (b) Corresponding exponent $\beta$ as a function of $\theta$. (c) $\beta - \beta_0$ as a function of $r\theta\sqrt{Re_p}$ with $\beta_0 = 0.565$. Lines correspond to 2/3 power-law functions with scale factor 1 (dotted line) and 0.85 (dashed line).</figcaption></figure>
of $\beta - \beta_0$ as a function of $r\theta\sqrt{\text{Re}_p}$ with $\beta_0 = 0.565$ an offset constant value. One finally obtains

$$ \beta = \begin{cases} \beta_0 + a (r \theta \sqrt{\text{Re}_p})^{2/3}, \\ \beta_0 + a r^{2/3} \theta \text{Ga}^{2/3}, \\ \beta_0 + a r^{2/3} \text{Re}_p \text{Ga}^{-4/3}, \end{cases} \quad (51) $$

depending on the set of the three independent dimensionless numbers chosen. In Eq. (51), $a$ is found to be $a = 0.85$ for the explored values of $\theta$ and $\text{Re}_p$ [dashed line in Fig. 8(c), dotted line corresponds to $a = 1$]. Note in particular that this scaling is not expected to hold at large $\text{Re}_p$ and $\theta$ for which different mechanisms probably become dominant. To understand the obtain scaling law (51), one can discuss its implication at small $\phi$. If $\phi \ll 1$, combining Eqs. (32) and (51) leads to

$$ \begin{aligned} \eta_{\text{eq}}^f / \eta_f &= 1 + \phi [\beta_0 + a (r \theta \sqrt{\text{Re}_p})^{2/3}], \\ &\equiv 1 + \beta_0 \phi [1 + l_g^2 \rho_f \dot{\Gamma} / \eta_f], \\ &= 1 + \beta_0 \phi \left[1 + \left(\frac{l_g}{d}\right)^2 \text{Re}_p\right], \end{aligned} \quad (52) $$

where $l_g = \sqrt{a/\beta_0} r^{1/3} \text{Ga}^{-2/3} d$ is a length scale. Note that we found $l_g$ in the range $0.4 \lesssim l_g/d \lesssim 3$ in our simulations. Green dashed lines in the insert of 8(a) highlight the relevance of such small-$\phi$ expansion, again for $(\text{Re}_p, \theta, \text{Ga}) = (0.1, 0.22, 0.67)$ and $(\text{Re}_p, \theta, \text{Ga}) = (10, 0.67, 3.9)$. By analogy with the eddy viscosity in turbulent flow, in which $l_g$ would be a mixing length, here it characterizes a viscous length scale associated with the vertical motion of grains induced by their own weight along the vertical direction during transport. This shall therefore somehow be linked to the thickness of the viscous boundary layer developing at the surface of these grains because of relative motion with the fluid along the vertical direction. In other words, this “eddy” viscosity is induced by vertical agitation of the grains moving at the surface of the granular bed and subjected to the characteristic shear flow $\dot{\Gamma}$. This contribution to the fluid-phase dissipation is therefore associated with grain weight. This probably explain why the obtained $\text{Re}_p$ contribution differs from Ref. [48], which found a $\text{Re}_p^{3/2}$ correction in the case of dilute and homogeneous neutrally buoyant suspension in a simple shear flow, i.e., in which inertial contribution comes from other mechanisms.

The value of $\beta_0$ then specifies the viscous limit for the effective viscosity of the equivalent fluid phase in the dilute regime, $\phi \ll 1$. It is observed to be quite smaller than 2.5 which would be the expected value obtained for homogeneous suspension. However, it shall be noted that such viscous limit has usually been obtained for nearly neutrally buoyant suspension, i.e., for $\text{Ga} \to 0$ and thus $l_g/d \to \infty$ here. Such limit is thus undefined with the obtained scaling Eqs. (51) and (52), as this would actually be the $\theta \to \infty$ limit, which is not relevant for bedload. Here, the $\phi$-dependency at small $\phi$ is a very different process, and is therefore not comparable to a dilute neutrally buoyant suspension. The dissipation in the fluid phase remains associated with the specific dynamics of heavy grains in close contact, which can modify a purely 1D-vertical diffusive process due to vertical agitation of grains.

Moreover, in the finite range of parameters $(\theta, \text{Ga}, \text{Re}_p)$ covered here for bedload application, the transition from the top fluid layer towards the granular bed remains a sharp interface at the grain scale (as the bed remains relatively dense, even if $\phi$ varies within this bed layer as discussed in the next section). The present spatial upscaling is performed across the entire system and hence it leads to an equivalent diffuse interface as grains remain microscale. The smooth $\phi$-transition observed in Fig. 8(a) is therefore the consequence of this diffuse interface between the bed and the fluid layer at the continuous scale. Such diffuse interface includes the dissipative mechanisms associated with heavy grains mentioned above. $\theta \to 0$ corresponds to the sharpest interface, i.e., for a fluid flow over a porous medium built of solid grains. In the range of $\text{Re}_p$ covered here, and according to (51),
of $\beta - \beta_0$ as a function of $r\theta\sqrt{\text{Re}_p}$ with $\beta_0 = 0.565$ an offset constant value. One finally obtains

$$ \beta = \begin{cases} \beta_0 + a (r \theta \sqrt{\text{Re}_p})^{2/3}, \\ \beta_0 + a r^{2/3} \theta \text{Ga}^{2/3}, \\ \beta_0 + a r^{2/3} \text{Re}_p \text{Ga}^{-4/3}, \end{cases} \qquad (51) $$

depending on the set of the three independent dimensionless numbers chosen. In Eq. (51), $a$ is found to be $a = 0.85$ for the explored values of $\theta$ and $\text{Re}_p$ [dashed line in Fig. 8(c), dotted line corresponds to $a = 1$]. Note in particular that this scaling is not expected to hold at large $\text{Re}_p$ and $\theta$ for which different mechanisms probably become dominant. To understand the obtain scaling law (51), one can discuss its implication at small $\phi$. If $\phi \ll 1$, combining Eqs. (32) and (51) leads to

$$ \begin{aligned} \eta_{\text{eq}}^f / \eta_f &= 1 + \phi [\beta_0 + a (r \theta \sqrt{\text{Re}_p})^{2/3}], \\ &\equiv 1 + \beta_0 \phi [1 + l_g^2 \rho_f \dot{\Gamma} / \eta_f], \\ &= 1 + \beta_0 \phi \left[1 + \left(\frac{l_g}{d}\right)^2 \text{Re}_p\right], \end{aligned} \qquad (52) $$

where $l_g = \sqrt{a/\beta_0} r^{1/3} \text{Ga}^{-2/3} d$ is a length scale. Note that we found $l_g$ in the range $0.4 \lesssim l_g/d \lesssim 3$ in our simulations. Green dashed lines in the insert of 8(a) highlight the relevance of such small-$\phi$ expansion, again for $(\text{Re}_p, \theta, \text{Ga}) = (0.1, 0.22, 0.67)$ and $(\text{Re}_p, \theta, \text{Ga}) = (10, 0.67, 3.9)$. By analogy with the eddy viscosity in turbulent flow, in which $l_g$ would be a mixing length, here it characterizes a viscous length scale associated with the vertical motion of grains induced by their own weight along the vertical direction during transport. This shall therefore somehow be linked to the thickness of the viscous boundary layer developing at the surface of these grains because of relative motion with the fluid along the vertical direction. In other words, this “eddy” viscosity is induced by vertical agitation of the grains moving at the surface of the granular bed and subjected to the characteristic shear flow $\dot{\Gamma}$. This contribution to the fluid-phase dissipation is therefore associated with grain weight. This probably explain why the obtained $\text{Re}_p$ contribution differs from Ref. [48], which found a $\text{Re}_p^{3/2}$ correction in the case of dilute and homogeneous neutrally buoyant suspension in a simple shear flow, i.e., in which inertial contribution comes from other mechanisms.

The value of $\beta_0$ then specifies the viscous limit for the effective viscosity of the equivalent fluid phase in the dilute regime, $\phi \ll 1$. It is observed to be quite smaller than 2.5 which would be the expected value obtained for homogeneous suspension. However, it shall be noted that such viscous limit has usually been obtained for nearly neutrally buoyant suspension, i.e., for $Ga \to 0$ and thus $l_g/d \to \infty$ here. Such limit is thus undefined with the obtained scaling Eqs. (51) and (52), as this would actually be the $\theta \to \infty$ limit, which is not relevant for bedload. Here, the $\phi$-dependency at small $\phi$ is a very different process, and is therefore not comparable to a dilute neutrally buoyant suspension. The dissipation in the fluid phase remains associated with the specific dynamics of heavy grains in close contact, which can modify a purely 1D-vertical diffusive process due to vertical agitation of grains.

Moreover, in the finite range of parameters $(\theta, \text{Ga}, \text{Re}_p)$ covered here for bedload application, the transition from the top fluid layer towards the granular bed remains a sharp interface at the grain scale (as the bed remains relatively dense, even if $\phi$ varies within this bed layer as discussed in the next section). The present spatial upscaling is performed across the entire system and hence it leads to an equivalent diffuse interface as grains remain microscale. The smooth $\phi$-transition observed in Fig. 8(a) is therefore the consequence of this diffuse interface between the bed and the fluid layer at the continuous scale. Such diffuse interface includes the dissipative mechanisms associated with heavy grains mentioned above. $\theta \to 0$ corresponds to the sharpest interface, i.e., for a fluid flow over a porous medium built of solid grains. In the range of $\text{Re}_p$ covered here, and according to (51),
<figure><img src="image_11.png" /><figcaption>FIG. 9. Vertical distribution of the normal components of the equivalent granular phase ($\text{Re}_p = 1, \theta = 0.67$, Ga = 1.22): (solid black) $\Sigma_{yy}^p$; (blue) $\Sigma_{xx}^p$; (red) $\Sigma_{zz}^p$; and (dashed black) $P^p = (\Sigma_{xx}^p + \Sigma_{yy}^p + \Sigma_{zz}^p)/3$.</figcaption></figure>

$\beta \to \beta_0$ when $\theta \to 0$ for nonzero $\text{Re}_p$ flow. In the latter case, the bed surface is a sharp transition between the pure fluid phase and the maximum bed packing, $\phi_m \approx 0.6$ here (see next section). For simplicity, one can assume the flow in the limit $\phi = 0$ to be dominated by the pure fluid (top layer) and to be zero within the porous medium, comparatively. With this in mind and using our present upscaling approach, we can show that, in this case, the effective viscosity of the diffuse interface is simply linear in $\phi$ as $\eta_{eq}^f/\eta_f = 1 + \phi/\phi_m$ (not detailed here). Here, one obtains $\beta_0 < \phi_m^{-1}$ indicating a more complex flow condition at the bed surface, such as slip velocity which would actually lead to $0 < \beta_0 < \phi_m^{-1} < 2.5$.

For now, the interpretation of the model proposed for the fluid-phase effective viscosity is only consistent with the obtained scaling laws mentioned previously. This should deserve specific attention in future works for validation and extension on other configurations and parameter ranges.

## 2. The equivalent granular phase

We now consider the apparent rheology of the equivalent granular phase. We compute the vertical distribution of the inertial number $\mathcal{I}(y) = |\dot{\gamma}^p|d/\sqrt{P^p/\rho_p}$ and the viscous number $\mathcal{J}(y) = |\dot{\gamma}^p|\eta_f/P^p$ to obtain $\mathcal{K}(y) = \mathcal{J} + \alpha\mathcal{I}^2$ (with $\alpha = 0.635$, as Ref. [13]).

It has been checked that the diagonal contribution of the stress is nearly isotropic, i.e., $\Sigma_{xx}^p \approx \Sigma_{yy}^p \approx \Sigma_{zz}^p$ (see Fig. 9 obtained for $\text{Re}_p = 1, \theta = 0.67$, Ga = 1.22). However, the y component $\Sigma_{yy}^p$ (solid black line) is smoother and closer to $P^p$ (dashed black line), while the other contributions show more oscillations in the granular bed. We therefore use $-\Sigma_{yy}^p$ instead of $P^p$.

As done in the previous section, the vertical coordinate $y$ can be eliminated to obtain relationship between state variables, still assuming a single-valued relationship. We first focus on the case ($\text{Re}_p = 1, \theta = 0.67$, Ga = 1.22). Then, for instance, the solid volume fraction $\phi$ is plotted as a function of $\mathcal{K}$ in Fig. 10(a) (dot symbols). $\phi(\mathcal{K})$ is well represented by model (37) up to $\mathcal{K} \simeq 25$ as shown in Fig. 10(a). In particular, the black thin dash line corresponds to $b = 1$ as proposed by Ref. [12] and the red solid line is for $b = 0.7$ which allows to provide a more accurate prediction of $\phi(\mathcal{K})$ over a large range of $\mathcal{K}$. In both cases, $\phi_c$ is found to be $\phi_c = 0.615$. Note that this range of validity is larger than intervals explored by Boyer et al. (2011) [12] and Trulsson et al. (2012) [13] for which $\mathcal{K} \le 0.1$, i.e., for $\phi \ge 0.45$. Figure 10(a) then shows that the validity of Eq. (37) seems to be extended to the range $\phi \ge 0.1$, which therefore includes the upper part of the moving granular layer and more surprisingly most of the diffuse interface at the mesoscale (as discussed in the previous section).
<figure><img src="image_11.png" /><figcaption>FIG. 9. Vertical distribution of the normal components of the equivalent granular phase ($\text{Re}_p = 1, \theta = 0.67$, Ga = 1.22): (solid black) $\Sigma_{yy}^p$; (blue) $\Sigma_{xx}^p$; (red) $\Sigma_{zz}^p$; and (dashed black) $P^p = (\Sigma_{xx}^p + \Sigma_{yy}^p + \Sigma_{zz}^p)/3$.</figcaption></figure>

$\beta \rightarrow \beta_0$ when $\theta \rightarrow 0$ for nonzero $\text{Re}_p$ flow. In the latter case, the bed surface is a sharp transition between the pure fluid phase and the maximum bed packing, $\phi_m \approx 0.6$ here (see next section). For simplicity, one can assume the flow in the limit $\phi = 0$ to be dominated by the pure fluid (top layer) and to be zero within the porous medium, comparatively. With this in mind and using our present upscaling approach, we can show that, in this case, the effective viscosity of the diffuse interface is simply linear in $\phi$ as $\eta_{eq}^f/\eta_f = 1 + \phi/\phi_m$ (not detailed here). Here, one obtains $\beta_0 < \phi_m^{-1}$ indicating a more complex flow condition at the bed surface, such as slip velocity which would actually lead to $0 < \beta_0 < \phi_m^{-1} < 2.5$.

For now, the interpretation of the model proposed for the fluid-phase effective viscosity is only consistent with the obtained scaling laws mentioned previously. This should deserve specific attention in future works for validation and extension on other configurations and parameter ranges.

## 2. The equivalent granular phase

We now consider the apparent rheology of the equivalent granular phase. We compute the vertical distribution of the inertial number $\mathcal{I}(y) = |\dot{\gamma}^p|d/\sqrt{P^p/\rho_p}$ and the viscous number $\mathcal{J}(y) = |\dot{\gamma}^p|\eta_f/P^p$ to obtain $\mathcal{K}(y) = \mathcal{J} + \alpha\mathcal{I}^2$ (with $\alpha = 0.635$, as Ref. [13]).

It has been checked that the diagonal contribution of the stress is nearly isotropic, i.e., $\Sigma_{xx}^p \approx \Sigma_{yy}^p \approx \Sigma_{zz}^p$ (see Fig. 9 obtained for $\text{Re}_p = 1, \theta = 0.67$, Ga = 1.22). However, the y component $\Sigma_{yy}^p$ (solid black line) is smoother and closer to $P^p$ (dashed black line), while the other contributions show more oscillations in the granular bed. We therefore use $-\Sigma_{yy}^p$ instead of $P^p$.

As done in the previous section, the vertical coordinate y can be eliminated to obtain relationship between state variables, still assuming a single-valued relationship. We first focus on the case ($\text{Re}_p = 1, \theta = 0.67$, Ga = 1.22). Then, for instance, the solid volume fraction $\phi$ is plotted as a function of $\mathcal{K}$ in Fig. 10(a) (dot symbols). $\phi(\mathcal{K})$ is well represented by model (37) up to $\mathcal{K} \simeq 25$ as shown in Fig. 10(a). In particular, the black thin dash line corresponds to $b = 1$ as proposed by Ref. [12] and the red solid line is for $b = 0.7$ which allows to provide a more accurate prediction of $\phi(\mathcal{K})$ over a large range of $\mathcal{K}$. In both cases, $\phi_c$ is found to be $\phi_c = 0.615$. Note that this range of validity is larger than intervals explored by Boyer et al. (2011) [12] and Trulsson et al. (2012) [13] for which $\mathcal{K} \le 0.1$, i.e., for $\phi \ge 0.45$. Figure 10(a) then shows that the validity of Eq. (37) seems to be extended to the range $\phi \ge 0.1$, which therefore includes the upper part of the moving granular layer and more surprisingly most of the diffuse interface at the mesoscale (as discussed in the previous section).
<figure><img src="image_12.png" /><figcaption>FIG. 10. Constitutive law of the equivalent granular phase ($Re_p = 1$, $\theta = 0.67$, Ga = 1.22). (a) Solid volume fraction $\phi$ versus $K$, (b) apparent friction coefficient $\mu_{eq}^p$ versus $K$, (c) apparent friction coefficient $\mu_{eq}^p$ versus $\phi$: (●) present IBM-DEM simulation; [red solid lines in panels (a)–(c)] $\phi(K)$, $\mu_{eq}^p(K)$, and $\mu_{eq}^p(\phi)$ given by Eqs. (37), (38), and (39), respectively, with $\phi_c = 0.615$, $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.8$, and $\sqrt{K_0} = 0.15$; (black thin dash line) in panel (a), same as red line but for $b = 1$; (green dashed lines) in panels (b), (c), same as red line but for fitting parameters $\mu_2$ and $\sqrt{K_0}$ as in Fig. 11. Confidence intervals of 95% are computed using the geometric standard deviation of the solid fraction (a) or the apparent friction coefficient (b) for each uncorrelated time sample.</figcaption></figure>

We now compute the equivalent friction coefficient of the granular medium $\mu_{eq}^p$ as

$$ \mu_{eq}^p = \frac{\frac{1}{2}(\Sigma_{yy}^p + \Sigma_{yx}^p)}{-\Sigma_{yy}^p}. \qquad (53) $$

Again, we here used $-\Sigma_{yy}^p$ as an estimation of $P^p$. The effective friction coefficient is shown in Fig. 10(b) for $Re_p = 1$ and $\theta = 0.67$ (dot symbols). The model (38) allows to capture the present numerical results on the range of $K$ covered (red solid line), with $\mu_1 = 0.23$, $\mu_2 = 0.8$ and $\sqrt{K_0} = 0.15$ which are, not surprisingly, slightly different from the values obtained by Ref. [13] in 2D. To finish with, $\mu_{eq}^p(\phi)$ can also be extracted with the same procedure [Fig. 10(c)]. Previous fitting parameters are used to provide model (39), shown as red solid line.

<figure><img src="image_13.png" /><figcaption>FIG. 11. Constitutive laws $\phi(K)$ (a), $\mu_{eq}^p(K)$ (b), and $\mu_{eq}^p(\phi)$ (c) of the equivalent granular phase for all the cases ($Re_p$, $\theta$, Ga) considered in the present work: see Table I for detail of the various symbols and Fig. 10 for legend. Green solid lines are Eqs. (37), (38), and (39), respectively, with $\phi_c = 0.615$, $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.75$, and $\sqrt{K_0} = 0.18$. Note that the cases corresponding to the smaller Shield numbers considered here, namely $\theta \le 0.22$, are not shown because they are too close to the threshold of motion of the granular phase.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Constitutive law of the equivalent granular phase ($Re_p = 1$, $\theta = 0.67$, Ga = 1.22). (a) Solid volume fraction $\phi$ versus $K$, (b) apparent friction coefficient $\mu_{eq}^p$ versus $K$, (c) apparent friction coefficient $\mu_{eq}^p$ versus $\phi$: (●) present IBM-DEM simulation; [red solid lines in panels (a)–(c)] $\phi(K)$, $\mu_{eq}^p(K)$, and $\mu_{eq}^p(\phi)$ given by Eqs. (37), (38), and (39), respectively, with $\phi_c = 0.615$, $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.8$, and $\sqrt{K_0} = 0.15$; (black thin dash line) in panel (a), same as red line but for $b = 1$; (green dashed lines) in panels (b), (c), same as red line but for fitting parameters $\mu_2$ and $\sqrt{K_0}$ as in Fig. 11. Confidence intervals of 95% are computed using the geometric standard deviation of the solid fraction (a) or the apparent friction coefficient (b) for each uncorrelated time sample.</figcaption></figure>

We now compute the equivalent friction coefficient of the granular medium $\mu_{eq}^p$ as

$$ \mu_{eq}^p = \frac{\frac{1}{2}(\Sigma_{yy}^p + \Sigma_{yx}^p)}{-\Sigma_{yy}^p} \quad (53) $$

Again, we here used $-\Sigma_{yy}^p$ as an estimation of $P^p$. The effective friction coefficient is shown in Fig. 10(b) for $Re_p = 1$ and $\theta = 0.67$ (dot symbols). The model (38) allows to capture the present numerical results on the range of $K$ covered (red solid line), with $\mu_1 = 0.23$, $\mu_2 = 0.8$ and $\sqrt{K_0} = 0.15$ which are, not surprisingly, slightly different from the values obtained by Ref. [13] in 2D. To finish with, $\mu_{eq}^p(\phi)$ can also be extracted with the same procedure [Fig. 10(c)]. Previous fitting parameters are used to provide model (39), shown as red solid line.

<figure><img src="image_13.png" /><figcaption>FIG. 11. Constitutive laws $\phi(K)$ (a), $\mu_{eq}^p(K)$ (b), and $\mu_{eq}^p(\phi)$ (c) of the equivalent granular phase for all the cases ($Re_p$, $\theta$, Ga) considered in the present work: see Table I for detail of the various symbols and Fig. 10 for legend. Green solid lines are Eqs. (37), (38), and (39), respectively, with $\phi_c = 0.615$, $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.75$, and $\sqrt{K_0} = 0.18$. Note that the cases corresponding to the smaller Shield numbers considered here, namely $\theta \le 0.22$, are not shown because they are too close to the threshold of motion of the granular phase.</figcaption></figure>
Similar analysis can be obtained for all cases ($Re_p$, $\theta$, Ga) considered in this study. Figure 11 presents these results, similarly to Fig. 10. One can observe that the rheology of the granular phase is not much dependent on $\theta$ and $Re_p$, compared to the rheology of the fluid phase. It is therefore assumed to be captured by a unique model of the form of Eqs. (37), (38), and (39). As best fit, one obtains $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.75$, $\sqrt{K_0} = 0.18$ shown with green solid lines in Figs. 11(a)–11(c). Note that obtained values of $\mu_2$ and $\sqrt{K_0}$ are slightly different from the one obtained for the specific case $Re_p = 1$, $\theta = 0.67$, Ga = 1.22) considered previously. This highlights the small discrepancy from case to case; this new fitting model is also shown in Figs. 10(b) and 10(c) for comparison.

### 3. The equivalent mixture phase

Let us finally consider the equivalent mixture phase. As shown in Eq. (27), the forces acting on this phase correspond to the divergence of the total apparent stress tensor and the weight of the mixture. Hence, the vertical distribution of the forces and stresses acting on this phase are equivalent to the sum of the forces and stresses on the fluid and granular equivalent phases presented in Figs. 3–6, respectively. As done for the apparent fluid and granular phases, we now investigate the rheological behavior of the mixture phase.

First, one considers the apparent viscosity $\eta_{eq}^m$ of the mixture phase. The ‘natural’ way for computing $\eta_{eq}^m$ is to use $\eta_{eq}^m = \Sigma_{xy}^{\text{tot}}/(\partial\langle u_x \rangle^m/\partial y) \equiv \Sigma_{xy}^{\text{tot}}/\dot{\gamma}^m$, with $\langle u \rangle^m$ defined as in Eq. (25). Here, however, and for the same reasons as those mentioned in Sec. IV A 1, we compute the velocity gradient using $\langle u \rangle = \phi\langle u \rangle^p + \epsilon\langle u \rangle^f$. In the present work, $\eta_{eq}^m$ is thus computed as

$$ \eta_{eq}^m = \frac{\Sigma_{xy}^{\text{tot}}}{\dot{\gamma}} = \frac{\frac{1}{2}(\Sigma_{xy}^p + \Sigma_{yx}^p + \Sigma_{xy}^f + \Sigma_{yx}^f)}{\frac{\partial \langle u_x \rangle}{\partial y}}, \quad (54) $$

with $\dot{\gamma}$ defined as in Eq. (50)

From the vertical profiles of both $\eta_{eq}^m(y)/\eta_f$ and $\phi(y)$, we extract the evolution of the apparent viscosity of the equivalent mixture phase $\eta_{eq}^m/\eta_f$ as a function of $\phi$ for $Re_p = 1$, $\theta = 0.67$, Ga = 1.22 (see Fig. 12). For comparison, various models are also plotted in Fig. 12, namely that of Einstein [8,9] (black dotted line), Batchelor and Green (1972) [10] (black dashed line), Krieger and Dougherty (1959) [49] (red dashed line), and Boyer et al. (2011) [12] (blue dashed line).

One can notice, quite surprisingly as will be discussed later on, a good agreement with Einstein's model, Batchelor and Green's model, and Krieger and Dougherty's model up to $\phi \lesssim 0.1$, 0.3, and 0.5, respectively. Recall that the range $0.5 \le \phi \le \phi_m$ corresponds the the static bed region for which the samples used for the time averaging are correlated. Nevertheless, the agreement between the present IBM-DEM results and Krieger and Dougherty's model for $0.5 \le \phi \le \phi_m$ is reasonable. Good agreement is also observed with Boyer et al. (2011)'s model for the whole range of $\phi$. Moreover, the apparent viscosity of the mixture obtained in the fully resolved simulations of Ref. [26] is shown for comparison in the insert of Fig. 12. The three curves are in fair agreement with our simulation at $Re_p = 1$, $\theta = 0.67$, Ga = 1.22 as they fall within the 95% confidence range. It is worth noting a slight increase of the apparent viscosity as the particle Reynolds number is increased in the range $0.1 \lesssim Re_p \lesssim 0.3$ for their simulations.

Figure 13(a) presents the evolution of the apparent friction coefficient for the mixture phase $\mu_{eq}^m$ (defined in Sec. IV B 4) as a function of the viscous number $J$. A general good agreement is observed between the $\mu_{eq}^m$ obtained by experiments and DNS (ours and others). A comparison with model proposed by Ref. [12] shows a relatively good prediction of the total friction. However, the local overestimation of this model around $J \approx 10^{-2}$ indicates over-dissipation due to the choice for the fluid viscosity model. An extension of this model using the present results is discussed in the next section (solid line in Fig. 13). Note that the present DNS results show a deviation of the total friction at high $J$, i.e., close to the pure fluid layer, compared to other studies. The transition toward the pure fluid layer in our specific configuration will be discussed in the next section. However, the
Similar analysis can be obtained for all cases ($Re_p$, $\theta$, Ga) considered in this study. Figure 11 presents these results, similarly to Fig. 10. One can observe that the rheology of the granular phase is not much dependent on $\theta$ and $Re_p$, compared to the rheology of the fluid phase. It is therefore assumed to be captured by a unique model of the form of Eqs. (37), (38), and (39). As best fit, one obtains $b = 0.7$, $\mu_1 = 0.23$, $\mu_2 = 0.75$, $\sqrt{K_0} = 0.18$ shown with green solid lines in Figs. 11(a)–11(c). Note that obtained values of $\mu_2$ and $\sqrt{K_0}$ are slightly different from the one obtained for the specific case $Re_p = 1$, $\theta = 0.67$, Ga = 1.22) considered previously. This highlights the small discrepancy from case to case; this new fitting model is also shown in Figs. 10(b) and 10(c) for comparison.

### 3. The equivalent mixture phase

Let us finally consider the equivalent mixture phase. As shown in Eq. (27), the forces acting on this phase correspond to the divergence of the total apparent stress tensor and the weight of the mixture. Hence, the vertical distribution of the forces and stresses acting on this phase are equivalent to the sum of the forces and stresses on the fluid and granular equivalent phases presented in Figs. 3–6, respectively. As done for the apparent fluid and granular phases, we now investigate the rheological behavior of the mixture phase.

First, one considers the apparent viscosity $\eta_{eq}^m$ of the mixture phase. The ‘natural’ way for computing $\eta_{eq}^m$ is to use $\eta_{eq}^m = \Sigma_{xy}^{\text{tot}}/(\partial\langle u_x \rangle^m/\partial y) \equiv \Sigma_{xy}^{\text{tot}}/\dot{\gamma}^m$, with $\langle u \rangle^m$ defined as in Eq. (25). Here, however, and for the same reasons as those mentioned in Sec. IV A 1, we compute the velocity gradient using $\langle u \rangle = \phi\langle u \rangle^p + \epsilon\langle u \rangle^f$. In the present work, $\eta_{eq}^m$ is thus computed as

$$ \eta_{eq}^m = \frac{\Sigma_{xy}^{\text{tot}}}{\dot{\gamma}} = \frac{\frac{1}{2}(\Sigma_{xy}^p + \Sigma_{yx}^p + \Sigma_{xy}^f + \Sigma_{yx}^f)}{\frac{\partial \langle u_x \rangle}{\partial y}}, \quad (54) $$

with $\dot{\gamma}$ defined as in Eq. (50)

From the vertical profiles of both $\eta_{eq}^m(y)/\eta_f$ and $\phi(y)$, we extract the evolution of the apparent viscosity of the equivalent mixture phase $\eta_{eq}^m/\eta_f$ as a function of $\phi$ for $Re_p = 1$, $\theta = 0.67$, Ga = 1.22 (see Fig. 12). For comparison, various models are also plotted in Fig. 12, namely that of Einstein [8,9] (black dotted line), Batchelor and Green (1972) [10] (black dashed line), Krieger and Dougherty (1959) [49] (red dashed line), and Boyer et al. (2011) [12] (blue dashed line).

One can notice, quite surprisingly as will be discussed later on, a good agreement with Einstein’s model, Batchelor and Green’s model, and Krieger and Dougherty’s model up to $\phi \lesssim 0.1$, 0.3, and 0.5, respectively. Recall that the range $0.5 \le \phi \le \phi_m$ corresponds the the static bed region for which the samples used for the time averaging are correlated. Nevertheless, the agreement between the present IBM-DEM results and Krieger and Dougherty’s model for $0.5 \le \phi \le \phi_m$ is reasonable. Good agreement is also observed with Boyer et al. (2011)’s model for the whole range of $\phi$. Moreover, the apparent viscosity of the mixture obtained in the fully resolved simulations of Ref. [26] is shown for comparison in the insert of Fig. 12. The three curves are in fair agreement with our simulation at $Re_p = 1$, $\theta = 0.67$, Ga = 1.22 as they fall within the 95% confidence range. It is worth noting a slight increase of the apparent viscosity as the particle Reynolds number is increased in the range $0.1 \lesssim Re_p \lesssim 0.3$ for their simulations.

Figure 13(a) presents the evolution of the apparent friction coefficient for the mixture phase $\mu_{eq}^m$ (defined in Sec. IV B 4) as a function of the viscous number $J$. A general good agreement is observed between the $\mu_{eq}^m$ obtained by experiments and DNS (ours and others). A comparison with model proposed by Ref. [12] shows a relatively good prediction of the total friction. However, the local overestimation of this model around $J \approx 10^{-2}$ indicates over-dissipation due to the choice for the fluid viscosity model. An extension of this model using the present results is discussed in the next section (solid line in Fig. 13). Note that the present DNS results show a deviation of the total friction at high $J$, i.e., close to the pure fluid layer, compared to other studies. The transition toward the pure fluid layer in our specific configuration will be discussed in the next section. However, the
<figure><img src="image_14.png" /></figure>

FIG. 12. Apparent viscosity $\eta_{eq}^m/\eta_f$ of the equivalent mixture phase vs solid volume fraction $\phi$: (●) present IBM-DEM simulation for $Re_p = 1$, $\theta = 0.67$, $Ga = 1.22$ [Eq. (54)]; (dotted black) Einstein's model [8,9] $\eta_{eq}^m/\eta_f = 1 + \frac{5}{2}\phi$; (dashed black) Batchelor and Green's model [10] $\eta_{eq}^m/\eta_f = 1 + \frac{5}{2}\phi + 7.6\phi^2$; (dashed red) Krieger and Dougherty's model [49] $\eta_{eq}^m/\eta_f = (1 - \phi/\phi_m)^{-\frac{5}{2}}\phi_m$ with $\phi_m = 0.615$; (dashed blue) Boyer et al.'s model [12] $\eta_{eq}^m/\eta_f = 1 + \frac{5}{2}\phi(1 - \phi/\phi_c)^{-1} + \mu_{eq}^p(\phi)[\phi/(\phi_c - \phi)]^2$ using fittings parameters as obtained for the granular phase (see legend of Fig. 11). Insert: linear representation. Confidence intervals of 95% are computed using the geometric standard deviation of the apparent viscosity for each uncorrelated time sample. Insert: comparison with the resolved simulations of Vowinckel et al. (2021) [26] (colored circles, see the legend of Fig. 13 for key).

Euler-Lagrange simulations slightly underestimate the total friction over the entire range of $\mathcal{J}$. This is probably due to the choice of the viscosity model used in the equivalent fluid phase usually found for these methods. As different models of fluid stress can be found in the literature [20,55] for this method, it would be interesting to deserve a specific attention on that point in future works. All

<figure><img src="image_15.png" /></figure>

FIG. 13. (a) Apparent friction coefficient for the mixture phase $\mu_{eq}^m$ as a function of the viscous number $\mathcal{J}$: a comparison between the experiments of Houssais et al. (2016) [16] here denoted H16 (blue dots), the Euler-Lagrange simulations of Pahtz and Duran (2018) [20] here denoted PD18 (gray triangles), the resolved simulations of Vowinckel et al. (2021) [26] here denoted V21 (gray dots) and the present simulations (black diamonds). The dot line corresponds to the apparent friction coefficient of a suspension flow proposed by Ref. [12], $\mu^{\text{BGP}} = \mathcal{J} + \frac{5}{2}\phi_m\sqrt{\mathcal{J}}$, with $\phi_m = 0.585$. Red line corresponds to the combination of effective friction models for granular and fluid phases, $\mathcal{F}^{-1}\eta_{eq}^f/\eta_f + \mu_{eq}^p$, as developed in the present paper, using Eqs. (32) and (51) to model the equivalent fluid phase, and Eqs. (37) and (38) and fitting parameters from Fig. 10 to model the equivalent granular phase. $\mathcal{F}^{-1}$ is computed using Eq. (57). (b) Constitutive law $\mu_{eq}^m(\mathcal{K})$ of the equivalent mixture phase for all the cases ($Re_p, \theta, Ga$) considered in the present work: see Table I for detail of the various symbols and Fig. 10 for legend.
<figure><img src="image_14.png" /></figure>

FIG. 12. Apparent viscosity $\eta_{eq}^m/\eta_f$ of the equivalent mixture phase vs solid volume fraction $\phi$: (●) present IBM-DEM simulation for $Re_p = 1$, $\theta = 0.67$, $Ga = 1.22$ [Eq. (54)]; (dotted black) Einstein's model [8,9] $\eta_{eq}^m/\eta_f = 1 + \frac{5}{2}\phi$; (dashed black) Batchelor and Green's model [10] $\eta_{eq}^m/\eta_f = 1 + \frac{5}{2}\phi + 7.6\phi^2$; (dashed red) Krieger and Dougherty's model [49] $\eta_{eq}^m/\eta_f = (1 - \phi/\phi_m)^{-\frac{5}{2}}\phi_m$ with $\phi_m = 0.615$; (dashed blue) Boyer et al.'s model [12] $\eta_{eq}^m/\eta_f = 1 + \frac{5}{2}\phi(1 - \phi/\phi_c)^{-1} + \mu_{eq}^p(\phi)[\phi/(\phi_c - \phi)]^2$ using fittings parameters as obtained for the granular phase (see legend of Fig. 11). Insert: linear representation. Confidence intervals of 95% are computed using the geometric standard deviation of the apparent viscosity for each uncorrelated time sample. Insert: comparison with the resolved simulations of Vowinckel et al. (2021) [26] (colored circles, see the legend of Fig. 13 for key).

Euler-Lagrange simulations slightly underestimate the total friction over the entire range of $\mathcal{J}$. This is probably due to the choice of the viscosity model used in the equivalent fluid phase usually found for these methods. As different models of fluid stress can be found in the literature [20,55] for this method, it would be interesting to deserve a specific attention on that point in future works. All

<figure><img src="image_15.png" /></figure>

FIG. 13. (a) Apparent friction coefficient for the mixture phase $\mu_{eq}^m$ as a function of the viscous number $\mathcal{J}$: a comparison between the experiments of Houssais et al. (2016) [16] here denoted H16 (blue dots), the Euler-Lagrange simulations of Pahtz and Duran (2018) [20] here denoted PD18 (gray triangles), the resolved simulations of Vowinckel et al. (2021) [26] here denoted V21 (gray dots) and the present simulations (black diamonds). The dot line corresponds to the apparent friction coefficient of a suspension flow proposed by Ref. [12], $\mu^{\text{BGP}} = \mathcal{J} + \frac{5}{2}\phi_m\sqrt{\mathcal{J}}$, with $\phi_m = 0.585$. Red line corresponds to the combination of effective friction models for granular and fluid phases, $\mathcal{F}^{-1}\eta_{eq}^f/\eta_f + \mu_{eq}^p$, as developed in the present paper, using Eqs. (32) and (51) to model the equivalent fluid phase, and Eqs. (37) and (38) and fitting parameters from Fig. 10 to model the equivalent granular phase. $\mathcal{F}^{-1}$ is computed using Eq. (57). (b) Constitutive law $\mu_{eq}^m(\mathcal{K})$ of the equivalent mixture phase for all the cases ($Re_p, \theta, Ga$) considered in the present work: see Table I for detail of the various symbols and Fig. 10 for legend.
<figure><img src="image_16.png" /></figure>

FIG. 14. (a) Apparent viscosity $\eta_{eq}^m/\eta_f$ of the mixture phase vs solid volume fraction $\phi$ for all the cases considered in the present work: see Table I for detail of the various symbols (solid gray line correspond to Einstein's model (28) at small $\phi$). (b) Values of $\phi_m$ from a best fit of Krieger and Dougherty's model [49] $\eta_{eq}^m/\eta_f = (1 - \phi/\phi_m)^{-5/2}\phi_m$; (▼), Re$_p = 0.1$; (▲), Re$_p = 1$; (▲), Re$_p = 10$; for all $\theta$ with the same markers as in Fig. 8. The dashed line indicates the corresponding mean value: $\phi_m \simeq 0.622$.

this indicates that the contribution of fluid phase and accounting for the transition from the granular bed to the pure fluid layer in mixture models for localized-shear and inertial configurations remain crucial (Sec. IV B 4).

Now, Figs. 13(b) and 14(a) present the apparent friction $\mu_{eq}^m$ of the mixture phase versus the inertial number $K$ and the apparent viscosity $\eta_{eq}^m/\eta_f$ of the mixture phase versus solid volume fraction $\phi$, respectively, for all the cases considered in the present work. Even if all curves are close to collapse to a single curve on a large range of $\phi$ and $K$ in both cases, some deviation is observed. From a granular perspective, i.e., $\mu_{eq}^m(K)$, one notes that deviation occurs when increasing $K$. This could seem in contradiction with the single law $\phi(K)$ reported in Sec. IV B 2. However, this is in line with recent results of rheological laws of granular suspension obtained from rheometer measurements [56]. Moreover, as the granular law induced by contact, $\mu_{eq}^p(K)$, was also found unique in Sec. IV B 2, we state that the deviation observed on the mixture law comes from the equivalent fluid contribution. As reported in Sec. IV B 1, fluid inertia influences the fluid-particle dissipation law of Einstein-type, which then explain the deviation in mixture rheological law from case-to-case reported here. Note that surprisingly, quasi invariance of $\phi(K)$ would suggests that the $\phi$ law is mostly controlled by the contact law between particles and not on details of the fluid flow regime around the particles. At the other end, from a suspension perspective, i.e., $\eta_{eq}^m(\phi)/\eta_f$, the previously mentioned inertial effect is more clearly observed at relatively small $\phi$ when linking to the pure fluid phase $\phi \to 0$ [see Fig. 14(a) and insert]. As a first attempt, the obtained results for the apparent viscosity are fitted by the Krieger and Dougherty's model (30) in which only $\phi_m$ can be adjusted. Figure 14(b) displays the value of $\phi_m$ obtained as a function of $\theta$ and for the different values of Re$_p$ considered. The range of variation of $\phi_m$ is rather small, i.e., $\phi_m = 0.62 \pm 0.02$. This is in line with the variations of $\phi_c$ for the equivalent granular phase. However, we can see in Fig. 14(a) that $\eta_{eq}^m/\eta_f$ may vary significantly depending on $\theta$ and Re$_p$, for $\phi \gtrsim 0.1$ (corresponding to $\eta_{eq}^m/\eta_f \gtrsim 1.25). This suggests that such model can not fully capture the mixture rheology in the present configuration of bedload transport. This is actually clearly depicted in Fig. 16(a) in which the obtained effective viscosity $\eta_{eq}^m/\eta_f$ is compared to the Krieger and Dougherty's model (30) at the corresponding $\phi$. More surprisingly, it is worth mentioning that $\eta_{eq}^m/\eta_f$ seems to converge to the same trend as $\phi \lesssim 0.1$, which is actually close to the Einstein viscosity (28). This observation previously reported in Fig. 12, and generalized in Fig. 14(a) (gray line) suggests that part of the contribution, missing in the equivalent fluid phase as $\eta_{eq}^f/\eta_f$ was shown to deviate from Einstein model [Fig. 8(a)], is now balanced by the granular contact contribution.
<figure><img src="image_16.png" /></figure>

FIG. 14. (a) Apparent viscosity $\eta_{eq}^m/\eta_f$ of the mixture phase vs solid volume fraction $\phi$ for all the cases considered in the present work: see Table I for detail of the various symbols (solid gray line correspond to Einstein's model (28) at small $\phi$). (b) Values of $\phi_m$ from a best fit of Krieger and Dougherty's model [49] $\eta_{eq}^m/\eta_f = (1 - \phi/\phi_m)^{-5/2}\phi_m$; (▼), Re$_p = 0.1$; (▲), Re$_p = 1$; (▲), Re$_p = 10$; for all $\theta$ with the same markers as in Fig. 8. The dashed line indicates the corresponding mean value: $\phi_m \simeq 0.622$.

this indicates that the contribution of fluid phase and accounting for the transition from the granular bed to the pure fluid layer in mixture models for localized-shear and inertial configurations remain crucial (Sec. IV B 4).

Now, Figs. 13(b) and 14(a) present the apparent friction $\mu_{eq}^m$ of the mixture phase versus the inertial number $K$ and the apparent viscosity $\eta_{eq}^m/\eta_f$ of the mixture phase versus solid volume fraction $\phi$, respectively, for all the cases considered in the present work. Even if all curves are close to collapse to a single curve on a large range of $\phi$ and $K$ in both cases, some deviation is observed. From a granular perspective, i.e., $\mu_{eq}^m(K)$, one notes that deviation occurs when increasing $K$. This could seem in contradiction with the single law $\phi(K)$ reported in Sec. IV B 2. However, this is in line with recent results of rheological laws of granular suspension obtained from rheometer measurements [56]. Moreover, as the granular law induced by contact, $\mu_{eq}^p(K)$, was also found unique in Sec. IV B 2, we state that the deviation observed on the mixture law comes from the equivalent fluid contribution. As reported in Sec. IV B 1, fluid inertia influences the fluid-particle dissipation law of Einstein-type, which then explain the deviation in mixture rheological law from case-to-case reported here. Note that surprisingly, quasi invariance of $\phi(K)$ would suggests that the $\phi$ law is mostly controlled by the contact law between particles and not on details of the fluid flow regime around the particles. At the other end, from a suspension perspective, i.e., $\eta_{eq}^m(\phi)/\eta_f$, the previously mentioned inertial effect is more clearly observed at relatively small $\phi$ when linking to the pure fluid phase $\phi \to 0$ [see Fig. 14(a) and insert]. As a first attempt, the obtained results for the apparent viscosity are fitted by the Krieger and Dougherty's model (30) in which only $\phi_m$ can be adjusted. Figure 14(b) displays the value of $\phi_m$ obtained as a function of $\theta$ and for the different values of Re$_p$ considered. The range of variation of $\phi_m$ is rather small, i.e., $\phi_m = 0.62 \pm 0.02$. This is in line with the variations of $\phi_c$ for the equivalent granular phase. However, we can see in Fig. 14(a) that $\eta_{eq}^m/\eta_f$ may vary significantly depending on $\theta$ and Re$_p$, for $\phi \gtrsim 0.1$ (corresponding to $\eta_{eq}^m/\eta_f \gtrsim 1.25). This suggests that such model can not fully capture the mixture rheology in the present configuration of bedload transport. This is actually clearly depicted in Fig. 16(a) in which the obtained effective viscosity $\eta_{eq}^m/\eta_f$ is compared to the Krieger and Dougherty's model (30) at the corresponding $\phi$. More surprisingly, it is worth mentioning that $\eta_{eq}^m/\eta_f$ seems to converge to the same trend as $\phi \lesssim 0.1$, which is actually close to the Einstein viscosity (28). This observation previously reported in Fig. 12, and generalized in Fig. 14(a) (gray line) suggests that part of the contribution, missing in the equivalent fluid phase as $\eta_{eq}^f/\eta_f$ was shown to deviate from Einstein model [Fig. 8(a)], is now balanced by the granular contact contribution.
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$.