# A modified kinetic theory for frictional-collisional bedload transport valid
from dense to dilute regime.

Rémi Chassagne<EMAIL_ADDRESS>Univ. Grenoble Alpes,
LEGI, CNRS UMR 5519 - Grenoble, France  Julien Chauchat  Univ. Grenoble Alpes,
LEGI, CNRS UMR 5519 - Grenoble, France  Cyrille Bonamy  Univ. Grenoble Alpes,
LEGI, CNRS UMR 5519 - Grenoble, France

(August 29, 2024)

###### Abstract

Modelling sediment transport is still a challenging problem and is of major
importance for the study of particulate geophysical flows. In this work, the
modelling of sediment transport in the collisional regime is investigated with
a focus on the continuum modelling of the granular flow. For this purpose, a
frictional-collisional approach, combining a Coulomb model with the kinetic
theory of granular gases, is developed. The methodology is based on a
comparison with coupled fluid-discrete simulations, that classical kinetic
theory model fails to reproduce. In order to improve the continuum model, the
fluctuating energy balance is computed in the discrete simulations and
systematically compared with the kinetic theory closure laws. Inter-particle
friction is shown to affect the radial distribution function and to increase
the energy dissipation, in agreement with previous observations in the
literature. Due to saltating particles, whose motion can not be captured by
the kinetic theory, departure from the viscosity and diffusivity laws are
observed in the dilute part of the granular flow. Finally, the quadratic
nature of the drag force is shown to increase the granular fluctuating energy
dissipation. Based on these observations, modifications of the kinetic theory
closure laws are proposed. The modified model reproduces perfectly the
discrete simulations in the entire depth structure of the granular flow. These
modifications are shown to impact the rheological properties of the flow and
they make it possible to recover the $\mu(I)$ rheology in the dense regime.

##  I Introduction

Sediment transport is one of the main processes that shapes river beds and
coastal regions and is of major importance for their management and
interactions with human infrastructures. Modelling the transport of particles
by a flowing fluid is therefore a great challenge for the prediction of the
morphological evolution of our environment. According to Berzi and
Fraccarollo, [3], sediment transport regimes may be differentiated using the
bed slope $\alpha$ and the Shields number $\theta=\tau_{b}/(\Delta\rho gd)$,
where $\tau_{b}$ is the fluid bed shear stress, $\Delta\rho$ is the density
difference between particles and fluid, $g$ is the gravity acceleration and
$d$ is the particle diameter. At bed slopes above approximately $5^{\circ}$
debris flows are observed for all Shields numbers. Below this value,
increasing the Shields number from just above the critical value, sediment
transport takes place as ordinary bedload corresponding to a single layer of
grains moving on top of the fixed sediment bed. For $\theta\geq 0.2$ a
transition is observed from collisional sediment transport to collisional
transport with turbulent suspension at larger Shields numbers. While ordinary
bedload may be modelled using statistical approaches (e.g. [22, 1])
collisional transport and turbulent suspension are usually modelled in the
continuum framework. The challenge in these approaches lies in the complexity
of the particle-particle and turbulence-particle interactions modelling. Over
the last two decades, extensive research has been carried out on this problem.
Most of the time, Reynolds-averaged turbulence modelling have been used such
as mixing length model (e.g. [36, 54]) or two-equations models such as
$k-\epsilon$ or $k-\omega$ (e.g. [30, 41, 12, 28]). More recently, Large Eddy
Simulation has been used to model turbulence in these two-phase flows (e.g.
[14, 46]). Concerning the granular stress modelling, two approaches have been
used, $\mu(I)$ rheology (e.g. [54, 41, 12]) and Kinetic Theory of granular
flows (e.g. [36, 31, 4, 14]). This work investigates the collisional transport
regime with a particular focus on the modelling of the granular part of the
flow.

In order to model granular flows continuously, two approaches are generally
considered depending on the granular flow regime. Dense dry granular flows can
be described with the phenomenological $\mu(I)$ rheology [26, 17, 51], which
relates the stress state, represented by the shear to normal stress ratio or
effective friction coefficient $\mu=\tau^{p}/P^{p}$, to the kinetic state of
the granular flow, represented by the inertial number
$I=d|\dot{\gamma}|/\sqrt{P^{p}/\rho^{p}}$, where $d$ is the particle diameter,
$\rho^{p}$ the particle density and $\dot{\gamma}$ is the granular velocity
shear rate. This phenomenological law, derived by fitting data, gives
predictive results for dense granular flows but it fails in the dilute regime.
In the bedload configuration, Maurin et al., [49] studied the rheology of
monodisperse beds with coupled fluid-Discrete Element Method (DEM). Despite
the presence of water, they showed that the dry inertial number is still the
controlling parameter. They found the $\mu(I)$ rheology to be valid in the
dense part of the granular flow and extended the classical laws for dry
granular flows to the bedload case as follows

$\mu(I)=\mu_{s}+\dfrac{\mu_{2}-\mu_{s}}{I_{0}/I+1},$ (1)
$\phi(I)=\dfrac{\phi_{c}}{1+bI},$ (2)

with $\mu_{s}=0.35$ is the static friction coefficient below which no motion
is a priori possible, $\mu_{2}=0.97$, $I_{0}=0.69$, $\phi_{c}=0.61$ and
$b=0.31$. These laws persist for inertial number of the order of unity but
depart from data at higher inertial numbers in the dilute regime.

Dilute granular flows are generally described with the Kinetic Theory (KT) of
granular gases, which is based on a statistical description of the granular
flow and binary collisions between particles [9]. It relies on the Enskog-
Boltzmann equation for the velocity distribution function and differs from the
classical KT of molecular gases by taking into account the finite size of
particles and energy dissipation at contacts due to inelasticity. It provides
hydrodynamics conservation equations for mass, momentum and granular
temperature. The granular temperature is defined as the averaged granular
fluctuating kinetic energy and is not at all related to thermal temperature.
These equations need closures for the stress tensor, granular temperature
diffusivity and rate of dissipation. In the framework of the KT they can be
expressed as complex integrals of the particle distribution function (solution
of the Enskog-Boltzmann equation) and of the radial distribution function
$g_{0}$, which characterizes the degree of spatial correlation between two
particles that are going to collide. Under certain assumptions, analytical
expressions for the coefficients may be obtained by application of the
Chapmann-Enskog method [9], leading to different models. The main difficulty
is to express the particle distribution function. The simplest model were not
solving the Enskog-Boltzmann equation and assumed to be at equilibrium (no
temporal nor spatial gradients) leading to the classical Maxwell distribution
function [56, 37]. Lun et al., [44] considered a small perturbation around the
Maxwell distribution and obtained constitutive relations for the coefficients
that are valid for weakly dissipative granular flows, i.e. $1-e^{2}<<1$ with
$e$ the restitution coefficient. More recently, performing expansion to first
order in spatial gradients, Garzo and Dufty [25] derived relations for the
whole range of restitution coefficient. Considering the large spatial
variations of all quantities in the bedload configuration, i.e. exponential
velocity profile, this latter model is certainly the most relevant one for the
present study.

Considering bedload, and environmental flows over an erodible bed in general,
the difficulty arises from the coexistence of all granular flow regimes with
transition from a dilute granular flow at the bed surface to dense and quasi-
static regime inside the bed. In terms of continuum modeling strategy, two
options are possible: extend the $\mu(I)$ rheology to the dilute regime or
extend the KT to the dense regime. In the dilute regime, non-local effects are
important features of the flow and it seems therefore hardly feasible to model
all regimes with a local rheology. Consequently the second approach is adopted
hereafter.

There are various reasons to observe departure from the KT even in the dilute
part of the flow. First because natural particles are frictional while the
classical KT models such as the one of Garzó and Dufty, [25] have been derived
for frictionless particles. Interparticle friction makes it possible to
transfer translational kinetic energy into rotational kinetic energy during
collisions and add another source of energy dissipation by sliding motion.
There have been some attempts to generalize KT to frictional particles
including the resolution of rotational kinetic energy and a tangential
restitution coefficient [53]. This however makes models even more complex and
some authors proposed a more pragmatic approach in which the restitution
coefficient is empirically reduced to account for interparticle friction
dissipation [39, 16].

At the top of the sediment bed particles are transported by saltation with
ballistic trajectories. Their motion is controlled by gravity and fluid drag
force and is therefore out of the scope of the KT [6]. Therefore, the KT is
not expected to reproduce bedload in the very dilute part of the flow. There
have been few attemps to model saltation in the continuous framework [38] and
this is still a challenge for the modelling of geophysical particulate flows.

Finally, the KT is not expected to predict the behavior of the granular flow
in the dense regime. This is because the motion of particles starts to be
correlated and the molecular chaos assumption breaks down while particles are
more concentrated. To model dense flows, Jenkins, [34, 35] and Berzi and
Jenkins, [4] proposed the introduction of a correlation length in the
dissipation term of the temperature equation, showing good results when
comparing with DEM shear flow simulations [2, 5]. In addition to these
considerations, the origin of stresses changes when going into the bed. While
stresses are due to fluctuating motions and very short collisions in the
dilute regime, force chains start to emerge with long lasting contacts in the
dense part of the bed leading to elastic stresses that can not be captured by
the KT. Johnson and Jackson, [40] proposed that the elastic stresses follow a
Coulomb-like behavior and assumed that they can simply be summed up with the
kinetic stresses to compute the total stresses. With this approach, Berzi et
al., [6] were able to reproduce DEM simulations of dry granular flows on an
erodible bed.

The literature review presented above highlights the duality between the
$\mu(I)$ rheology, valid in the dense regime, and the KT a priori valid in the
dilute regime. A two-fluid model based on the KT that would be able to
reproduce quantitatively the entire depth structure of the bedload layer would
represent a major contribution for the sediment transport community. In
particular, reproducing the $\mu(I)$ rheology with a KT based model is still
an open question. In order to focus on the granular flow modelling and on the
particle-particle interactions, collisional transport only is considered
without any turbulent suspension and fluid-particle interactions will be
reduced to the simplest ones, i.e. buoyancy and drag forces.

In this context, the aim of the present work is to propose a two-fluid model
for bedload transport using the KT for granular flows and based on the
frictional-collisional approach first proposed by Johnson and Jackson, [40].
The model is implemented in SedFoam [12], which is an open-source three
dimensional (3D) solver based on the OpenFOAM toolbox. Following the same
approach as Chialvo and Sundaresan, [16] in simple shear cell granular flows,
comparison with coupled fluid-DEM simulations will be performed. Keeping in
mind the difficulties highlighted previously, empirical modifications to the
classical KT model of Garzó and Dufty, [25] will be proposed in order to
reproduce quantitatively the discrete simulations.

First the two phase flow model (section II) and the coupled fluid-DEM model
(section II.2) are presented. The results of both models will be compared and
some modifications to the Garzó and Dufty, [25] model will be proposed in
section III and implemented in SedFoam to quantitatively reproduce the
discrete simulations. Finally, results are discussed in section IV.

##  II Two-phase flow models for bedload transport

###  II.1 Two-fluid model

The 3D two-fluid model for sediment transport described in Chauchat et al.,
[12], based on the equations derived by Jackson, [33], is used in this
contribution. Although it has been written in 3D, the model will only be used
in 1D, with all variables depending only on the vertical axis $z$. Herein,
only the simplified 1D equations will be presented [11], and the interested
reader is referred to Chauchat et al., [12] for more details on the 3D model.
The release of the code used in this work can be found in Bonamy et al., [7].
In the following, overscript f (resp. p) denotes for fluid (resp. particle)
phase quantity. The mass and momentum conservation equations aims at solving
the fluid and solid phase volume fraction, respectively $\epsilon$ and
$\phi=1-\epsilon$, as well as the Favre-averaged fluid and particle phase
velocities, respectively $\bm{u}^{f}=(u_{x}^{f},u_{z}^{f})$ and
$\bm{v}^{p}=(v_{x}^{p},v_{z}^{p})$. The Favre-averaged operators are detailed
in appendix A.

####  II.1.1 Mass and momentum conservation equations

The mass conservation equations of the fluid and particle phase reduce to

$\dfrac{\partial\rho^{f}\epsilon}{\partial t}+\dfrac{\partial\rho^{f}\epsilon
u_{z}^{f}}{\partial z}=0,$ (3) $\dfrac{\partial\rho^{p}\phi}{\partial
t}+\dfrac{\partial\rho^{p}\phi v_{z}^{p}}{\partial z}=0,$ (4)

where $\rho^{f}$ and $\rho^{p}$ are respectively the fluid and particle
density. The momentum conservation equations in the streamwise direction are
written as

$\dfrac{\partial\epsilon\rho^{f}u_{x}^{f}}{\partial t}+\dfrac{\partial
u_{z}^{f}\epsilon\rho^{f}u_{x}^{f}}{\partial
z}=\epsilon\rho^{f}g\sin(\alpha)+\dfrac{\partial\tau_{xz}^{f}}{\partial
z}-nf_{Dx},$ (5) $\dfrac{\partial\phi\rho^{p}v_{x}^{p}}{\partial
t}+\dfrac{\partial v_{z}^{p}\phi\rho^{p}v_{x}^{p}}{\partial
z}=\phi\rho^{p}g\sin(\alpha)+\dfrac{\partial\tau_{xz}^{p}}{\partial
z}+nf_{Dx},$ (6)

and in the wall-normal direction

$\dfrac{\partial\epsilon\rho^{f}u_{z}^{f}}{\partial t}+\dfrac{\partial
u_{z}^{f}\epsilon\rho^{f}u_{z}^{f}}{\partial
z}=-\epsilon\rho^{f}g\cos(\alpha)-\epsilon\dfrac{\partial p^{f}}{\partial
z}-nf_{Dz},$ (7) $\dfrac{\partial\phi\rho^{p}v_{z}^{p}}{\partial
t}+\dfrac{\partial v_{z}^{p}\phi\rho^{p}v_{z}^{p}}{\partial
z}=-\phi\rho^{p}g\cos(\alpha)-\dfrac{\partial p^{p}}{\partial
z}-\phi\dfrac{\partial p^{f}}{\partial z}+nf_{Dz},$ (8)

where $p^{f}$ (resp. $p^{p}$) is the fluid (resp. particle) pressure,
$\tau_{xz}^{f}$ (resp. $\tau_{xz}^{p}$) is the fluid (resp. particle) shear
stress, $g=9.81$ $m^{2}.s^{-1}$ is the gravity acceleration, $\alpha$ is the
bed slope angle and $n\bm{f_{D}}=n(f_{Dx},f_{Dz})$ is the drag force between
the fluid and particle phase.

####  II.1.2 Fluid phase closures

To solve the fluid equations (3), (5) and (7), it is necessary to have
closures for the fluid shear stress $\tau_{xz}^{f}$ and drag force
$n\bm{f_{D}}$. The fluid shear stress is the sum of a viscous shear stress and
turbulent shear stress. The latter stress, also called Reynolds shear stress,
is computed based on an eddy viscosity concept and the total shear stress
expresses

$\tau_{xz}^{f}=\rho^{f}\epsilon(\nu^{f}+\nu_{t})\dfrac{\partial
u_{x}^{f}}{\partial z},$ (9)

where $\nu^{f}$ is the fluid kinematic viscosity. The turbulent viscosity
$\nu_{t}$ follows a mixing length approach that depends on the integral of the
solid volume fraction to account for the presence of the particles [42]

$\nu_{t}=l_{m}^{2}|\dfrac{\partial u_{x}^{f}}{\partial z}|,\quad
l_{m}(z)=\kappa\int_{0}^{z}1-\dfrac{\phi(\xi)}{\phi_{max}}d\xi,$ (10)

with $\kappa=0.41$ the Von-Karman constant. This simple formulation of mixing
length extends the Prandtl, [52] law for flow inside and over a mobile
sediment bed, and has been widely used for sheet-flow and bedload applications
[42, 21, 54, 48].

The drag force is computed as

$n\bm{f_{D}}=\phi(1-\phi)K\left(\bm{u}^{f}-\bm{u}^{p}\right),\quad
K=0.75C_{D}\dfrac{\rho^{f}}{d}||\bm{u}^{f}-\bm{u}^{p}||\left(1-\phi\right)^{-\zeta-1}.$
(11)

The $\left(1-\phi\right)^{-\zeta-1}$ term accounts for hindrance effects due
to the collective presence of particles and the exponent is fixed to
$\zeta=3.1$ [36]. The drag coefficient $C_{D}$ is computed following Dalla
Valle, [18]

$C_{D}=C_{D}^{\infty}+\dfrac{24.4}{Re_{p}},\quad
Re_{p}=\dfrac{||\bm{u}^{f}-\bm{u}^{p}||d}{\nu^{f}},$ (12)

where $C_{D}^{\infty}=0.4$ is the value of the drag coefficient in the limit
of infinite particle Reynolds number.

####  II.1.3 Particle phase closures

As already pointed out in the introduction, granular stresses have two
physical origins. First, particles endure elastic stresses, resulting from
enduring contacts, which do not depend on shear-rate. These elastic stresses
appear at high volume fractions and lead to a frictional Coulomb-like behavior
of the granular bed. Secondly, particles also endure kinetic stresses
resulting from transfer of momentum due to fluctuating motions and at
collisions between particles. In order to account for both these stresses in
our model, the approach of Johnson and Jackson, [40], which consists in adding
both contributions to the total granular stress, is adopted,

$\displaystyle\tau_{xz}^{p}$ $\displaystyle=\tau_{xz}^{el}+\tau_{xz}^{kin},$
(13) $\displaystyle p^{p}$ $\displaystyle=p^{el}+p^{kin}.$ (14)

##### Elastic stresses

Models for elastic stresses originate from soil mechanics or poro-elasticity
and are, for most of them, empirical. In this paper, the model proposed by
Johnson and Jackson, [40] is used,

$p^{el}=\left\\{\begin{array}[]{lr}0,&\phi<\phi_{rlp},\\\
P_{0}\dfrac{(\phi-\phi_{rlp})^{3}}{(\phi_{max}-\phi)^{5}},&\text{otherwise}\end{array}\right.$
(15)

with $P_{0}=0.05$ ($kg.m^{-1}.s^{-2}$), $\phi_{rlp}=0.57$ the random loose
packing fraction. The shear stress is then computed following a Coulomb model,

$\tau^{el}=\mu_{s}p^{el},$ (16)

where $\mu_{s}=0.35$ is the effective static friction coefficient computed by
Maurin et al., [49] with discrete element simulations of bedload transport.
Note that the elastic pressure only exists for volume fraction larger than
$\phi_{rlp}$ and the elastic stresses are therefore restricted to the very
dense part of the granular flow.

##### Kinetic stresses

As already mentioned, the kinetic stresses are computed in the framework of
the KT. One needs to solve an energy balance equation for the particle phase
in addition to the set of equations (4), (6) and (8). For multiphase flows,
the rate of change of particle fluctuating kinetic energy can be written as
[20],

$\dfrac{3}{2}\left(\dfrac{\partial\rho^{p}\phi T}{\partial t}+\dfrac{\partial
u_{z}^{p}\rho^{p}\phi T}{\partial z}\right)=\tau_{xz}^{kin}\dfrac{\partial
v_{x}^{p}}{\partial z}+\dfrac{\partial q}{\partial z}+\Gamma+J_{int},$ (17)

where $T=1/3\left<v_{i}^{p\prime}v_{i}^{p\prime}\right>^{s}$ is the granular
temperature and $\left<\cdot\right>^{s}$ denotes for the particle phase
average. The first term in the Right Hand Side (RHS) is a production term of
temperature by the work of the mean granular flow. The second one in the RHS
is a diffusion term with $q$ the granular temperature flux,

$q=-\kappa\dfrac{\partial T}{\partial z}.$ (18)

The third term in the RHS is dissipation of granular temperature due to
contact inelasticity. The last term is specific to multiphase flows and
represents the work done by the drag force due to fluctuating motion. In order
to close the granular temperature equation, the Garzó and Dufty, [25] model is
used. It has been derived for inelastic hard spheres to first order in the
spatial gradients. The transport coefficients are expressed as a function of
the granular temperature as

$\displaystyle p^{kin}$ $\displaystyle=\rho^{p}F_{1}(\phi)T,$ (19)
$\displaystyle\eta^{kin}$ $\displaystyle=\rho^{p}dF_{2}(\phi)\sqrt{T},$ (20)
$\displaystyle\kappa$ $\displaystyle=\rho^{p}dF_{3}(\phi)\sqrt{T},$ (21)
$\displaystyle\Gamma$ $\displaystyle=\rho^{p}/dF_{4}(\phi)T^{3/2},$ (22)

where $\eta^{kin}=\tau_{xz}^{kin}/|\dot{\gamma}|$ is the viscosity associated
with the kinetic stress. The expressions of $F_{1}$, $F_{2}$, $F_{3}$ and
$F_{4}$ are reported in the left column of table 1 and depend on the
restitution coefficient $e$. Note that the dimensionless viscosity and
diffusivity are expressed as the sum of a kinetic (due to fluctuations)
collisional and bulk contributions.

$\begin{array}[]{rcl|rcl}\lx@intercol\hfil\text{\cite[cite]{[\@@bibref{Number}{garzo1999}{}{}]}
model}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Proposed
model}\hfil\lx@intercol\\\ \hline\cr F_{1}(\phi)&=&\phi\left(1+2(1+e)\phi
g_{0}(\phi)\right)&\lx@intercol\hfil\text{--}\hfil\lx@intercol\\\
F_{2}(\phi)&=&\dfrac{5\sqrt{\pi}}{96}\eta^{*}=\dfrac{5\sqrt{\pi}}{96}\left(\eta_{k}^{*}+\eta_{c}^{*}+\eta_{b}^{*}\right)&\lx@intercol\hfil\text{--}\hfil\lx@intercol\\\
F_{3}(\phi)&=&\dfrac{225\sqrt{\pi}}{1152}\kappa*=\dfrac{225\sqrt{\pi}}{1152}\left(\kappa_{k}^{*}+\kappa_{c}^{*}+\kappa_{b}^{*}\right)&\lx@intercol\hfil\text{--}\hfil\lx@intercol\\\
F_{4}(\phi)&=&\dfrac{12}{\sqrt{\pi}}(1-e^{2})\phi^{2}g_{0}(\phi)&F_{4}^{\prime}(\phi)&=&\dfrac{12}{\sqrt{\pi}}(1-e_{eff}^{2})\phi^{2}g_{0}(\phi)\\\
\eta_{k}^{*}&=&\dfrac{1-2/5(1+e)(1-3e)\phi
g_{0}(\phi)}{\left(1-1/4(1-e)^{2}-5/24(1-e^{2})\right)g_{0}(\phi)}&\eta_{k}^{*}&=&\dfrac{48/(5\sqrt{\pi})\phi-2/5(1+e)(1-3e)\phi
g_{0}(\phi)}{\left(1-1/4(1-e)^{2}-5/24(1-e^{2})\right)g_{0}(\phi)}\\\
\eta_{c}^{*}&=&\dfrac{4}{5}(1+e)\phi
g_{0}(\phi)\eta_{k}^{*}&\lx@intercol\hfil\text{--}\hfil\lx@intercol\\\
\eta_{b}^{*}&=&\dfrac{384}{25\pi}(1+e)\phi^{2}g_{0}(\phi)&\lx@intercol\hfil\text{--}\hfil\lx@intercol\\\
\kappa_{k}^{*}&=&\dfrac{2\left[1+3/5(1+e)^{2}(2e-1)\phi
g_{0}(\phi)\right]}{(1-7/16(1-e))(1+e)g_{0}(\phi)}&\kappa_{k}^{*}&=&\dfrac{2\left[576/(225\sqrt{\pi})\phi+3/5(1+e)^{2}(2e-1)\phi
g_{0}(\phi)\right]}{(1-7/16(1-e))(1+e)g_{0}(\phi)}\\\
\kappa_{c}^{*}&=&\dfrac{6}{5}(1+e)\phi
g_{0}(\phi)\kappa_{k}^{*}&\lx@intercol\hfil\text{--}\hfil\lx@intercol\\\
\kappa_{b}^{*}&=&\dfrac{2304}{225\pi}(1+e)\phi^{2}g_{0}(\phi)&\lx@intercol\hfil\text{--}\hfil\lx@intercol\end{array}$

Table 1: Expression of the coefficients in the constitutive relations of the
kinetic theory given by the [25] model (left column) and the proposed
modifications (right column).

In these relations, a key parameter is the radial distribution function
$g_{0}(\phi)$. It characterizes the degree of spatial correlation of two
particles that are going to collide. It is therefore expected that $g_{0}$ is
one in the dilute regime (no spatial correlation) and diverges when $\phi$
approaches $\phi_{max}$ (full spatial correlations). The first functional form
has been proposed by Carnahan and Starling, [8], as
$g_{0}(\phi)=(2-\phi)/(2(1-\phi)^{3})$. It gives very good results for low
values of packing fraction but does not diverge in the dense limit. To correct
this behavior, other functional forms have been proposed theoretically [45,
59] and empirically [43, 16, 60]. In the following, the form of Chialvo and
Sundaresan, [16] proposed based on discrete element simulations of simple
shear flows, is adopted,

$g_{0}(\phi)=\dfrac{2-\phi}{2(1-\phi)^{3}}+\dfrac{a\phi^{2}}{(\phi_{max}-\phi)^{3/2}},$
(23)

where $a=0.58$ is an empirical coefficient. This expression has the advantage
of matching the Carnahan and Starling, [8] radial distribution function for
low packing fraction and to diverge close to $\phi_{max}$.

The last term for which a closure is required is $J_{int}$. It corresponds to
the work done by the drag force due to fluctuating particle motion. Assuming
that the drag is linear in relative velocity (i.e. $K\sim Const$), it can be
expressed as [20, 23]

$J_{int}=\phi(1-\phi)K(\left<u_{i}^{f\prime}v_{i}^{p\prime}\right>^{s}-3T).$
(24)

The second term in the parenthesis of equation (24) is a dissipation term. The
first term is a source term representing a transfer of fluctuating energy from
the fluid to particles through turbulence. It depends on the degree of
correlation between fluctuating velocities of both phases and typically varies
from $0$ (uncorrelated motions) to $2k$ (fully correlated motions), where $k$
is the fluid fluctuating kinetic energy [19, 13]. In the DEM model that is
used in this study (see section II.2), no fluid velocity fluctuations are
computed and there is therefore no correlation between the fluid and granular
phase fluctuation velocities. This is equivalent to consider that particles
are very inertial and that their motion is not at all influenced by the
turbulent structures of the fluid flow. The drag term therefore simplifies as

$J_{int}=-\phi(1-\phi)K\times 3T,$ (25)

and it represents the dissipation of granular temperature by the drag force.

###  II.2 Euler-Lagrange model

The continuum two-fluid model will be tested against coupled fluid-DEM
simulations. This is a three-dimensional DEM model, based on the open-source
code YADE [58], coupled with a one-dimensional turbulent fluid flow model. It
is presented in [48, 47] and has been validated with experiments [24]. It will
be briefly presented herein but the interested reader is referred to Maurin et
al., [48] and Maurin, [47] for a complete description.

Frictional spheres of density $\rho^{p}$ and diameter $d$ submitted to a
gravity acceleration $g=9.81$ $m.s^{-2}$ are considered. A linear spring-
dashpot model [57] composed of a spring of stiffness $k_{n}$, computed to stay
in the rigid limit of grains [55], in parallel with a viscous damper of
coefficient $c_{n}$ in the normal direction; and a spring of stiffness
$k_{s}=k_{n}$ associated with a slider of friction coefficient $\mu^{p}$ in
the tangential direction. The normal stiffness, in parallel with the viscous
damper, defines a normal restitution coefficient $e$. In addition, each
particle is submitted to a buoyancy and a fluid drag force. The same drag
model as for the two-fluid model is considered (eq.(11) and (12)).

For the fluid phase, a one-dimensional (1D) vertical model is used, in which
the fluid velocity is only a function of the wall-normal component, $z$, and
is aligned with the streamwise direction. This is the exact same fluid model
as presented previously for the two phase flow model (eq. (5)). The only
difference lies in the drag interaction term. Indeed as a fluid drag force
applies to each individual particle, it is necessary to compute the average
momentum transferred by the fluid to the particles through the drag force,

$n\left<\bm{f_{D}}\right>^{s}=\dfrac{3\rho^{f}}{4d}\phi\left(1-\phi\right)^{-\zeta}\left<C_{D}||\bm{u}^{f}-\bm{u}^{p}||(\bm{u}^{f}-\bm{u}^{p})\right>^{s},$
(26)

where $\left<\cdot\right>^{s}$ denotes the solid phase average. The averaging
procedure proposed by Jackson, [32, 33] is adopted. It is based on a weighting
function, representing the volume in which the averaging is performed.
Considering the symmetry of the problem, a cuboid weighting function
$\mathcal{H}$ with the same length and width as the 3D domain is applied. In
the vertical direction, in order to capture the strong gradient of the mean
flow, the vertical thickness of the box is taken as $d_{z}=d/30$. The volume
fraction is computed as

$\phi(z)=\sum_{p}\int_{V_{p}}\mathcal{H}(|z-z^{\prime}|)dV,$ (27)

where $V_{p}$ represents the volume of particle $p$. The average over the
solid phase of any scalar quantity $\gamma$ is computed as

$\left<\gamma\right>^{s}(z)=\dfrac{1}{\phi(z)}\sum_{p}\int_{V_{p}}\gamma(z^{\prime})\mathcal{H}(|z-z^{\prime}|)dV.$
(28)

This averaging procedure makes it possible to compute the total transfer of
momentum between fluid and particles through the drag force and to obtain
vertical profile of all particle quantity (i.e. particle velocity profile).

Goldhirsch, [27] and later Pähtz et al., [50] proposed a rigorous procedure to
compute the DEM granular stresses $\sigma^{p}$, granular temperature flux $q$,
dissipation during collision $\Gamma$ and drag dissipation $J_{int}$. For the
same weighting function as previously, they are expressed as

$\displaystyle\sigma^{p}_{ij}(z)$
$\displaystyle=-\rho^{p}\phi\left<v_{i}^{p\prime}v_{j}^{p\prime}\right>^{s}-\dfrac{1}{V}\sum_{c}f_{i}^{c}b_{j}^{c},$
(29) $\displaystyle q(z)$
$\displaystyle=-\dfrac{1}{2}\rho^{p}\phi\left<v_{z}^{p\prime}v_{i}^{p\prime}v_{i}^{p\prime}\right>^{s}-\dfrac{1}{V}\sum_{c}f_{i}^{c}v_{i}^{p}b_{z}^{c},$
(30) $\displaystyle\Gamma(z)$
$\displaystyle=-\dfrac{1}{2V}\sum_{p,q|z^{p}\in[z-d_{z}/2,z+d_{z}/2]}f_{i}^{c}(v_{i}^{p}-v_{i}^{q}),$
(31) $\displaystyle J_{int}$
$\displaystyle=n\left<f_{Di}v_{i}^{p\prime}\right>^{s},$ (32)

where the sums are performed over the ensemble of particles $p$, particles $q$
in contact with $p$ and over contacts $c$. $\bm{f^{c}}$ is the interaction
force at contact $c$ applied on particle $p$ by particle $q$ and
$\bm{b}^{c}=\bm{x}^{q}-\bm{x}^{p}$ is the branch vector. In these expressions,
the Einstein summation has been used. Note that the granular stress tensor and
granular temperature flux are the sum of kinetic (similar to a Reynolds stress
for a fluid) and contact contributions.

##  III Results

In this section, coupled fluid-discrete simulations will be used to compare
with two-fluid simulations, showing important departures between both models.
The origins of the discrepancies will be physically investigated and
corrections to the KT model will be proposed in order to reproduce
quantitatively the DEM simulations.

###  III.1 Euler-Lagrange simulation results

The numerical setup is presented in figure 1a. Particles of diameter $d=6$ mm
and density $\rho^{p}=2500$ kg.m-3 are initially deposited by gravity over a
rough bed made of fixed particles. Particles are immersed in a fluid of
density $\rho^{f}=1000$ kg.m-3. The free surface is fixed to
$H_{f}=H_{bed}+h_{w}$, where $H_{bed}=12.5d$ represents the bed height at rest
and $h_{w}$ represents the water depth. The slope is fixed to
$\sin(\alpha)=5\%$ ($2.85^{\circ}$). This defines the Shields number as
$\theta=\rho^{f}gh_{w}\sin(\alpha)/\left[(\rho^{p}-\rho^{f})gd\right]$. At
initial time, the fluid flows by gravity and sets particles into motion. After
a short transient, during which fluid and particles are accelerating, a steady
state takes place at transport equilibrium. A set of simulations have been
performed in which the Shields number is varied from $\theta=0.2$, close to
the threshold of motion, until $\theta=1$ corresponding to intense bedload
transport. Three restitution coefficients $e=0.5$, $e=0.7$ and $e=0.9$ have
been explored and frictional particles are considered with $\mu^{p}=0.4$. The
parameters of all simulations are summarized in table 2.

(a)

(b)

(c)

Figure 1: Fluid-DEM simulation example for $\theta=0.6$, $e=0.7$,
$\mu^{p}=0.4$. (a) Numerical setup of the DEM simulations. $N=12$ layers of
particles are deposited by gravity above a rough bed made of fixed particles.
The fluid free surface is set to $H_{f}=Nd+h_{w}$. (b) Packing fraction
profile normalized for visibility by one tenth of $\phi_{max}=0.635$ and
dimensionless particle and fluid velocity profiles. (c) Granular temperature
balance in Pa.s-1.

Figure 1b shows typical profiles of packing fraction, particle and fluid
velocity for the case $\theta=0.6$, $e=0.7$. The sediment phase shows a very
complex phenomenology, with continuous transition from a pure fluid phase
($z>17d$) to a quasi-static creeping regime ($z<8d$). This transition zone is
called the bedload layer in which most of the transport is taking place. At
the same time, both the fluid and particle velocity decrease when going down
into the bed. The quasi-static regime, below $z<8d$, is characterized by a
roughly constant maximal packing fraction $\phi=\phi_{max}$ and very slow
fluid and granular motions with exponentially decreasing velocity profiles
with depth [29, 10]. The quasi-static regime has a negligible contribution to
the overall transport and will be neglected in this work.

Figure 1c shows the profiles of granular temperature budget (eq. (17)):
production, diffusion, dissipation through collisions and drag dissipation
term. At steady state, the left hand side of the granular temperature equation
(17) vanishes and the production of temperature should be balanced by
diffusion, contact dissipation and drag dissipation. As it can be observed in
figure 1c, the black dashed line representing the sum of these three terms is
not superimposed to the production term indicating that the balance is not
completely closed. Pähtz et al., [50] had the same issue in their DEM
simulations and attributed it to the dissipation term $\Gamma$ due to the lack
of scale separation in the averaging procedure, i.e. the vertical
discretization $d_{z}=d/30<<d$. For this reason in the following, the
dissipation term will be computed in order to close the granular temperature
balance equation, $\Gamma=-\tau_{xz}^{p}\dot{\gamma}-\partial q/\partial
z-J_{int}$.

$\begin{array}[]{|c|c|c|c|c|}\hline\cr\theta&h_{w}&H_{f}&N_{cell}&\Delta z\\\
\hline\cr 0.2&6d&18.5d&200&0.0925d\\\ 0.4&12d&24.5d&120&0.204d\\\
0.6&18d&30.5d&120&0.254d\\\ 0.8&24d&36.5d&200&0.183d\\\
1&30d&42.5d&200&0.213d\\\ \hline\cr\end{array}$

Table 2: Physical and numerical parameters for each configuration. Shields
number, fluid depth, free surface position, number of cells and vertical
discretization in sedFoam simulations.

###  III.2 Direct comparison with the two-fluid simulations

Figure 2: Comparison between DEM (symbols) and two phase flow results with
initial model (full line) for $e=0.5$ (first row), $e=0.7$ (second row) and
$e=0.9$ (third row). (a-d-g) Packing fraction profiles, (b-e-h) particle
velocity and transport (inset) profiles and (c-f-i) granular temperature.

The two-fluid model can now be compared with the results of the coupled fluid-
DEM model. The fluid equations being the same in both models, a particular
focus will be given to the granular behavior of the particle bed. The same
setup as presented in section III.1 is simulated using the two-fluid model,
with parameters reported in table 2. The mesh is composed of $120$ or $200$
cells depending on the configuration (see table 2) in order to maintain a grid
resolution between $0.1d$ and $0.25d$. The numerical implementation is fully
described in Chauchat et al., [12], consisting of a Total Variation Diminution
(TVD) central scheme, with Euler implicit first order scheme for time
derivatives, Gauss linear scheme for spatial gradient operators, Gauss
limitedLinear for divergence operators and Gauss linear corrected for
laplacian operators. The details of these schemes can be found in the OpenFOAM
documentation and they are all second order in space. The boundary conditions
are set to cyclic in the stream-wise direction and empty in the vertical
direction. At the top boundary, the pressure of both phases is fixed to $0$
and a zero gradient is imposed to all other fields. The bottom boundary
condition is set to fixedFluxPressure for the pressure and to noslip for the
velocities. The CFL is fixed to $0.1$ with a maximal time step $\Delta
t_{max}=10^{-3}$s.

Figure 2 compares, for all Shields numbers and for the three restitution
coefficients, the profiles of packing fraction, particle velocity and granular
temperature computed with the DEM (symbols) and with the two-fluid model (full
lines). The two-fluid model mispredicts the shape of the packing fraction
profile and overestimates strongly the particle velocity and the granular
temperature for all Shields numbers. The observed errors are more important
for larger restitution coefficients. The bed flows deeper in the two-fluid
simulations than in the DEM results and a plateau of velocity in the dilute
part of the granular flow is predicted which is not observed in the DEM
results. As a consequence the width of the bedload layer is larger in the
continuum model than in the DEM (see inset of figure 2b). Similarly, the
granular temperature is overestimated by the continuum model, in particular
inside the granular bed, suggesting that granular fluctuating energy
dissipation is underestimated in the continuum simulations. As discussed in
the introduction, this is because particles are frictional and dissipate more
energy than what is predicted in the frictionless KT model.

In order to improve the results of the two-fluid model, some corrections to
the Garzó and Dufty, [25] model are needed. In particular this model does not
account for interparticle friction while, as it will be shown, it plays an
important role in the dynamics of the granular behavior. In addition, the
departure between both models in the dilute part of the granular flow
indicates that saltation plays an important role in the dynamics of the flow.
Following Chialvo and Sundaresan, [16], the DEM results will be used as a
guideline to determine modified coefficients of the KT to account for
frictional interactions and saltation.

###  III.3 Corrections to the Garzó and Dufty, [25] model

Following equations (29), (30), (31) and (32), all terms of the granular
temperature equation can be computed in the DEM and the KT laws (eq (19),
(20), (21), (22) and (25)) can be compared with the DEM results. Additionally
to previous simulations, a DEM simulation with frictionless particles
($\mu^{p}=0$, $\theta=0.4$, $e=0.7$) have also been performed to investigate
the role of friction on the granular behavior.

(a)

(b)

Figure 3: Radial distribution function computed in the DEM simulations :
$\theta=0.4$ ($\circ$), $\theta=0.6$ ($+$), $\theta=0.8$ ($\times$) and
$\theta=1$ ($\bigtriangleup$) for three restitution coefficients $e=0.5$ (dark
blue), $e=0.7$ (medium blue) and $e=0.9$ (light blue). Green square
corresponds to frictionless particles $\theta=0.4$ and $e=0.7$. Lines
correspond to equation (23) for : $a=0.58$, $\phi_{max}=0.635$ ( ); $a=0.58$,
$\phi_{max}=0.59$ ( ); $a=2.71$, $\phi_{max}=0.635$ ( )

.

As mentioned previously, the radial distribution function is a key parameter
of the KT. From the DEM simulations, it can be estimated using the expression
of $F_{1}$ and by reversing the pressure law (eq. (19)) as,

$g_{0}(\phi)=\dfrac{1}{2\phi(1+e)}\left(\dfrac{p^{kin}}{\rho^{p}\phi
T}-1\right).$ (33)

In the DEM simulations, it is not possible to isolate the kinetic contribution
from the elastic one in the granular pressure. The radial distribution
function is therefore computed in equation (33) based on the total pressure.
The elastic pressure is however predominant only for packing fraction very
close to the random close packing $\phi_{max}\sim 0.635$ and it is assumed
that assimilating the kinetic pressure to the total pressure has a negligible
effect on the results. Figure 3b shows the computed radial distribution
function in the DEM simulations (symbols) compared with classical expressions
of the literature (black lines). The DEM data show two distinct trends for
frictional and frictionless particles. The radial distribution function for
frictionless particles is in perfect agreement with the form proposed by
Chialvo and Sundaresan, [16] (eq. (23), dashed dotted line). This is quite
remarkable because expression (23) were obtained in very different
configuration – dry granular flow in simple shear cell – showing that our DEM
simulations are in line with the literature on this point. Vescovi et al.,
[60] proposed another expression of the radial distribution function for
frictionless spheres that is almost superimposed with equation (23) (not shown
in figure 3b for clarity).

The DEM data show that inter-particle friction impacts the radial distribution
function. As observed by da Cruz et al., [17], Chialvo et al., [15] and
Chialvo and Sundaresan, [16] in shear cell numerical simulations, the radial
distribution function diverges at lower values of packing fraction. To account
for this effect, Chialvo and Sundaresan, [16] simply suggests to let the
critical packing fraction depend on interparticle friction in equation (23)
with $\phi_{max}(\mu^{p}=0.4)=0.59$. Figure 3b (dashed line) however shows
that it is not sufficient to capture the increase of $g_{0}$ for intermediate
values of $\phi$. In addition, in the bedload configuration, the packing
fraction will always exceed $0.59$, even for frictional particles (see figure
1b). This is because in the bed, there is a transition from kinetic
(collisional) stresses to elastic stresses. Fitting the DEM data, and based on
the form of Chialvo and Sundaresan, [16], it is proposed to modify the
empirical coefficient in equation (23) to $a=2.71$ and letting
$\phi_{max}=0.635$. This is plotted in figure 3b (full line) showing good
agreements with the DEM data.

With the obtained radial distribution functions for frictionless ($a=0.58$)
and frictional ($a=2.71$) particles, it is possible to compare the KT closure
laws listed in the table 1 with the DEM results. This is done in figure 4,
where coefficients $F_{2}$, $F_{3}$, $F_{4}$ and $J_{int}$ are plotted as a
function of $\phi$. For all coefficients, the data points are superimposed at
low packing fraction but separated into two branches at high packing fraction
whether particles are frictional or not. This separation is mainly due to the
difference of radial distribution function. It indicates that interparticle
friction has a stronger effect than the restitution coefficient on the KT
closure laws.

The previous discussion has shown that interparticle friction impacts strongly
modifies the radial distribution function. It is also expected to influence
granular temperature dissipation. Indeed, at contact between particles, energy
is dissipated through inelasticity and friction. The Garzó and Dufty, [25]
model was derived for frictionless spheres and therefore does not account for
this second source of dissipation. This is clearly evidenced in figure 4a
showing the dimensionless dissipation computed in the DEM simulations.
Dissipation is clearly higher for frictional (blue symbols) than for
frictionless particles (green squares). Dissipation is a bit higher for lower
restitution coefficient (compare dark blue to light blue symbols) but it is
more sensitive to interparticle friction. The dissipation for frictionless
particles is very well predicted by the Garzó and Dufty, [25] model (green
line) but not for frictional particles (dashed line underestimates light blue
symbols). To account for inter-particle friction, Jenkins and Zhang, [39]
proposed to treat it similarly to inelasticity by a modification of the
restitution coefficient $e_{eff}(\mu^{p})\leq e$ where $e_{eff}$ is the
effective restitution coefficient. Based on discrete simulations, Chialvo and
Sundaresan, [16] proposed the following empirical expression,

$e_{eff}=e-f(\mu^{p})=e-\dfrac{3}{2}\mu^{p}\exp{(-3\mu^{p})}.$ (34)

The modified coefficient $F_{4}^{\prime}=F_{4}(1-e_{eff}^{2})/(1-e^{2})$ is
plotted in full colored lines for the three restitution coefficient and is now
in very good agreements with the DEM results.

(a)

(b)

(c)

(d)

Figure 4: KT coefficients computed in the DEM (symbols, same as in figure 3b)
compared with the modified KT coefficients (colored full line) for
frictionless (green) and frictional case with $e=0.9$ (dark blue), $e=0.7$
(medium blue) and $e=0.5$ (light blue). Laws with initial model and $e=0.9$
are plotted for comparison ( ). (a) dimensionless collisional dissipation, (b)
dimensionless viscosity, (c) dimensionless heat transport coefficient. (d)
Dimensionless granular temperature dissipation by the drag force compared with
equation (25) ( ) and with corrected model (eq. 39) ( ).

The dimensionless viscosity and heat transport coefficients are plotted in
figures 4b and c. The Garzó and Dufty, [25] coefficients are plotted in black
dashed lines for frictional particles and $e=0.9$. They are quantitatively
good for large packing fraction but fails at reproducing the DEM data at low
packing fraction. In this limit, the kinetic contributions $\eta_{k}^{*}$ and
$\kappa_{k}^{*}$ of both coefficients are predominant and the observed
discrepancies therefore originate from these contributions. Indeed, both
$\eta_{k}^{*}$ and $\kappa_{k}^{*}$ contain a non dependent term with $\phi$
(see table 1) and the dimensionless viscosity and heat transfer coefficients
are non-vanishing in the very dilute limit contrary to what is observed in the
DEM simulations. This discrepancies between theory and the DEM coefficients is
at the origin of the difference of behavior of particle velocity in the dilute
regime (see figure 2b). In particular, the constant viscosity for $\phi$ close
to zero leads to the formation of the velocity plateau.

We attribute this departure from the Garzó and Dufty, [25] model to saltation
which is the main mode of transport in the dilute part of the flow. Particles
have ballistic trajectories that are controlled by gravity and drag force and
therefore can not be described by the KT. Jenkins and Valance, [38] proposed a
continuum theory for saltation in aeolian transport in which the granular
viscosity is linear with the packing fraction and therefore goes to zero in
the dilute regime, as observed in our DEM simulations. In the DEM data, it is
possible to isolate the kinetic contributions (first term in the right hand
side of equations (29) and (30)) and, following Jenkins and Valance, [38], the
constant parts of $\eta_{k}^{*}$ and $\kappa_{k}^{*}$ are made linear with
$\phi$ as

$\displaystyle\eta_{k}^{*}$
$\displaystyle=\dfrac{48/(5\sqrt{\pi})\phi-2/5(1+e)(1-3e)\phi
g_{0}(\phi)}{\left(1-1/4(1-e)^{2}-5/24(1-e^{2})\right)g_{0}(\phi)},$ (35)
$\displaystyle\kappa_{k}^{*}$
$\displaystyle=\dfrac{2\left[576/(225\sqrt{\pi})\phi+3/5(1+e)^{2}(2e-1)\phi
g_{0}(\phi)\right]}{(1-7/16(1-e))(1+e)g_{0}(\phi)}.$ (36)

Note that these modifications also impact the contact contributions
$\eta_{c}^{*}$ and $\kappa_{c}^{*}$ as they are proportional to the kinetic
contributions. The modified coefficients are plotted in figures 4b,c for
frictional (blue lines) and frictionless (green line) particles and they
perfectly reproduce the DEM data. The proposed empirical modifications should
be verified in other configurations and parameterized following physical
arguments. This would represent a complete study that lie outside of the
present work objectives. However, it suggests that the modelling of saltation,
a priori out of the scope of the KT, could be included in this framework.

The work done by the drag force due to granular fluctuating motion
$J_{int}=n\left<f_{Di}v_{i}^{p\prime}\right>^{s}$ is plotted in figure 4d. As
expected, no dependency on the restitution coefficient nor on interparticle
friction is observed. The black dashed line corresponds to the closure of Ding
and Gidaspow, [20] in equation (25). This closure underestimates the effective
dissipation of fluctuating energy by the drag force for the whole range of
packing fraction. This is because it was assumed a linear drag force to derive
closure (25), while it is not the case in this configuration (or in any
turbulent flow application). To understand how the quadratic nature of the
drag force influences the dissipation of fluctuating energy, let us replace by
the drag force expression (11) in $J_{int}$,

$J_{int}=n\left<f_{Di}v_{i}^{p\prime}\right>^{s}=\dfrac{3}{4}\dfrac{\phi\rho^{f}(1-\phi)^{-\zeta}}{d}\left<C_{D}||\bm{u^{f}}-\bm{v^{p}}||(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}\right>^{s}.$
(37)
$J_{int}\propto\left<C_{D}||\bm{u^{f}}-\bm{v^{p}}||(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}\right>^{s}.$
(38)

In order to provide a closure, it is necessary to discuss the term inside the
averaging operator $\left<\right>^{s}$. Indeed, depending on the value of the
particle Reynolds number $Re_{p}$,
$C_{D}||\bm{u^{f}}-\bm{v^{p}}||(u_{i}^{f}-v_{i}^{p})$ can be linear with the
relative velocity (laminar regime), quadratic (turbulent regime) or in
between. While it is usually assumed that the drag is linear [20, 23], this
assumption can not be done in the bedload transport configuration. Indeed, the
particle Reynolds number ranges between $Re_{p}\sim 100$ in the bed to
$Re_{p}\sim 2500$ at the bed surface. The drag force is therefore quadratic at
the surface and in a transitional regime in the bed. Splitting the particle
velocity into a mean and fluctuating component and performing a Taylor
expansion to first order with respect to the mean values, the following
expression of the granular temperature dissipation through the drag force is
obtained (see appendix):

$J_{int}\sim-\phi(1-\phi)K\left(3+2\dfrac{C_{D}^{\infty}}{C_{D}}\right)T,$
(39)

where $C_{D}^{\infty}$ is the value of the drag coefficient for an infinite
particle Reynolds number. Note that this term is negative and corresponds to a
loss of energy for the granular phase. The term $2C_{D}^{\infty}/C_{D}$
represents an additional source of dissipation due to the quadratic nature of
the drag force. In order to derive this expression it has been necessary to
assume that
$||\bm{u^{f}}-\left<\bm{v^{p}}\right>^{s}||\sim|u_{x}^{f}-\left<v_{x}^{p}\right>^{s}|$
and that $\left<v_{x}^{p\prime}v_{x}^{p\prime}\right>^{s}\sim 2T$. These
assumptions are specific to configurations with a preferential shearing
direction with anisotropic velocity fluctuations. However computations
presented in appendix can be easily adapted to other cases. The corrected
closure (39) is plotted in colored full lines and compared with DEM data in
figure 4. It improves quantitatively the prediction of dissipated energy. It
shows a very good agreement with the DEM in the dilute regime but
underestimates slightly at high packing fractions. In this limit, the
influence of the drag is however very small and this underestimation should
not affect the balance of granular temperature.

All corrections to the Garzó and Dufty, [25] model proposed in this section
are summarized in the right column of table 1. The corrected model can now be
used to compare with the DEM simulations.

###  III.4 Evaluation of the corrected KT model

Figure 5: Comparison between DEM (symbols) and two phase flow results with
corrected model (full lines) for $e=0.5$ (first row), $e=0.7$ (second row) and
$e=0.9$ (third row). (a-d-g) Packing fraction profiles, (b-e-h) particle
velocity and transport (inset) profiles and (c-f-i) granular temperature.

The two-fluid model with corrections presented in the previous section is
tested against the DEM simulations and the results are presented in figure 5.
Both the packing fraction and velocity profiles are quantitatively improved
compared with the initial model (see figure 2). The results are better for
lower restitution coefficients and very acceptable even for $e=0.9$
considering the complexity of the problem.

The shape of the packing fraction profiles is nicely reproduced, in particular
the apparition of the "shoulder" of packing fraction around $\phi\sim 0.3$ for
large Shields numbers. The velocity profiles are quantitatively predicted by
the continuum model for all Shields numbers and restitution coefficients. In
particular, the depth of the flowing layer is perfectly captured. This is
quite remarkable as the depth at which the flow stops is a response of the
stress and energy balance between both phases and not due to a bottom boundary
condition. The plateau of velocity observed previously in the dilute regime
(figure 2) is not present anymore thanks to the correction applied on the
viscosity and diffusion laws. As a consequence of the well predicted packing
fraction and velocity profiles, the transport profiles are very well
reproduced (see inset of figure 5b and e). The thickness of the transport
layer, and its increase with the Shields number, are particularly well
predicted. The modified continuum model improves quantitatively the predicted
granular temperature (figure 5c, f and i). It is perfectly predicted for
$e=0.5$ but underestimated for larger restitution coefficient. Nevertheless,
the results for $e=0.9$ are still acceptable.

##  IV Discussions

The corrected model implemented in SedFoam has been shown to reproduce very
accurately the DEM results. The proposed corrections are essentially empirical
but the improvements they provide can be analyzed in order to interpret the
physics of the granular flow.

First, inter-particular friction has been shown to play an important role in
the granular behavior. It modifies the radial distribution function and
increases dissipation at contact. In our DEM simulations with frictional
particles, almost no influence of the restitution coefficient is observed. The
packing fraction and velocity profiles are completely superimposed whatever
the value of $e$ ranging from $e=0.5$ to $e=0.9$ (compare figure 5a, d and g
or figure 5b, e and h). It indicates that the dissipation of granular
temperature is governed by friction rather than inelasticity of collisions.
This is what is observed in figure 4 where the KT coefficients computed with
the DEM are more affected by interparticle friction rather than the
restitution coefficient. It is interesting to note that the Garzó and Dufty,
[25] model is able to reproduce this property of the flow, with KT
coefficients that are more affected by the change of $g_{0}$ for frictional
particles than by the restitution coefficient. This is particularly true for
the viscosity and granular temperature diffusivity laws (superimposed full
blue lines of figure 4b and c), but less evident for the dissipation law.
Indeed, $\Gamma$ is proportional to $1-e^{2}$ and therefore highly dependent
on the restitution coefficient when it is large $e>0.7$. When applying the
correction proposed by Jenkins and Zhang, [39] and Chialvo and Sundaresan,
[16] which consists in the reduction of the restitution coefficient to account
for friction ($e_{eff}<0.7$), the dissipation law is made less dependent on
$e$ as illustrated in figure 4a (full lines very close for frictional
particles). The KT model is therefore applied in a range of parameters that
barely depends on the restitution coefficient and is more affected by $g_{0}$.
This highlights the role played by inter-particular friction in the sediment
bed dynamics and shows the capacity of the KT model to represent such granular
flow.

(a)

(b)

Figure 6: (a) $\mu(I)$ and (b) $\phi(I)$ relations computed in the DEM
(symbols, same as in figure 3b) and compared with: the rheological law
proposed by Maurin et al., [49] (equations (1) and (2)) for bedload transport
( ) ; two phase-flow simulation with the initial model before any correction
for the case $e=0.9$, $\mu^{p}=0.4$ ( ); two phase-flow simulations with the
corrected model for the three restitution coefficients and $\theta=0.6$ ( ).

Kinetic theory is often criticized for its inability to predict dense granular
flows. Our simulations however show very good results even in the dense
regime. To be more precise, the rheology of the granular flow is investigated.
Figure 6 shows the $\mu$ versus $I$ and $\phi$ versus $I$ relations, where it
is recalled that $\mu=\tau^{p}/p^{p}$ is the effective friction coefficient
and $I=d\dot{\gamma}/\sqrt{P^{p}/\rho^{p}}$ is the inertial number. The DEM
data are barely dependent on the restitution coefficients and the frictional
nature of particles influence strongly the rheological properties of the
granular flow. The DEM data for frictional particles are superimposed for the
whole range of inertial number. The $\mu(I)$ rheology for bedload transport
(equations (1) and (2)) derived by Maurin et al., [49] (dashed-dotted line) is
retrieved in the dense part of the granular flow ($10^{-2}<I<5$). The
departure between the data points and the $\mu(I)$ rheology for small $I$
(resp. large $I$) is characteristic of the quasi-static regime (resp. the
dilute regime) that cannot be captured by a local rheology such as $\mu(I)$
and $\phi(I)$ rheology. The same observations can be made on the scaling
between the packing fraction and the inertial number (figure 6b) with
departure from the $\phi(I)$ law at low and large values of $I$. Both these
scaling laws can also be computed in the two-fluid simulations, in which the
effective friction coefficient is computed with the total stresses, i.e. the
sum of the elastic and the kinetic contributions. They are plotted in colored
full lines in figure 6. Surprisingly, the $\mu(I)/\phi(I)$ rheology is
retrieved in the dense part of the granular flow. The elastic stresses make
the effective friction coefficient to asymptotically join the static friction
coefficient $\mu_{s}$ in the limit of vanishing inertial number. At large
inertial number, the two-fluid model over estimates the friction coefficient
but as in the DEM, a deviation to the master curve is observed with a similar
trend. The $\phi(I)$ rheology is also very well reproduced in the entire range
of inertial number. The rheology obtained with the initial model, before
application of the corrections proposed in section III.3, is plotted in dashed
line and highlights the role of these corrections in the rheological behavior.
The dashed line shows two branches. In the dense regime, corresponding to the
lower branch, the model fits with the DEM data for frictionless particles. The
corrections to account for interparticle friction, including modification of
$g_{0}$ and restitution coefficient in the dissipation law, make the granular
flow less inertial improving quantitatively the results. In the dilute regime,
corresponding to the higher branch, the behavior is completely different with
the initial model for both scaling laws, with turning-point of the curves
leading to an hysteresis effect. The improvements with the corrected model are
clearly due to the modification of the kinetic contribution of the viscosity
law, i.e. $\eta_{k}^{*}$ in $F_{2}(\phi)$. Indeed, in the dilute regime, the
elastic stresses are zero, the effective friction coefficient may be expressed
as $\mu=F_{2}/\sqrt{F_{1}}I$, and the modification of $F_{2}$ makes it
possible to recover a qualitative behavior of the scaling law. A better
understanding of the kinetic contribution of the viscosity law would certainly
improve the rheological results in the dilute regime.

In the literature, the recent attempts to model dense granular flows with the
kinetic theory have used the extended kinetic theory, consisting in the
introduction of a correlation length in the dissipation term $\Gamma$. Our
results suggest that the reproduction of the dense granular flow regime over
an erodible bed is dominated by the competition between elastic-frictional and
collisional-kinetic stresses rather than the development of velocity
correlation at high volume fraction [34, 35, 4, 6]. Here this is the
hybridization of the kinetic theory with a frictional model which makes it
possible to recover the correct behavior in the dense regime.

##  V Conclusions

In this paper, the modelling of sediment transport in the collisional regime
has been studied. In order to focus on the sediment phase modelling, simple
fluid-particle interactions have been considered without turbulent
interactions. The sediment phase is modelled with a frictional-collisional
approach [40]. Comparisons with coupled fluid-DEM simulations have highlighted
that the classical Kinetic Theory model of Garzó and Dufty, [25] is not able
to represent correctly the behavior of the granular flow. It is shown that
inter-particle friction influences the radial distribution function $g_{0}$
and increases energy dissipation. In the dilute regime, saltation modifies the
viscosity and diffusivity laws. Finally the quadratic nature has to be taken
into account because it increases the granular energy dissipation. Based on
these observations, a modified KT model has been proposed that reproduces
quantitatively the DEM results in the entire depth structure of the granular
flow.

This work has also shown that, in the bedload configuration, the frictional-
collisional approach is consistent with the $\mu(I)$ rheology and represents a
first success to reproduce the dense granular flow rheology with the Kinetic
Theory. The approach presented in this paper could be easily adapted to other
granular configurations, in particular to any granular flows over an erodible
bed, and should also be able to reproduce the rheological properties of the
flow.

The proposed modifications on the viscosity and granular temperature
diffusivity laws to account for saltation are empirical and should be tested
in other configurations. The success of the corrected model to reproduce
particle velocity, packing fraction and rheology in the dilute regime however
shows that the KT is a relevant framework to model saltation continuously.

Finally, the proposed corrected model represents a strong basis to further
study turbulence-particle interactions and their modelling in more complex
configurations such as real turbulent sediment transport with coherent
structures or in presence of hydraulic structures such as pipes or bridges
piles.

## Acknowledgments

Most of the computations presented in this paper were performed using the
GENCI infrastructure under allocations A0060107567 and A0080107567 and the
GRICAD infrastructure. LEGI is member of Labex TEC21 (Investissements d’Avenir
Grant Agreement ANR-11-LABX-0030) and Labex Osug@2020 (Investissements
d’Avenir Grant Agreement ANR-10-LABX-0056). The authors would like to
acknowledge the financial support from the Agence Nationale de la Recherche
(ANR) through the project SheetFlow (ANR-18-CE01-0003).

The authors would also like to acknowledge the fruitful discussions on
granular flows with P. Frey and R. Maurin.

##  Appendix A Definition of Favre-averaged operators and computation of the
granular temperature dissipation by the drag force

The average operator of the two-fluid equations is presented here. Equations
3, 4, 5, 6, 7 and 8 describe the evolution of Favre-averaged quantities. It
involves the ensemble averaged operator defined as

$\widetilde{\gamma}=\dfrac{1}{N}\lim_{N\rightarrow\infty}\sum_{k=1}^{N}\gamma_{k},$
(40)

where $\gamma_{k}$ is a given realization $k$ of a quantity $\gamma$ and $N$
the total number of realizations. The Favre-averaged is then defined as a
concentration-weighted averaged and is therefore different depending on the
phase considered. For a quantity $\gamma$, the Favre-averaged over the fluid
and solid phases are defined respectively as

$\displaystyle\left<\gamma\right>^{f}=\dfrac{\widetilde{(1-\phi)\gamma}}{1-\widetilde{\phi}},$
(41)
$\displaystyle\left<\gamma\right>^{s}=\dfrac{\widetilde{\phi\gamma}}{\widetilde{\phi}}.$
(42)

In the main text, the brackets and tildes were dropped for clarity.
Considering the fluid and particle velocities, they can therefore be
decomposed into a Favre-averaged and a fluctuating component. In the present
configuration, the averaged velocities are non-zero only in the streamwise
direction. In addition, as the models used in this paper do not consider any
fluid velocity fluctuations, they are considered to be zero. The fluid and
particles velocities can therefore be expressed as

$\displaystyle\bm{u^{f}}$
$\displaystyle=\left(\left<u_{x}^{f}\right>^{f},0,0\right),$ (43)
$\displaystyle\bm{v^{p}}$
$\displaystyle=\left(\left<v_{x}^{p}\right>^{s}+v_{x}^{p\prime},v_{y}^{p\prime},v_{z}^{p\prime}\right).$
(44)

The aim of this section is also to compute the phase-averaged dissipation of
granular temperature by the drag force. Following [50] and replacing by the
form of drag force expression (11) and drag coefficient (12), this term is
computed as

$\displaystyle J_{int}$
$\displaystyle=n\left<f_{Di}v_{i}^{p\prime}\right>^{s}=\dfrac{3}{4}\dfrac{\phi\rho^{f}(1-\phi)^{-\zeta}}{d}\left<C_{D}||\bm{u^{f}}-\bm{v^{p}}||(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}\right>^{s},$
(45)
$\displaystyle=\dfrac{3}{4}\dfrac{\phi\rho^{f}(1-\phi)^{-\zeta}}{d}\left[C_{D}^{\infty}\left<||\bm{u^{f}}-\bm{v^{p}}||(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}\right>^{s}+\dfrac{24.4\nu^{f}}{d}\left<(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}\right>^{s}\right].$
(46)

Let us consider the second term into the brackets of equation (46). Replacing
with the decomposition of velocities (43) and (44), it can be expressed as,

$(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}=\left(\left<u_{x}^{f}\right>^{f}-\left<v_{x}^{p}\right>^{s}\right)v_{x}^{p\prime}-v_{i}^{p\prime}v_{i}^{p\prime},$
(47)

and applying the Favre-averaged over the solid phase, the first term vanishes
and the second is three times the granular temperature,

$\left<(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}\right>^{s}=-3T.$ (48)

Let us consider now the first term into the brackets of equation (46). The
norm of the relative velocity can be expressed as

$||\bm{u^{f}}-\bm{v^{p}}||=||\left<\bm{u^{f}}\right>^{f}-\left<\bm{v^{p}}\right>^{s}||\sqrt{1-2\dfrac{v_{x}^{p\prime}}{\left<u_{x}^{f}\right>^{f}-\left<v_{x}^{p}\right>^{s}}+\dfrac{v_{i}^{p\prime}v_{i}^{p\prime}}{\left(\left<u_{x}^{f}\right>^{f}-\left<v_{x}^{p}\right>^{s}\right)^{2}}},$
(49)

and performing a Taylor expansion at first order in fluctuating velocity, it
becomes

$||\bm{u^{f}}-\bm{v^{p}}||=||\left<\bm{u^{f}}\right>^{f}-\left<\bm{v^{p}}\right>^{s}||\left(1-\dfrac{v_{x}^{p\prime}}{\left<u_{x}^{f}\right>^{f}-\left<v_{x}^{p}\right>^{s}}\right)+o(v^{p\prime
2}).$ (50)

Multiplying by equation (47) yields

$||\bm{u^{f}}-\bm{v^{p}}||(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}=||\left<\bm{u^{f}}\right>^{f}-\left<\bm{v^{p}}\right>^{s}||\left[\left(\left<u_{x}^{f}\right>^{f}-\left<v_{x}^{p}\right>^{s}\right)v_{x}^{p\prime}-v_{x}^{p\prime}v_{x}^{p\prime}-v_{i}^{p\prime}v_{i}^{p\prime}+\dfrac{v_{i}^{p\prime}v_{i}^{p\prime}v_{x}^{p\prime}}{\left<u_{x}^{f}\right>^{f}-\left<v_{x}^{p}\right>^{s}}\right]+o(v^{p\prime
3}).$ (51)

The last term is a third order term and can therefore be neglected. Applying
the Favre-averaged over the solid phase, the first term into the brackets
vanishes and it becomes

$\left<||\bm{u^{f}}-\bm{v^{p}}||(u_{i}^{f}-v_{i}^{p})v_{i}^{p\prime}\right>^{s}=||\left<\bm{u^{f}}\right>^{f}-\left<\bm{v^{p}}\right>^{s}||\left(-\left<v_{x}^{p\prime}v_{x}^{p\prime}\right>^{s}-3T\right).$
(52)

Replacing equations (48) and (52) in the drag term expression (46), it is
expressed as

$\displaystyle J_{int}$
$\displaystyle=\dfrac{3}{4}\dfrac{\phi\rho^{f}(1-\phi)^{-\zeta}}{d}||\left<\bm{u^{f}}\right>^{s}-\left<\bm{v^{p}}\right>^{s}||\left[-C_{D}^{\infty}\left<v_{x}^{p\prime}v_{x}^{p\prime}\right>^{s}-C_{D}^{\infty}3T-\dfrac{24.4\nu^{f}}{d||\left<\bm{u^{f}}\right>^{s}-\left<\bm{v^{p}}\right>^{s}||}3T\right],$
(53)
$\displaystyle=\dfrac{3}{4}\dfrac{\phi\rho^{f}(1-\phi)^{-\zeta}}{d}||\left<\bm{u^{f}}\right>^{s}-\left<\bm{v^{p}}\right>^{s}||C_{D}\left(-\dfrac{C_{D}^{\infty}}{C_{D}}\left<v_{x}^{p\prime}v_{x}^{p\prime}\right>^{s}-3T\right),$
(54)

where $C_{D}$ is the drag coefficient computed with phase averaged velocities.

In the bedload configuration, velocity fluctuations are not isotropic and, as
shown in figure 7, the streamwise fluctuating energy can be approximately
estimated as $\left<v_{x}^{p\prime}v_{x}^{p\prime}\right>^{s}\sim 2T$. Finally
the drag interaction term becomes

$J_{int}=-\phi(1-\phi)K(3+2\dfrac{C_{D}^{\infty}}{C_{D}})T.$ (55)

The drag coefficient $K$ is expressed using the Favre-averaged velocities and
the term $2\dfrac{C_{D}^{\infty}}{C_{D}}$ represents a correction to account
for the quadratic nature of the drag force.

Figure 7: Streamwise fluctuating kinetic energy as a function of the granular
temperature computed in the DEM simulations. Same symbols as in figure 3b.
Dashed line is $2T$.

## References

  * Ancey, [2020] Ancey, C. (2020).  Bedload transport: A walk between randomness and determinism. Part 1. The state of the art.  Journal of Hydraulic Research, 58(1):1–17.
  * Berzi, [2014] Berzi, D. (2014).  Extended kinetic theory applied to dense, granular, simple shear flows.  Acta Mechanica, 225(8):2191–2198.
  * Berzi and Fraccarollo, [2013] Berzi, D. and Fraccarollo, L. (2013).  Inclined, collisional sediment transport.
  * Berzi and Jenkins, [2011] Berzi, D. and Jenkins, J. T. (2011).  Surface flows of inelastic spheres.  Physics of Fluids, 23(1):013303.
  * Berzi and Jenkins, [2015] Berzi, D. and Jenkins, J. T. (2015).  Steady shearing flows of deformable, inelastic spheres.  Soft Matter, 11(24):4799–4808.
  * Berzi et al., [2020] Berzi, D., Jenkins, J. T., and Richard, P. (2020).  Extended kinetic theory for granular flow over and within an inclined erodible bed.  Journal of Fluid Mechanics, 885.
  * Bonamy et al., [2021] Bonamy, C., Chauchat, J., Hsu, T.-J., Puig-Montella, E., Chassagne, R., Mathieu, A., Cheng, Z., Nagel, T., Salimi-Tarazouj, A., Tsai, B., and Higuera, P. (2021).  SedFoam/sedfoam: Release 3.2.  Zenodo.
  * Carnahan and Starling, [1969] Carnahan, N. F. and Starling, K. E. (1969).  Equation of State for Nonattracting Rigid Spheres.  The Journal of Chemical Physics, 51(2):635–636.
  * Chapman and Cowling, [1970] Chapman, S. and Cowling, T. (1970).  The Mathematical Theory of Non-Uniform Gases.  Cambridge: Cambrige university press edition.
  * Chassagne et al., [2020] Chassagne, R., Maurin, R., Chauchat, J., Gray, J. M. N. T., and Frey, P. (2020).  Discrete and continuum modelling of grain size segregation during bedload transport.  Journal of Fluid Mechanics, 895.
  * Chauchat, [2018] Chauchat, J. (2018).  A comprehensive two-phase flow model for unidirectional sheet-flows.  Journal of Hydraulic Research, 56(1):15–28.
  * Chauchat et al., [2017] Chauchat, J., Cheng, Z., Nagel, T., Bonamy, C., and Hsu, T.-J. (2017).  SedFoam-2.0: A 3-D two-phase flow numerical model for sediment transport.  Geoscientific Model Development, 10(12):4367–4392.
  * Chen and Wood, [1985] Chen, C. P. and Wood, P. E. (1985).  A turbulence closure model for dilute gas-particle flows.  The Canadian Journal of Chemical Engineering, 63(3):349–360.
  * Cheng et al., [2018] Cheng, Z., Chauchat, J., Hsu, T.-J., and Calantoni, J. (2018).  Eddy interaction model for turbulent suspension in Reynolds-averaged Euler–Lagrange simulations of steady sheet flow.  Advances in Water Resources, 111:435–451.
  * Chialvo et al., [2012] Chialvo, S., Sun, J., and Sundaresan, S. (2012).  Bridging the rheology of granular flows in three regimes.  Physical Review E, 85(2):021305.
  * Chialvo and Sundaresan, [2013] Chialvo, S. and Sundaresan, S. (2013).  A modified kinetic theory for frictional granular flows in dense and dilute regimes.  Physics of Fluids, 25(7):070603.
  * da Cruz et al., [2005] da Cruz, F., Emam, S., Prochnow, M., Roux, J.-N., and Chevoir, F. (2005).  Rheophysics of dense granular materials: Discrete simulation of plane shear flows.  Physical Review E, 72(2):021309.
  * Dalla Valle, [1943] Dalla Valle, J. M. (1943).  Micromeritics: The Technology of Fine Particles.  Pitman.
  * Danon et al., [1977] Danon, H., Wolfshtein, M., and Hetsroni, G. (1977).  Numerical calculations of two-phase turbulent round jet.  International Journal of Multiphase Flow, 3(3):223–234.
  * Ding and Gidaspow, [1990] Ding, J. and Gidaspow, D. (1990).  A bubbling fluidization model using kinetic theory of granular flow.  AIChE Journal, 36(4):523–538.
  * Dong and Zhang, [1999] Dong, P. and Zhang, K. (1999).  Two-phase flow modelling of sediment motions in oscillatory sheet flow.  Coastal Engineering, 36(2):87–109.
  * Einstein, [1937] Einstein, H. A. (1937).  Der Geschiebetrieb als Wahrscheinlichkeitsproblem (Bedload transport as a probability problem).  (English translation by W.W. Sayre, in Sedimentation (Symposium to honour H.A. Einstein), edited by H.W. Shen, Fort collins, Colorado, 1972, C1-C105). Zürich: ETHZ.
  * Fox, [2014] Fox, R. O. (2014).  On multiphase turbulence models for collisional fluid–particle flows.  Journal of Fluid Mechanics, 742:368–424.
  * Frey, [2014] Frey, P. (2014).  Particle velocity and concentration profiles in bedload experiments on a steep slope.  Earth Surface Processes and Landforms, 39(5):646–655.
  * Garzó and Dufty, [1999] Garzó, V. and Dufty, J. W. (1999).  Dense fluid transport for inelastic hard spheres.  Physical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, 59(5 Pt B):5895–5911.
  * GDR MiDi, [2004] GDR MiDi (2004).  On dense granular flows.  The European Physical Journal E, 14(4):341–365.
  * Goldhirsch, [2010] Goldhirsch, I. (2010).  Stress, stress asymmetry and couple stress: From discrete particles to continuous fields.  Granular Matter, 12(3):239–252.
  * Gonzalez-Ondina et al., [2018] Gonzalez-Ondina, J. M., Fraccarollo, L., and Liu, P. L.-F. (2018).  Two-level, two-phase model for intense, turbulent sediment transport.  Journal of Fluid Mechanics, 839:198–238.
  * Houssais et al., [2015] Houssais, M., Ortiz, C. P., Durian, D. J., and Jerolmack, D. J. (2015).  Onset of sediment transport is a continuous transition driven by fluid shear and granular creep.  Nature Communications, 6:6527.
  * Hsu et al., [2004] Hsu, T.-J., Jenkins, J. T., and Liu, P. L.-F. (2004).  On two-phase sediment transport: Sheet flow of massive particles.  Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 460(2048):2223–2250.
  * Hsu and Liu, [2004] Hsu, T.-J. and Liu, P. L.-F. (2004).  Toward modeling turbulent suspension of sand in the nearshore.  Journal of Geophysical Research: Oceans, 109:C06018.
  * Jackson, [1997] Jackson, R. (1997).  Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid.  Chemical Engineering Science, 52(15):2457–2469.
  * Jackson, [2000] Jackson, R. (2000).  The Dynamics of Fluidized Particles.  Cambridge University Press.
  * Jenkins, [2006] Jenkins, J. T. (2006).  Dense shearing flows of inelastic disks.  Physics of Fluids, 18(10):103307.
  * Jenkins, [2007] Jenkins, J. T. (2007).  Dense inclined flows of inelastic spheres.  Granular Matter, 10(1):47–52.
  * Jenkins and Hanes, [1998] Jenkins, J. T. and Hanes, D. M. (1998).  Collisional sheet flows of sediment driven by a turbulent fluid.  Journal of Fluid Mechanics, 370:29–52.
  * Jenkins and Savage, [1983] Jenkins, J. T. and Savage, S. B. (1983).  A theory for the rapid flow of identical, smooth, nearly elastic, spherical particles.  Journal of Fluid Mechanics, 130:187–202.
  * Jenkins and Valance, [2018] Jenkins, J. T. and Valance, A. (2018).  Two-phase continuum theory for windblown sand.  Physical Review Fluids, 3(3):034305.
  * Jenkins and Zhang, [2002] Jenkins, J. T. and Zhang, C. (2002).  Kinetic theory for identical, frictional, nearly elastic spheres.  Physics of Fluids, 14(3):1228–1235.
  * Johnson and Jackson, [1987] Johnson, P. C. and Jackson, R. (1987).  Frictional–collisional constitutive relations for granular materials, with application to plane shearing.  Journal of Fluid Mechanics, 176:67–93.
  * Lee et al., [2016] Lee, C.-H., Low, Y. M., and Chiew, Y.-M. (2016).  Multi-dimensional rheology-based two-phase model for sediment transport and applications to sheet flow and pipeline scour.  Physics of Fluids, 28(5):053305.
  * Li and Sawamoto, [1995] Li, L. and Sawamoto, M. (1995).  Multi-Phase Model on Sediment Transport in Sheet-Flow Regime Under Oscillatory Flow.  Coastal Engineering in Japan, 38(2):157–178.
  * Lun and Savage, [1986] Lun, C. K. K. and Savage, S. B. (1986).  The effects of an impact velocity dependent coefficient of restitution on stresses developed by sheared granular materials.  Acta Mechanica, 63(1):15–44.
  * Lun et al., [1984] Lun, C. K. K., Savage, S. B., Jeffrey, D. J., and Chepurniy, N. (1984).  Kinetic theories for granular flow: Inelastic particles in Couette flow and slightly inelastic particles in a general flowfield.  Journal of Fluid Mechanics, 140:223–256.
  * Ma and Ahmadi, [1986] Ma, D. and Ahmadi, G. (1986).  An equation of state for dense rigid sphere gases.  The Journal of Chemical Physics, 84(6):3449–3450.
  * Mathieu et al., [2021] Mathieu, A., Chauchat, J., Bonamy, C., Balarac, G., and Hsu, T.-J. (2021).  A finite-size correction model for two-fluid large-eddy simulation of particle-laden boundary layer flow.  Journal of Fluid Mechanics, 913.
  * Maurin, [2015] Maurin, R. (2015).  Investigation of Granular Behavior in Bedload Transport Using a Eulerian-Lagragian Model. PhD thesis, Univ. Grenoble Alpes.
  * Maurin et al., [2015] Maurin, R., Chauchat, J., Chareyre, B., and Frey, P. (2015).  A minimal coupled fluid-discrete element model for bedload transport.  Physics of Fluids, 27(11):113302.
  * Maurin et al., [2016] Maurin, R., Chauchat, J., and Frey, P. (2016).  Dense granular flow rheology in turbulent bedload transport.  Journal of Fluid Mechanics, 804:490–512.
  * Pähtz et al., [2015] Pähtz, T., Durán, O., Ho, T.-D., Valance, A., and Kok, J. F. (2015).  The fluctuation energy balance in non-suspended fluid-mediated particle transport.  Physics of Fluids, 27(1):013303.
  * Pouliquen et al., [2006] Pouliquen, O., Cassar, C., Jop, P., Forterre, Y., and Nicolas, M. (2006).  Flow of dense granular material: Towards simple constitutive laws.  Journal of Statistical Mechanics: Theory and Experiment, 2006(07):P07020.
  * Prandtl, [1926] Prandtl, L. (1926).  Bericht über neuere Turbulenzforschung.  Hydraulische Probleme. Vorträge Hydrauliktagung Göttingen, 5:1–13.
  * Rao and Nott, [2008] Rao, K. K. and Nott, P. R. (2008).  An Introduction to Granular Flow.  Cambridge Series in Chemical Engineering. Cambridge University Press, Cambridge.
  * Revil-Baudard and Chauchat, [2013] Revil-Baudard, T. and Chauchat, J. (2013).  A two-phase model for sheet flow regime based on dense granular flow rheology.  Journal of Geophysical Research: Oceans, 118(2):619–634.
  * Roux and Combe, [2002] Roux, J. and Combe, G. (2002).  Quasistatic rheology and the origins of strain.  Comptes Rendus Physique, 3(2):131–140.
  * Savage and Jeffrey, [1981] Savage, S. B. and Jeffrey, D. J. (1981).  The stress tensor in a granular flow at high shear rates.  Journal of Fluid Mechanics, 110:255–272.
  * Schwager and Poschel, [2007] Schwager, T. and Poschel, T. (2007).  Coefficient of restitution and Linear–Dashpot model revisited.  Granular Matter, 9(6):465–469.
  * Smilauer et al., [2015] Smilauer et al., V. (2015).  Yade Documentation 2nd Ed. The Yade Project. Zenodo.
  * Torquato, [1995] Torquato, S. (1995).  Nearest-neighbor statistics for packings of hard spheres and disks.  Physical Review E, 51(4):3170–3182.
  * Vescovi et al., [2014] Vescovi, D., Berzi, D., Richard, P., and Brodu, N. (2014).  Plane shear flows of frictionless spheres: Kinetic theory and 3D soft-sphere discrete element method simulations.  Physics of Fluids, 26(5):053305.