# Exner diagnosis method for two-fluid morphodynamics simulations

Alban Gilletta, Cyrille Bonamy, Julien Chauchat

► To cite this version:

Alban Gilletta, Cyrille Bonamy, Julien Chauchat. Exner diagnosis method for two-fluid morphodynamics simulations. International Journal of Multiphase Flow, 2026, 194, pp.105457. 10.1016/j.ijmultiphaseflow.2025.105457. hal-05299771

HAL Id: hal-05299771

https://hal.science/hal-05299771v1

Submitted on 6 Oct 2025

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
<figure><img src="image_2.png" /></figure>

Research Paper

# Exner diagnosis method for two-fluid morphodynamics simulations

Alban Gilletta<sup>a,b ID,*</sup>, Cyrille Bonamy<sup>b ID</sup>, Julien Chauchat<sup>c ID</sup>

## ARTICLE INFO

Dataset link: https://github.com/SedFoam/exnerDiagnosis.git

**Keywords:**

Bedload

Erosion processes

Sedimentation

Suspended sediments

Water-sediment interface interactions

## ABSTRACT

Mass balance in sediment transport is often associated with single-phase flow model through the resolution of the Exner equation. In this contribution, the mass balance is derived from the two-phase flow, two-fluid, model equations and compared with the conventional single-phase flow approach. Five terms are identified from the two-phase flow formulation: the bed evolution, the storage evolution, the lateral flux, the bed flux and the top flux. Three different flow configurations with increasing complexity are tested, namely a one dimensional pure sedimentation, a two dimensional scour downstream an apron and a three dimensional scour around a vertical square cylinder. The mass balance formulation identifies the driving process of an erosion or a deposition process in each case. This approach can be applied on other configurations related with sediment flux and morphology bed evolution such as scour process.

## 1. Introduction

Morphological processes occur across various aquatic environments, including rivers, submarine dunes, and localized areas around structures affected by scour. The sediment bed constantly seeks an equilibrium, influenced by factors such as flow rate, sediment load, and slope. This dynamic balance is famously illustrated by the equilibrium concept introduced by Lane (1955) and further detailed by Borland (1960). Changes in any of these input factors lead to shifts in morphological evolution.

The prediction of the bed elevation evolution $z_b(x, y, t)$ has long been grounded in the Exner equation, introduced a century ago (Exner, 1925). As numerous studies have shown, this equation represents the sediment mass balance in the bed and it is essential for understanding the fundamental processes driving morphological evolution of the sediment bed (e.g., Paola and Voller, 2005; Coleman and Nikora, 2009; Garcia, 2008). Large-scale numerical simulations apply the Exner equation to forecast bed evolution; however, these predictions are often unreliable due to uncertainties in assessing one key input: the solid discharge. In these large-scale models, it is calculated using empirical formulae based on hydrodynamic flow parameters (e.g., Meyer-Peter and Müller, 1948; Wong and Parker, 2006; Engelund and Fredsoe, 1976). Another important quantity is the vertical, or bed-normal, sediment flux, the so-called pick-up flux for which bed-shear stress based empirical formulae are used (e.g. Van Rijn, 1984a). The main weakness

of this approach is that it relies on empirical laws, which require parameter adjustment on case-by-case basis to match observations.

In recent decades a new approach has emerged, the two-phase flow approach, which owes its appeal to its ability to resolve the physical processes essential to sediment transport, often driven with turbulent interactions in natural environment. The Eulerian-Lagrangian method, which treats the fluid phase as a continuum phase (Eulerian) while tracking the trajectory of each particle within the fluid (Lagrangian), plays a key role in defining sediment-turbulence interaction regimes (Finn and Li, 2016; Balachandar, 2009). This method combined with DNS or LES has also proven effective in resolving erosion and deposition phenomena (e.g. Zhang et al., 2015; Mazzuoli et al., 2024b; Kidanemariam and Uhlmann, 2014). However, such an approach is computationally demanding, and despite advances in computing power, simulations remain impractical for modeling larger-scale morphological processes. The Eulerian-Eulerian method or two-fluid approach consists in considering the fluid and the sediment phase as two interpenetrating continua for which momentum and continuity equations are solved for both the fluid and solid phases (e.g. Jenkins and Hanes, 1998; Hsu et al., 2003; Chauchat and Guillou, 2008; Lee et al., 2016; Chauchat et al., 2017), allows to circumvent some of the limitations associated with the Eulerian-Lagrangian approach. In this two-fluid approach, the sediment phase concentration $\phi$ is solved from the free surface down to the immobile bed. Therefore, the position of the sediment bed is not explicitly defined in the two-fluid approach.

<footnote><sup>a</sup> France Énergies Marines (FEM), 525 Avenue Alexis de Rochon, Plouzané, 29280, Finistère, France</footnote>

<footnote><sup>b</sup> Laboratoire des Écoulements Géophysiques et Industriels (LEGI), University of Grenoble Alpes, G-INP, CNRS, Grenoble, 38000, Isère, France</footnote>

<footnote><sup>c</sup> University of Grenoble Alpes, INRAE, CNRS, G-INP, IRD, IGE, Grenoble, 38000, Isère, France</footnote>

<footnote>* Corresponding author at: Laboratoire des Écoulements Géophysiques et Industriels (LEGI), University of Grenoble Alpes, G-INP, CNRS, Grenoble, 38000, Isère, France.</footnote>

<footnote>E-mail address: alban.gilletta-de-saint-joseph@univ-grenoble-alpes.fr (A. Gilletta).</footnote>
In Nagel et al. (2020) the authors suggested using the following proxy: $z_b = \{z|\phi(x, y, z, t) = \phi_b\}$ where $\phi_b = 0.57$ is taken as the critical volume fraction below which particles can be considered immobile. The main advantage of the two-fluid approach is that it explicitly resolves the sediment fluxes from first principles. Chauchat (2018) has shown that for steady uniform flow conditions the sediment fluxes predicted by two-phase flow models are within the experimental uncertainties, strongly supporting the idea that these models can be used to solve three-dimensional problems involving a bed surface such as scour (Nagel et al., 2020; Gilletta et al., 2025), wave-driven sediment transport (Mathieu et al., 2022), orbital wave ripples (Salimi-Tarazouj et al., 2021b,a, 2024) or underwater avalanches (Montellà et al., 2021, 2023).

An interesting perspective for the two-phase flow approach would be to use the numerical simulation results to perform an upscaling of the sediment fluxes in order to better understand the underlying physics and the key control parameters beyond the fluid bed shear stress. To achieve this goal, a clear methodology for establishing the links between the single-phase and two-phase flow models is needed. The aim of this paper is to fill this gap by proposing a mathematically rigorous derivation of the mass balance diagnostic, or Exner equation, derived from the sediment phase continuity equation of the two-phase flow model. This new methodology will be tested on a number of different configurations of increasing complexity. It should also be mentioned that the proposed methodology would be of benefit to the wider multiphase flow community working on similar physical problems where particles are deposited on surfaces or eroded from a substrate (e.g. Chou et al., 2014; Ozel et al., 2013; Chen et al., 2011). This is a ubiquitous situation in process and chemical engineering. It must be stated clearly that the goal of the present work is not to implement a single-phase flow model and as such no details about the closures for this class of model will be provided herein.

The formulations of the Exner equation for both approaches and the methodology chosen for the mass balance diagnostic are presented is Section 2. This methodology is tested on three numerical configurations with increasing complexity in Section 3. Each contribution of the mass balance equation are discussed to better understand their role in the simulation. Finally, conclusions are drawn in Section 4.

## 2. Mass balance diagnostic formulation

### 2.1. Single-phase flow formulation

In the single-phase flow model, hydrodynamics and morphodynamic evolution are usually solved separately. Flow is solved using three-dimensional (3D) Navier-Stokes equations, which are based on mass and momentum conservation. Sediment motion is described by the combined effects of advection by the mean flow, gravitational settling in the absence of flow and a turbulent diffusive flux, which accounts for turbulence-particle interactions. Sediment transport is coupled with fluid flow via an advection-diffusion equation. Therefore, this model is not strictly a single-phase flow approach; rather, it is a mixture model based on the passive scalar approximation. For the sake of simplicity, however, this class of model will be denoted as a single-phase flow formulation, in contrast to the two-phase flow formulation presented in the next subsection. Morphodynamic evolution is obtained from the Exner equation, which represents sediment mass conservation in the bed layer. This equation describes how the bed evolves over time in response to flow solicitation through bed shear stress and bed slope. Bed evolution is coupled with a moving mesh that influences the hydrodynamics. In most models, however, bed evolution is not directly coupled with hydrodynamics or sediment transport, which prevents the applicability of the approach in cases of fast bed evolution. The historical form of the Exner equation (Exner, 1925) is written as follows:

$$ \frac{\partial z_b}{\partial t}(x, y, t) + A \nabla \cdot \bar{u}(x, y, t) = 0, \quad (1) $$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Sketch of the control volume V subjected to mass balance equation. In red: $\Sigma$ = closed surface of V; in green: $S_l$ = lateral surface, $l$ = contour line of $S_l$ and $n_l$ = normal unit vector of $S_l$; in orange: $S_b$ = bed surface, $z_b$ = bed elevation and $n_b$ = normal unit vector of $S_b$; in purple: $z_{ref}$ = suspension reference elevation; in blue: $S_l$ = top surface.</figcaption></figure>

where $z_b$ is the bed elevation, $\bar{u}$ is the depth average flow velocity vector, $A$ is a coefficient depending on the river features in length unit (m) and $\nabla \cdot \bar{u} \equiv \frac{\partial n_l}{\partial x_l}$ stands for the divergence operator.

Nowadays, the Exner equation is often written with the sediment flux vector $q_{tl}$ that is related to the average flow velocity and could be rewritten as follows:

$$ (1 - p) \frac{\partial z_b}{\partial t} + \nabla \cdot q_{tl} = 0, \quad (2) $$

where $p$ is the porosity of the bed that can be written as $p = 1 - \phi_b$ with $\phi_b$ the solid volume fraction in the bed. Fig. 1 is a sketch depicting the different variables in the control volume over which the mass balance budget is performed.

For more recent single-phase flow models, the morphology of the bed includes the time evolution of the sediment concentration defined in Eq. (5) and separates the total-load $q_{tl}$ into two components, the bed-load $q_{bl}$ generally computed explicitly from the hydrodynamic quantities, and suspended-load $q_{sl}$ defined in Eq. (3).

$$ q_{sl}(x, y, t) = \int_{z_{ref}(x, y, t)}^{h} \phi(x, y, t) u^f(x, y, z, t) dz \quad (3) $$

where $u^f$ is the local flow velocity and $z_{ref}$ is a reference elevation associated with the interface between the bed-load and suspended-load. All variables introduced for the single-phase flow model are illustrated in the left panel of Fig. 3.

With these definitions, a first formulation of the Exner equation can be expressed as follows:

$$ \iint_{S_b} \phi_b \frac{\partial z_b}{\partial t} dS_b + \oint_l (q_{bl} + q_{sl}) \cdot n_l dl + \frac{d\Psi}{dt} = 0, \quad (4) $$

where $\Psi$ is the solid fraction integrated over the volume V and defined in Eq. (5), $l$ is the contour of the lateral surface of the control volume V and $n_l$ is the unit normal vector of this contour. Note that the Green theorem has been used to transform the divergence operator. $S_b$ is the bed surface corresponding to the bottom mesh boundary condition which distort with time. Eq. (4) comes down to a mass balance equation over the water column between $z_b$ and $h$ (see left panel of Fig. 3).

$$ \Psi(t) = \iint_{S_b} \int_{z_b(t)}^{h} \phi(x, y, z, t) dz dS_b = \iiint_{V(t)} \phi(x, y, z, t) dV \quad (5) $$

The following reasoning aims at writing the Exner equation in another formalism commonly used in the single-phase flow approach. For this purpose, an advection-diffusion equation is added to resolve the suspended-load:

$$ \frac{\partial \phi}{\partial t} + \left( u^f + w_s \frac{g}{\|g\|} \right) \cdot \nabla \phi - \nabla \cdot (\epsilon_s \nabla \phi) = 0, \quad (6) $$
where $w_s$ in m s⁻¹ is the settling velocity parameterized with an empirical law, $g$ in m s⁻² is the gravity vector and $\epsilon_s$ in m² s⁻¹ is the sediment concentration turbulent diffusivity. From this equation projected in the wall-normal direction, the bottom boundary condition gives:

$$ - \epsilon_s \frac{\partial \phi}{\partial z} (z_{ref}) - (w^f - w_s) \phi (z_{ref}) = E (z_{ref}) - D (z_{ref}), \quad (7) $$

where $w^f$ is the vertical component of the flow velocity, $\phi(z_{ref})$ is the reference solid volume fraction calculated at the reference elevation $z_{ref}$, $E$ is the erosion flux also called the pick-up flux and $D$ the deposition flux defined at $z_{ref}$. The deposition flux at the interface $z_{ref}$ is simply the product of the vertical velocity and the volume fraction: $D = (w_s - w^f) \phi(z_{ref})$. The reference elevation can be assumed to be proportional to the mean diameter such that $z_{ref} - z_b = \alpha d_{50}$ (e.g., Cekil and Rodi, 1988; Renaud et al., 2025). Either the distribution of $\phi$ is known at the elevation $z = z_{ref}$ with Eq. (6) if the grid resolution is fine enough or the value at $z = z_{ref}$ is parameterized with the first grid cell numerically resolved, or a Rouse profile is assumed to extrapolate the value.

The erosion flux has to be modeled, the method proposed by Wu et al. (2000) consists in assuming an equilibrium concentration $\phi_{eq}$ at $z_{ref}$ empirically determined (Van Rijn, 1984b) such that $E(z_{ref}) = w_s\phi_{eq}$. It corresponds to a concentration at the top of the bed-load layer that would compensate the deposition flux such that there is no influence of the suspended-load.

The net sediment flux $F$ corresponds to $E - D$ is sometimes introduced, the sign of which indicates if there is erosion or deposition. In other words, assuming that $w^f$ is negligible, $E - D = w_s(\phi_{eq} - \phi(z_{ref}))$. If $\phi^{eq} > \phi(z_{ref})$, the actual concentration $\phi(z_{ref})$ is lower than the equilibrium concentration $\phi_{eq}$, meaning that an erosion occurs. Conversely if $\phi_{eq} < \phi(z_{ref})$, the actual concentration $\phi(z_{ref})$ is greater than the equilibrium concentration $\phi_{eq}$, meaning that a deposition occurs. This term $E - D$ can be associated with $\oint_l q_{st} \cdot n_l dl + \frac{d}{dt} \int_{z_{ref}}^h \phi dz$ so that Eq. (4) can be rewritten as follows:

$$ \iint_{S_b} \phi_b \frac{\partial z_b}{\partial t} dS_b + \frac{d}{dt} \iint_{S_b} \int_{z_b(t)}^{z_{ref}(t)} \phi dz dS_b + \oint_l q_{bl} \cdot n_l dl + \iint_{S_b} (E - D) dS_b = 0, \quad (8) $$

where the second term on the left hand side corresponds to the storage evolution in the bed-load layer.

This equation is equivalent to a mass balance equation applied in a control volume defined by the bed-load layer between $z_b$ and $z_{ref}$ (see left panel Fig. 3).

In this section, we have provided a brief overview of the Exner equation within the context of single-phase flow. The next subsection will present a complete derivation of the equivalent Exner equation from the two-phase model governing equations.

## 2.2. Two-phase flow formulation

In the two-phase flow model, the mass and momentum conservation equations are resolved for the fluid and solid phases.

$$ \frac{\partial (1 - \phi)}{\partial t} + \nabla \cdot ((1 - \phi) \mathbf{u}^f) = 0 \quad (9) $$

$$ \frac{\partial \phi}{\partial t} + \nabla \cdot (\mathbf{u}^s \phi) = 0 \quad (10) $$

$$ \frac{\partial \rho_s \phi u^s}{\partial t} + \nabla \cdot (\rho_s \phi u^s \otimes u^s) = -\phi \nabla p^f + \phi f - \nabla \tilde{p}^s + \nabla \cdot \tilde{\tau}^s + \phi \rho_s g + f_D, \quad (11) $$

$$ \frac{\partial \rho_f (1 - \phi) \mathbf{u}^f}{\partial t} + \nabla \cdot (\rho_f (1 - \phi) \mathbf{u}^f \otimes \mathbf{u}^f) = -(1 - \phi) \nabla p^f + (1 - \phi) f + \nabla \cdot \tilde{\tau}^f + (1 - \phi) \rho_f g - f_D, \quad (12) $$

where $\phi$ is the solid volume fraction, $u^s$ and $u^f$ are the solid particle and fluid velocity vectors respectively. $\rho_s$ and $\rho_f$ are the solid particle and fluid density respectively, $g$ is the gravitational acceleration vector, $p^f$ is the fluid pressure, $f$ is an external force that drives the flow. The fluid stress $\tilde{\tau}^f$ is compounded by the fluid viscous stress and the fluid Reynolds stresses and $\tilde{p}^s$, $\tilde{\tau}^s$ are the solid normal stress and shear stress. Fluid–particle interactions are accounted for through drag $f_D$ and buoyancy $\phi \nabla p^f$ forces. With appropriate closures, the two-phase flow equations allow both dilute and dense regions of the flow to be handled. This framework enables the continuous evolution of the solid volume fraction from the water surface to the bottom of the porous medium. The interested reader is referred to the extensive literature on this subject (e.g. Jenkins and Hanes, 1998; Hsu et al., 2004; Revil-Baudard and Chauchat, 2013; Lee et al., 2016; Chauchat and Lacaze, 2024; Chauchat and Mazzuoli, 2024; Mazzuoli et al., 2024a). The Exner equation is not solved explicitely, instead the mass conservation equation for the solid phase (10) is resolved.

Since the Exner equation derives from a mass balance for the sediments, it should be verified in the two-phase flow model as well.

By integrating Eq. (10) over a control volume $V$ delimited by the bed interface corresponding to a sufficiently densely packed bed, usually arbitrary fixed by an iso-concentration $\phi_b$, and $h$ the water height considered fixed over time for simplicity, the following equation is obtained:

$$ \iiint_{V(t)} \frac{\partial \phi}{\partial t} dV + \iiint_{V(t)} \nabla \cdot (\mathbf{u}^s \phi) dV = 0 \quad (13) $$

The Leibniz's rule in the three dimensional case reads:

$$ \iiint_{V(t)} \frac{\partial \phi}{\partial t} dV = \frac{d}{dt} \iint_{V(t)} \phi dV - \iint_{\Sigma(t)} \phi \mathbf{u}^b \cdot \mathbf{n} d\Sigma \quad (14) $$

where $\mathbf{u}^b$ is the interface velocity vector and $\mathbf{n}$ is the normal unit vector facing outwards of the closed surface $\Sigma(t)$. As the mesh is fixed in this framework, the interface of the bed can be seen as an interface moving only locally through the vertical axis at each location. While it is quite obvious in a discrete environment, it can be proved in a continuous environment that: $\mathbf{u}^b \cdot \mathbf{n}_b dS = -\frac{\partial z_b}{\partial t} dS_b$ (see proof in Appendix A), where $dS$ is an infinitesimal element of bed surface defined by the function $f(x, y) = z_b$ and $dS_b$ is the projected surface of $dS$ on the horizontal plane. The positive sign of the first term is due to the normal unit vector of the bottom surface pointing downward. Combining Eqs. (13) and (14), the following form of the solid phase mass conservation equation is obtained:

$$ \iint_{S_b} \phi \frac{\partial z_b}{\partial t} dS_b + \frac{d}{dt} \iiint_{V(t)} \phi dV + \iiint_{V(t)} \nabla \cdot (\mathbf{u}^s \phi) dV = 0 \quad (15) $$

The Green–Ostrogradsky theorem enables to rewrite the last term as follows:

$$ \iiint_{V(t)} \nabla \cdot (\mathbf{u}^s \phi) dV = \iint_{\Sigma(t)} \phi \mathbf{u}^s \cdot \mathbf{n} d\Sigma \quad (16) $$

Leading to:

$$ \iint_{S_b} \phi \frac{\partial z_b}{\partial t} dS_b + \frac{d}{dt} \iiint_{V(t)} \phi dV + \iint_{\Sigma(t)} \phi \mathbf{u}^s \cdot \mathbf{n} d\Sigma = 0 \quad (17) $$

Because the closed surface $\Sigma(t)$ is defined by $S_l(t) \cup S_b(t) \cup S_t$ where $S_l$ is the lateral surface and $S_t$ is the top surface at the fixed water depth $h$ simplifying assumption, it gives:

$$ \underbrace{\iint_{S_b} \phi \frac{\partial z_b}{\partial t} dS_b}_{\textcircled{1} \text{ Bed evolution}} + \underbrace{\frac{d}{dt} \iiint_{V(t)} \phi dV}_{\textcircled{2} \text{ Storage evolution}} + \underbrace{\iint_{S_l} \phi u^s \cdot n_l dS_l}_{\textcircled{3} \text{ Lateral flux}} + \underbrace{\iint_{S_b} \phi u^s \cdot n_b dS_b}_{\textcircled{4} \text{ Bed flux}} + \underbrace{\iint_{S_l} \phi u^s \cdot z dS_l}_{\textcircled{5} \text{ Top flux}} = 0 \quad (18) $$

Introducing the following definition for the sediment transport flux, where $z_s$ corresponds with the interface between the suspension and
<figure><img src="image_4.png" /><figcaption>Fig. 2. Representation of the different term of Eq. (18). In orange: ① bed evolution; in blue: ② storage evolution; in green: ③ lateral flux; in magenta: ④ bed flux; in cyan: ⑤ top flux.</figcaption></figure>

bed-load layer and the storage term ②, Eq. (20) corresponds to the Exner equation in the two-phase flow formalism:

$$ q(x, y, t) = \int_{z_b(x, y, t)}^{h} \phi \mathbf{u}^s \, dz = \underbrace{\int_{z_b(x, y, t)}^{z_s(x, y, t)} \phi \mathbf{u}^s \, dz}_{q_{bl}(x, y, t)} + \underbrace{\int_{z_s(x, y, t)}^{h} \phi \mathbf{u}^s \, dz}_{q_{sl}(x, y, t)} \quad (19) $$

$$ \iint_{S_b} \phi \frac{\partial z_b}{\partial t} dS_b + \frac{d\Psi}{dt} + \oint_l (q_{bl} + q_{sl}) \cdot n_l dl + \iint_{S_b} \phi \mathbf{u}^s \cdot n_b dS_b + \iint_{S_l} \phi \mathbf{u}^s \cdot z dS_t = 0 \quad (20) $$

Fig. 2 depicts the different terms according to a control volume V. The three dimensional arrows represent the sediment flux vector. The right panel of Fig. 3 gives the classical conceptual sketch of the Eulerian–Eulerian two-phase flow approach and mainly differs with the single-phase flow conceptual sketch in the left panel of the same figure due to resolution of the bed under layer.

## 2.3. Models comparison: Key similarities and differences

Following this brief introduction to the single- and two-phase flow approaches, Eq. (21) summarizes the various forms that the Exner equation can take. The first applies to single-phase flow models, while the last applies to two-phase flow models when considering a control volume comprising the bed-load and suspended-load layers, i.e. between $z_b$ and $h$. Note that in the two-phase flow formalism, the additional terms represent the flux crossing the top surface, which is often negligible depending on the case and the boundary conditions set. For the sake of clarity the top flux in the two-phase flow approach will be neglected in the section. The second additional flux is the one crossing the bed interface and exchanging mass with the fixed porous media which is unresolved in the single-phase flow model.

$$ \left\{ \begin{array}{l} \iint_{S_b} \phi_b \frac{\partial z_b}{\partial t} dS_b + \frac{d\Psi}{dt} + \oint_l (q_{bl} + q_{sl}) \cdot n_l dl = 0 \\ \iint_{S_b} \phi \frac{\partial z_b}{\partial t} dS_b + \frac{d\Psi}{dt} + \oint_l (q_{bl} + q_{sl}) \cdot n_l dl + \iint_{S_b} \phi \mathbf{u}^s \cdot n_b dS_b = 0 \end{array} \right. \quad (21) $$

The term “bed elevation” in the single-phase flow approach does not exactly match its counterpart in the two-phase flow approach. While the solid volume fraction of the bed is fixed to $1 - p = \phi_b$ in the former, it is evaluated at $z_b$ in the latter (see Section 2.4). In the single-phase approach, the bed-load layer is not explicitly resolved. This means that the bed elevation, $z_b$, does not directly correspond to the level at which bed-load transport occurs. To differentiate the notations, the variable $\tilde{z}_b$ is introduced which represents the bed surface corresponding to the top of the bed-load layer in the single-phase flow approach. Since the bed-load layer is not resolved, its physical thickness $\delta_b$ may be neglected and $z_b \approx \tilde{z}_b$ is assumed. As the separation between the bed-load and the suspended-load layers is quite arbitrary in the single-phase flow approach, it is difficult, if not impossible, to provide a clear definition. In this manuscript, we suggest to define it as follows: $\tilde{\delta}_b = z_{ref} - z_b \approx z_{ref} - \tilde{z}_b$. On the contrary, in the two-phase flow approach, this separation is not required and the two layers are defined a posteriori based on a concentration criteria as follows, the bed-load layer corresponds to the region between two iso-concentration $\phi_b$ and $\phi_s$; $z_b = \{z \in [-\infty; h] | \phi = \phi_b\}$ and $z_s = \{z \in [-\infty; h] | \phi = \phi_s\}$, corresponding to the bed and the suspension interfaces; the suspension layer corresponds to the region between the iso-concentration $\phi_s$ and the free surface at $z = h$. In this framework, the bed-load flux is computed between $z_b$ and $z_s$, leading to the definition: $\delta_b = z_s - z_b$. The different notations introduced herein are illustrated in Fig. 3.

Both approaches yield similar velocity and concentration profiles within the suspension layer. The major difference lies in the existence of a velocity specific to the solid phase in the two-phase approach, whereas in the single-phase approach the sediments are transported at the velocity of the fluid flow, except in the vertical direction where a settling velocity is superimposed.

The bed-load layer thickness $\delta_b$ is on the order of a few grain diameters, a classical linear scaling with the Shields number $\theta$ was proposed by Wilson (1987): $\delta_b/d_{50} = 10 \times \theta$. Such scaling is recovered by the two-phase flow approach (e.g. Chauchat, 2018). In general, the bedload layer thickness $\tilde{\delta}_b$, is ignored in single-phase flow models. In some models, it is used to parameterize the bed roughness (e.g. Warner et al., 2008). The second main difference comes from the definition of the suspended-load layer, it is not obvious to claim that $z_s$ in the two-phase flow approach is equivalent to the reference elevation $z_{ref}$. Furthermore, in the single-phase flow approach, it is not entirely clear how $z_{ref}$ is defined and where the net flux, $\pm(E - D)$, shall be imposed. To the best of the authors knowledge, there is no publications on the net flux using the two-phase flow approach. This is a strong motivation to apply the Exner-vision between $z_b$ and $h$ in the two-phase flow approach and not between $z_s$ and $h$.

Nevertheless, it is possible to express the net flux in the two-phase flow approach by subtracting the mass balance between $z = z_b$ and $z = h$ and the mass balance between $z = z_s$ and $z = h$:

$$ \iint_{S_b} (E - D) dS_b \approx - \iint_{S_b(t)} \phi \mathbf{u}^s \cdot n_s dS_s - \iint_{S_s} \phi \frac{\partial z_s}{\partial t} dS_s, \quad (22) $$

where the net sediment flux consists of two contributions: the sediment flux crossing the suspension interface $S_s$, which represents the entrainment of sediments from the bed-load layer into the suspension layer and the temporal evolution of the suspension level. This is exactly how van Rijn, defined the net flux in his seminal work in 1984 (Van Rijn, 1984a). However it is not possible at this point to make a comparison with the single-phase flow approach as the definition of the suspended-load layer is arbitrary.

## 2.4. Definition

In this study, the bed elevation $z_b(x, y, t)$ is defined as the top interface of the cell having a solid volume fraction $\phi > \phi_b$. It can be visualized in red in the example given in Fig. 4(a) with $\phi_b = 0.59$. The threshold value will be analyzed for different examples. This arbitrary
<figure><img src="image_5.png" /><figcaption>Fig. 3. Representation of the single-phase (left panel) and two-phase (right panel) flow models concept.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Representation of the two-phase flow model concept with discrete approach.</figcaption></figure>

solid volume fraction is however chosen with a value larger than the random loose packing ($\phi_{min}^{Fric} \approx 0.57$) so that particles form a contact network responsible for dissipating the particle phase kinetic energy by friction.

The particle phase velocity is known inside the cell, however it is not known at the interface $z_b$ defined previously. To establish the velocity at the interface, a linear interpolation value is used from the two neighbor cells. It can be visualized with text color in cyan in the example given in Fig. 4(a). The same method is used to estimate the value of the velocity at the lateral interface of the control volume as can be seen with text color in green again in Fig. 4(a).

Fig. 4(b) represents the different components of the lateral flux and bed flux term when considering three dimensional configuration. As shown in the figure, the number of components of the lateral flux depends of the number of lateral faces, here four in this sketch. Concerning the bed flux term, it exists three components depending on the orientation of the faces representing the bed interface. These components are investigated in Section 3.3.

### 3. Results and discussions

Mass balance conservation is studied in three different configurations with an increasing complexity, first a one dimensional pure sedimentation case, second a two dimensional scour downstream an apron case and finally a three dimensional scour case. The numerical simulations have been performed using *sedFOAM*, an open-source two-phase flow Eulerian–Eulerian model (Chauchat et al., 2017; Bonamy et al., 2024) implemented in the CFD toolbox *openFOAM* and the post processing to apply Eq. (18) using Python with help of *fluidFOAM* package (Bonamy et al., 2025). The two-phase flow model has been successfully validated on numerous numerical test cases, notably in the estimation of the bed-load sediment flux (e.g., Chauchat, 2018; Lee et al., 2016; Hsu et al., 2004). It also has been applied on many configurations including morphodynamic processes such as sedimentation and scour (e.g., Nagel et al., 2020; Tsai et al., 2022).
<figure><img src="image_7.png" /><figcaption>Fig. 5. Sketch of the time evolution of the numerical experiment one dimensional pure sedimentation.</figcaption></figure>

### 3.1. One dimensional pure sedimentation

The first test case to validate this methodology is a pure one dimensional sedimentation case. The numerical simulation reproduces experiment presented in Pham Van Bang et al. (2008). It consists of the sedimentation of a suspension of mono-disperse spherical polystyrene beads having diameter $d_p = 0.25 \pm 0.03$ mm, and density $\rho_s = 1.05$ kg m⁻³ in Rhodorsil silicone oil having a viscosity $\eta_f = 20 \ 10^{-3}$ Pa s, density $\rho_f = 0.95$ kg m⁻³. This configuration has been successfully validated with sedFOAM in Chauchat et al. (2017).

A uniform grid made of 1200 cells ($\Delta z = 0.2 d_p = 5 \ 10^{-5}$ m) and a simulation time step of $\Delta t = 10^{-2}$ s have been used. The initial condition is a constant solid volume fraction equal to $\phi_0 = 0.5$ below $z = 0.06$ m and gradual decreasing to zero above. The fluid and particle phase vertical velocities are initially equal to zero. A Dirichlet boundary conditions set to zero at the top and bottom boundaries is implemented for the vertical velocities and Neumann boundary conditions taken to zero for the solid volume fraction. For the fluid pressure, a zero value is imposed at the top boundary while a simplified momentum balance is used at the bottom boundary (fixedFluxPressure). Fig. 5 describes schematically the numerical experiment.

In this case, there is no lateral flux term ⑆ and the top flux term ⑇ is negligible. Consequently, the mass balance equation (18) is then reduced to ⑆ + ⑇ + ⑉ = 0, which reads as follows in discrete form:

$$ \phi(z_b, t) \frac{\Delta z_b(t)}{\Delta t} \Delta x \Delta y + \frac{\Delta \Psi(t)}{\Delta t} - \phi(z_b, t) u_z^x(z_b, t) \Delta x \Delta y = 0, \quad (23) $$

where the notation $\Delta z_b(t)$ stands for $z_b(t + \Delta t) - z_b(t)$, $\Delta t = 20$ s is the simulation output time step, $\Delta x = \Delta y = 0.0006$ m are the constant cell length. In this discrete formulation it is assumed that the temporal derivative of term ⑆ and ⑇ is accurately estimated by an Euler forward scheme defined as $\frac{\Delta x}{\Delta t} \int_{t}^{t+\Delta t} (x_i^x - x_i^y)$ representing the derivative of $x$ at the $i$th cell at the $n$th output simulation time step.

Fig. 6 shows the time evolution of the residual defined as ⑆ + ⑇ + ⑉ plotted for different value of $\phi_b$. It illustrates the level of accuracy depending on the iso-concentration chosen, the greater $\phi_b$ the smaller the residual. Important fluctuations are observed for $\phi_b = 0.57$ while it is significantly smaller for larger values of $\phi_b$. Looking at the inset where the relative residual is compared with the mean value of the bed

<figure><img src="image_8.png" /><figcaption>Fig. 6. Time evolution of the residual of the mass balance equation in the one dimensional pure sedimentation case for different iso-concentration criteria. The inset depicts the relative residual of the average bed evolution: $Res = (\textcircled{1} + \textcircled{2} + \textcircled{3}) / \frac{\int_{t_{\text{final}}}^{t_{\text{initial}}} \textcircled{1} dt}{t_{\text{final}}}$.</figcaption></figure>

evolution term ⑆ of each criteria, the numerical error is smaller than 10% for $\phi_b = 0.57$ which is still reasonable. The fluctuations possibly originate from the error made on temporal derivative estimations that assume slow variation of each term between every output time step $\Delta t$. This is less true the case when particles are loose.

According to the previous figure, the choice of the criteria $\phi_b$ must be decided based on the analysis of the numerical error while maintaining sufficiently small iso-concentration so that a bed evolution occurs. Indeed, a too large iso-concentration for the bed interface would mean carrying out a mass balance on the entire domain and losing information on the most important term of the Exner equation $\frac{\partial z_b}{\partial t}$. Fig. 7 shows results from iso-concentration criteria $\phi_b = 0.59$ chosen as the best compromise for this configuration.

Fig. 7(a) represents the equilibrium between the bed evolution term ⑆ and the other terms, namely - (② + ④) in this case. As the symbols are on the straight line $y = x$, the mass balance is satisfied. They are colored from cold to warm representing the time evolution of the mass balance. At the beginning of the simulation, the symbols are located in the top right-hand corner meaning that the bed evolution is positive. Physically, sedimentation occurs and the bed elevation increases. At the end, the symbols get closer to zero meaning that the bed reaches its equilibrium position and the sedimentation process is complete.

Fig. 7(b) shows the time evolution of each term of Eq. (23). The bed evolution term ⑆ is initially equal to zero because no bed is formed yet with a solid volume fraction greater than $\phi_b = 0.59$. From $t \approx 120$ s, a bed exists and deposition occurs responsible for a rise of the bed and implies a positive bed evolution. It reaches a stable value meaning that the bed is linearly increasing with time until $t \approx 1400$ s when the bed evolution returns to zero as the sedimentation process is finished. The storage evolution term ⑇ evolves as the opposite of the bed evolution. As the sedimentation process takes place, the suspended solid particles are settling meaning that the number of solid particles in suspension decreases with time explaining the negative value of this term. Clearly, the storage evolution counterbalances the bed evolution term ⑆ and is associated with the main process that drives the bed evolution. The bed flux term ⑉ is slightly positive meaning that the sediment flux goes out from the bottom surface. This is not the principal term that drives the mass balance in this case but is essential to perfectly close the mass balance. It demonstrates that two-phase flow model perfectly resolved the Exner equation by itself in this case.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Results of the mass balance budget for the one dimensional sedimentation case with an iso-concentration criteria $\phi_b = 0.59$.</figcaption></figure>

## 3.2. Two dimensional scour downstream an apron

In order to demonstrate the capability of the method to deal with two dimensional situations, the second test case is the scour downstream an apron. Numerous studies (e.g., Amoudry and Liu, 2009; Cheng et al., 2017; Nagel, 2018) investigated this configuration and since the purpose of this study is to show the capability of the mass balance methodology, no comparison with previous results are shown. The sand has the following characteristics in the numerical experiments, a solid density $\rho_s = 2.65 \ 10^3 \text{ kg m}^{-3}$ (sand) with a particle diameter $d_{50} = 0.25 \text{ mm}$. The fluid is water with density $\rho_f = 1.0 \ 10^3 \text{ kg m}^{-3}$ and kinematic viscosity $\nu_f = 10^{-6} \text{ m}^2 \text{ s}^{-1}$.

The water depth is $h = 0.15 \text{ m}$, and the initial bed is $0.05 \text{ m}$ thick with a concentration profile initialized with a 1D simulation. The domain is 1 m long. A uniform grid is used in the streamwise direction $\Delta x = 10^{-3} \text{ m} = 4d_{50}$ while the mesh is refined vertically near the bed interface: $\Delta z \in [1.1 \ 10^{-4} \text{ m}, 1.1 \ 10^{-2} \text{ m}] = [0.4d_{50}, 42d_{50}]$ in the water column and $\Delta z \in [1.1 \ 10^{-4} \text{ m}, 4.6 \ 10^{-4} \text{ m}] = [0.4d_{50}, 1.8d_{50}]$ in the sediment bed. Front and back boundaries are set to empty for the two dimensional case. The bottom boundary, the lower parts of the inlet forming the step and of the outlet are set to wall boundaries. The upper part of the inlet is an inlet boundary where the velocity profiles are imposed according to the rough-wall-log-law profile. The upper part of the outlet, Neumann boundary conditions ($\frac{\partial}{\partial x_i} = 0$) are specified for all quantities, except for the pressure for which a homogeneous Dirichlet boundary condition ($p^* = p - \rho_f g z$, namely the pressure where the hydro-static contribution is removed. For the velocities, a homogeneous Neumann boundary condition is used when the velocity vector points outside of the domain at the outlet, and a homogeneous Dirichlet boundary condition is used otherwise (inletOutlet condition in openFOAM). Because the free surface is not resolved, the top boundary is set to a Neumann boundary condition for all quantities. Concerning the numerical schemes, a second order scheme for the discretization of the divergence operators (linearUpwind in openFOAM), a second order for the discretization of the gradient operator (linear) and a first order Euler temporal scheme are used in this simulation. As initial condition, the velocity of both phases, the sediment concentration, the TKE and the TKE dissipation are set based on one-dimensional simulation results using the funkysetFields tool and routines from the fluidFOAM python package (Bonamy et al., 2025). Nagel (2018) showed a sensitivity analysis using turbulence Reynolds-Averaged Navier-Stokes and granular stresses models. His results showed a higher dependency on the former

<figure><img src="image_10.png" /><figcaption>Fig. 8. Sketch of the two dimensional scour downstream an apron numerical experiment. The red dashed volume represents the control volume where the mass balance is applied.</figcaption></figure>

than the latter on the scour depth. The simulation performed herein uses a RANS $k-\omega$2006 turbulence model (Nagel et al., 2020) and the $\mu(I)$ rheology for the granular stresses model. Fig. 8 shows a sketch of the problem.

In this case, the lateral flux term ③ is non-zero and the mass balance Eq. (18) reduces to ①+②+③+④ = 0. The mass balance is applied to a control volume delimited by the red dashed lines in Fig. 8. This control volume is defined by fixed abscissa values: $x = 0.05 \text{ m}$ and $x = 0.15 \text{ m}$ and the water column defined between $z_b$ and $h$. The mass balance equation reads as follows in discrete form:

$$ \begin{aligned} & \sum_{i=l_{min}}^{l_{max}} \phi(x_i, z_b(x_i, t), t) \frac{\Delta z_b(x_i, t)}{\Delta t} \Delta x \Delta y + \frac{\Delta \Psi(t)}{\Delta t} \\ & - \sum_{k=b(l_{min})}^{n_z} \phi u_x^s(x_{l_{min}} - \frac{1}{2}, z_k, t) \Delta y \Delta z_k \\ & + \sum_{k=b(l_{max})}^{n_z} \phi u_x^s(x_{l_{max}} + \frac{1}{2}, z_k, t) \Delta y \Delta z_k \\ & - \sum_{i=l_{min}}^{l_{max}} \phi u_z^s(x_i, z_b(x_i, t), t) \Delta x \Delta y \\ & + \sum_{i=l_{min}}^{l_{max}-1} \sum_{k=b(i)}^{b(i+1)} \phi u_x^s(x_{i+1/2}, z_k, t) \Delta y (z_b(x_{i+1}, t) - z_b(x_i, t)) = 0, \quad (24) \end{aligned} $$
<figure><img src="image_11.png" /><figcaption>(a) Time evolution of the residual of the mass balance (b) Time evolution of vertical bed flux $\phi u_z^s(z_b)$ in equation. The inset depicts the relative residual of the solid lines and the horizontal bed flux $\phi u_x^s(z_b)$ in bed evolution : $Res = (①+②+③+④)/①$ dashed lines</figcaption></figure>

<caption>**Fig. 9.** Analysis of the iso-concentration criteria for the two dimensional configuration.</caption>

where the notation $\Delta z_b(x_i, t)$ and $\Delta \Psi(t)$ stands for $z_b(x_i, t + \Delta t) - z_b(x_i, t)$ and $\Psi(t + \Delta t) - \Psi(t)$, $i$ is the index of x axis, $b$ ($i$) is the index of the first cell above the bed interface at the $i$th cell (see Appendix B for details). The indices $i_{min-\frac{1}{2}}$ and $i_{max+\frac{1}{2}}$ corresponds to the interface at $x = 0.05$ and $x = 0.15$ respectively. The flux at this interface is estimated with a linear interpolation illustrated in Fig. 4(a). Note that the fifth term is denoted as the vertical bed flux and the sixth term as the longitudinal bed flux.

Fig. 9(a) depicts the time evolution of the mass balance residual defined as ①+②+③+④ plotted for different value of $\phi_b$. Similarly as in the one dimensional case, the larger $\phi_b$ the smaller the residual. The error is particularly important for iso-concentration $\phi_b = 0.57$ and $\phi_b = 0.58$. It is correlated with the vertical bed flux term illustrated in Fig. 9(b). The mass balance is completed for $\phi_b \ge 0.60$, when the vertical bed flux remains negligible.

As for the 1D sedimentation case, the results presented in Fig. 9(a) confirm that the bed criteria, $\phi_b$, must be large enough for the mass balance diagnostic to be closed. In order to better understand the reason for this issue Fig. 10 shows the bed flux longitudinal profile in the control volume. It is clearly seen that for the lowest values of $\phi_b$ (0.57; 0.58; 0.59) the bed flux is exhibiting significant longitudinal oscillations which are due to the fact that the volume fraction at the cell position (see Fig. B.19) is changing significantly from one cell to the other in the streamwise direction. These significant variations are due to the strong vertical gradients of the volume fraction profile at volume fractions below 0.6. The position $z_b$ deduced from a $\phi_b$ criterion that is too low, fluctuates greatly in space. This induces the observed oscillations, which correspond to the integral of $u_z^s\phi$ at $z = z_b$ over the bed surface $S_b$. These oscillations are due to the discretization in z direction and therefore to an approximation error in the position of $z_b$, which is linked to the $\phi_b$ criterion, and not to a numerical instability.

According to the previous sensitivity analysis to the iso-concentration criteria defining the bed $\phi_b$, in the following, the results will be shown for $\phi_b = 0.60$. Fig. 11 shows the instantaneous snapshots of the simulation. In the control volume region dashed red lines, the bed evolution is slowly decreasing in time as expected. The simulation will be further analyzed using the mass balance equation formulated for two-phase flow simulation.

Fig. 12(a) depicts similarly to Fig. 7(a) the equilibrium between the bed evolution term ① and the remaining terms of the mass balance, namely -(②+③+④). As the symbols are on the straight line $y = x$, the

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Bed vertical velocity profile inside the control volume considered for different iso-concentration criteria at $t = 50$ s.</figcaption></figure>

equilibrium is satisfied. All points are located on the bottom left-hand corner meaning that the bed evolution is negative inside the control volume and that an erosion occurs. At the beginning of the simulation, the symbols are located far from the origin before getting closer to it at the end. Indeed, the scour process is fast at the beginning before slowing down and stop.

Fig. 12(b) represents similarly to Fig. 7(b) the time evolution of each terms of this configuration inside a control volume delimited by the red dashed lines in Fig. 8. Looking at the bed evolution term ①, as previously observed in Fig. 12(a), it remains negative during all the simulation meaning that erosion takes place. The scour process is maximum at the early stages before decreasing gradually towards a stable negative value meaning that the scour is still active at the end of the simulation. The storage evolution term ② begins with a strong negative peak explained by a quick decrease of the solid volume fraction in the water column. This is due to the initialization that fixed an important solid volume fraction into the suspension, which was quickly washed away from the control volume. This also explains the
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Instantaneous snapshots of the two dimensional scour case with an iso-concentration criteria $\phi_b = 0.60$. The red dashed line represents the control volume. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

importance of the lateral flux term ③. Next, the storage term oscillates around zero with fluctuations with the bed evolution term ①. The oscillations observed in Fig. 12(b) differ from those in Fig. 9. We are convinced that these oscillations have a physical origin and are linked to the intermittent erosion/deposition process between the fixed bed (below $z_b$) and the storage term. This intermittent behavior is likely due to turbulent flow and the existence of a recirculation cell downstream of the apron. While the use of a URANS model does not guarantee that these oscillations are physically realistic, it is reasonable to assume that the recirculation cell intensity oscillates in time and that these oscillations influence the erosion process. The small residual of the mass balance tends to confirm that these oscillations are physical and not due to a numerical issue.

### 3.3. Three dimensional scour around a square cylinder

Finally, the methodology proposed in this manuscript is tested on a three dimensional case of scour around a vertical square cylinder. The third dimension adds an other component on the spanwise direction for terms ③ and ④. The numerical configuration is shown in Fig. 13, a channel filled with a sand layer of 0.05 m thick with a median diameter $d_{50} = 260$ µm and a solid phase density $\rho_s = 2650$ kg m$^{-3}$ and a water layer of $h = 0.2$ m. A square cylinder with pile diameter $D = 0.1$ m is buried and fixed in the sand layer. The configuration is live-bed meaning that sediment is supplied at the channel inlet with a Shields number of $\theta_0 = \frac{\tau_b}{(\rho_s-\rho_f)gd_{50}} = 0.19$. A numerical simulation is performed with the following computational domain: 12D long, 8D wide and 2.5D high. The center of the cylinder is placed in the symmetry plane of the computational domain at 4D from the inlet. Fig. 13 shows the location

of the control volume in red dashed line where the mass balance will be assessed. In terms of dimensionless numbers, two Reynolds numbers can be defined, the bulk Reynolds number $Re_b = \frac{\bar{u}h}{v_f} = 9.2 \ 10^4$ and the pile Reynolds number $Re_D = \frac{\bar{u}D}{v_f} = 4.6 \ 10^4$, with $h$ the water depth, $\bar{u} = 0.46$ m s$^{-1}$ the depth-averaged velocity defined as $\bar{u} = \frac{1}{h-z_b} \int_{z_b}^{h} u^f dz$ where $z_b = z(\phi_b)$ is the bed elevation computed with a solid volume fraction $\phi_b$ criteria.

Periodic boundary conditions are used on the sides whereas wall boundary conditions are applied on the cylinder and the bottom. A steady RANS velocity profile performed previously using a 1D simulation is imposed at the inlet. As in the 2D case, Neumann boundary conditions are specified at the outlet for all quantities except for the pressure for which a homogeneous Dirichlet boundary condition is imposed for the reduced pressure and for the velocities, the inletOutlet boundary condition is used. As the free surface is not resolved, Neumann boundary condition is set for all quantities at the top. At the walls (including the cylinder), zero velocity (no-slip) is imposed for the three components. The boundary condition for the turbulent kinetic energy and the dissipation frequency are specified using the classical wall function from openFOAM at the cylinder, using a Neumann boundary condition and constant value respectively at the bottom and conversely around the pit. Concerning the numerical schemes, a second order Total Variation Diminishing scheme for the discretization of the divergence operators (limitedLinear in openFOAM), a second order for the discretization of the gradient operator (linear) and a second order temporal scheme (backward) are used in this simulation. The simulation performed herein uses a RANS $k-\omega$2006 turbulence model as in Nagel et al. (2020) and the kinetic theory of granular flows
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Results of the mass balance budget for the two dimensional scour case with an iso-concentration criteria $\phi_b = 0.60$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Sketch of the three dimensional scour numerical experiment.</figcaption></figure>

for the granular stresses model (Chassagne et al., 2023). The simulation is initialized by mapping the fields (mapFields function in openFOAM) from a previous simulation having the same configuration lasting 5 s where the sediment phase below $z = 0$ is freed (initFreeze option in sedFOAM). This enables the flow to develop especially around the cylinder without morphological evolution.

Considering the same solid phase parameters as the two dimensional case, the iso-concentration criteria is fixed to $\phi_b = 0.60$. A sensitivity analysis to the grid is investigated in this configuration. Indeed, because of the additional third dimension, the number of cells increases drastically such that the mesh refinement in the horizontal directions is 1 to 16 times coarser than in the previous two-dimensional case. The mesh is refined close to the cylinder and Table 1 gives the range of cell size in the streamwise, spanwise and vertical direction, inside the control volume defined in the red dashed box in Fig. 13. This control volume is defined by fixed values of $x_{min} = -2D$, $x_{max} = -0.5D$, $y_{min} = -0.5D$, $y_{max} = 0.5D$ and the water column defined between $z_b(t)$ and $h$. Due to different meshes, the discrete control volume is delimited by the first cell centers located inside the limits defined previously.

Similarly to the two-dimensional case the mass balance equation (18) is reduced to $① + ② + ③ + ④ = 0$, the top flux is also negligible. The mass balance equation reads the following discrete form:

$$ \sum_{i=\min}^{l_{max}} \sum_{j=\min}^{J_{max}} \phi(x_i, y_j, z_b(x_i, y_j, t), t) \frac{\Delta z_b(x_i, y_j, t)}{\Delta t} \Delta x_i \Delta y_j + \frac{\Delta \psi(t)}{\Delta t} $$

<table><caption>Table 1
Mesh refinement for the scour simulation inside the control volume.</caption><thead><tr><th rowspan="3">Mesh name</th><th colspan="4">Mesh refinement (min-max)</th><th rowspan="3">Cells number</th><th rowspan="3">Δt (s)</th></tr><tr><th colspan="2">Δx (d<sub>50</sub>)</th><th colspan="2">Δy (d<sub>50</sub>)</th></tr><tr><th colspan="2">Bed (d<sub>50</sub>)</th><th colspan="2">Water (d<sub>50</sub>)</th></tr></thead><tbody><tr><td>M1</td><td>64–112</td><td>43</td><td>0.8–1.4</td><td>1.3–55</td><td>3.1 10<sup>5</sup></td><td>5 10<sup>-4</sup></td></tr><tr><td>M2</td><td>32–55</td><td>32</td><td>0.8–1.4</td><td>1.3–55</td><td>9.5 10<sup>5</sup></td><td>10<sup>-4</sup></td></tr><tr><td>M3</td><td>16–28</td><td>16</td><td>0.8–1.4</td><td>1.3–55</td><td>3.8 10<sup>6</sup></td><td>10<sup>-4</sup></td></tr><tr><td>M4</td><td>8–14</td><td>8</td><td>0.8–1.4</td><td>1.3–55</td><td>6.6 10<sup>6</sup></td><td>10<sup>-4</sup></td></tr><tr><td>M5</td><td>5–9</td><td>5</td><td>0.8–1.4</td><td>1.3–55</td><td>9.8 10<sup>6</sup></td><td>10<sup>-4</sup></td></tr></tbody></table>

$$ \begin{aligned} & -\sum_{j=j_{min}}^{j_{max}} \sum_{k=b(i_{jmin}, j)}^{n_z} \phi u_x^s (x_{i_{min}-\frac{1}{2}}, y_j, z_k, t) \Delta y_j \Delta z_k \\ & -\sum_{i=\min}^{i_{max}} \sum_{k=b(i, j_{min})}^{n_z} \phi u_y^s (x_i, y_{j_{min}-\frac{1}{2}}, z_k, t) \Delta x_i \Delta z_k \\ & +\sum_{i=\min}^{i_{max}} \sum_{k=b(i, j_{max})}^{n_z} \phi u_y^s (x_i, y_{j_{max}+\frac{1}{2}}, z_k, t) \Delta x_i \Delta z_k \\ & -\sum_{i=\min}^{i_{max}} \sum_{j=j_{min}}^{j_{max}} \phi u_z^s (x_i, y_j, z_b(x_i, y_j, t), t) \Delta x_i \Delta y_j \\ & +\sum_{j=j_{min}}^{j_{max}} \sum_{i=\min}^{i_{max}-1} \sum_{k=b(i, j)}^{b(i+1, j)} \phi u_x^s (x_{i+\frac{1}{2}}, y_j, z_k, t) \Delta y_j (z_b(x_{i+1}, y_j, t) - z_b(x_i, y_j, t)) \\ & +\sum_{i=\min}^{i_{max}} \sum_{j=j_{min}}^{j_{max}-1} \sum_{k=b(i, j)}^{b(i+1, j)} \phi u_y^s (x_i, y_{j+\frac{1}{2}}, z_k, t) \Delta x_i (z_b(x_i, y_{j+1}, t) \\ & \qquad -z_b(x_i, y_j, t)) = 0, \end{aligned} \quad (25) $$

where $b(i, j)$ is the indices of the first cell above the bed interface at the $i^{th}$, $j^{th}$ cell. The indices $i_{min}$ and $i_{max}$ corresponds with the closest interface between cells at $x = -2D$ and $x = -0.5D$ respectively and $j_{min}$ and $j_{max}$ with the closest interface at $y = -0.5D$ and $y = 0.5D$ respectively. Note that the sixth term is denoted as the vertical bed flux, the seventh term as the longitudinal bed flux and the eighth term as the spanwise bed flux, all together corresponding to the bed flux term ④.

Fig. 14(a) depicts the equilibrium between the bed evolution term ① and the other terms, namely -((②+③+④)) for five different meshes. Only negative range is shown on this figure because the control volume is subjected to the scour process. Apart from the first second of simulation time, the symbols are located in the farthest part from the origin at the beginning before getting closer to it and reach an equilibrium stage
<figure><img src="image_17.png" /><figcaption>Fig. 14. Results of the mass balance budget for the three dimensional scour case with an iso-concentration criteria $\phi_b$ = 0.60.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Instantaneous snapshots of the three dimensional scour case with an iso-concentration criteria $\phi_b$ = 0.60. The cloud represents the suspended-particles and the red line represents the control volume.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Instantaneous scour pattern with an iso-concentration criteria $\phi_b = 0.60$ for mesh M4 at $t = 20$ s.</figcaption></figure>

where the bed evolution remains constant. While the coarse mesh M1 (cross 'x' points) fails in satisfying the mass balance, the results start to be reasonable for mesh M4 as the symbols get closer to the straight line $y = x$ demonstrating the closure of the mass balance diagnostic. Looking at the residuals evolution in time shown in Fig. 14(b), the finer mesh, the smaller the residual. The results converge to a reasonable residual for mesh M4, however the error keeps increasing from $t = 10$ s. In terms of cell length criteria, it can be concluded that the mass balance two-phase flow formulation is satisfied for a horizontal grid spacing smaller than to $14d_{50}$.

According to the sensitivity analysis to the grid resolution presented above, in the following the results of the mass balance equation for mesh M4 are shown. Fig. 15 presents different instantaneous snapshots of the simulation. In the control volume, shown as in red lines, the bed evolution is decreasing in time as expected. The particles are suspended very early in the simulation in the wake and at a lesser extend inside the scour hole. Fig. 16 presents snapshots at time $t = 20$ s of the bed elevation $z_b$ defined as an iso-concentration $\phi = 0.60$ made dimensionless by the pile length $D$. As a result, the scour depth reach the computational domain bottom close to the pile at the upstream pile. At the downstream side, no scour are developing in the vicinity of the pile and two dunes are slowly migrating. The scour mark is perfectly symmetric in this configuration. The simulation process is going to be further analyzed in the next paragraphs thanks to the mass balance equation formulated for two-phase flow simulation.

Fig. 17(a) represents the time evolution of each term inside the control volume delimited by the red box upstream the pile shown in Fig. 13. The bed evolution term ① is negative corresponding to an erosion due to the scour process upstream the pile. Note that this evolution tends to slow down at $t = 10$ s due to the complete erosion of the initial sediment bed. The storage evolution term ② is slightly positive during the first 10 s suggesting an extraction of particles from the bed that are put into suspension leading to an increase of the solid volume fraction in the water column. In a second stage, it evolves around zero with opposite fluctuations of the bed evolution term ① suggesting a similar mechanism as the two-dimensional scour case. The lateral flux term ③ is positive all the simulation time meaning that more sediments are going out through the lateral surface closing the control volume. Fig. 17(b) describes the amount of flux crossing each lateral surfaces, while less sediments are going in through the inflow term ③a at $x = -2D$. The sediments are mainly going out equally through the outflow term ③b and ③c located on the sides of the control volume. It largely counterbalances the bed evolution during all the

simulation. Finally the bed flux term ④ is almost equal to zero. Giving a closer look to each term in Fig. 17(b), the vertical component ④b is positive suggesting that sediments are going out through the horizontal bed interface counterbalanced by the longitudinal component term ④a which is negative corresponding to sediments that are coming in through the bed interface. The spanwise component term ④c is negligible here due to the symmetry of the configuration. Overall, a zero bed flux can be associated to the sediment flux parallel to the bed slope.

The methodology has been successfully applied on a three dimensional case of scour process. Consequently, it can be used as a starting point to analyze the scour process in different configurations.

A key point for further analysis is the parameterization of sediment flux which drives morphological evolution. In the single-phase flow approach, the instantaneous sediment flux is typically modeled using empirical formulations (Meyer-Peter and Müller, 1948; Wong and Parker, 2006), which are valid under steady, uniform flow conditions. As a result, applying these formulations within the Exner equation to predict instantaneous and local bed evolution is questionable, since they are only reliable at larger spatial and temporal scales (Guevara et al., 2024). These relationships have been recovered in Eulerian-Lagrangian simulations using coarse-graining methods to obtain spatially averaged data (Guan et al., 2021); however, time-averaging appears to be less effective. Building on these post-processing tools, a key objective is to investigate, within an Eulerian-Eulerian framework, the validity of the empirical formulations used in the Exner equation (e.g., bedload flux, net flux) across different spatial and temporal scales.

## 4. Conclusion

A new methodology for the mass balance diagnostics has been introduced to analyze two-phase flow simulations with morphological evolution. This study targets the two-phase flow community and more specifically the multiphase sediment transport community, it is aiming to establish links between single- and two-phase flow approaches by computing the bed evolution. While the former must solve the century-old Exner equation using empirical parameterization for sediment fluxes, the latter is based on the solid continuity equation and momentum equations. The mass balance budget formulation in the two-phase flow Eulerian-Eulerian approach has been introduced, showing the links with the Exner equation. This mathematical demonstration can be assimilated to a mass balance budget performed within a time-dependent control volume.

Three different cases have demonstrated that applying the proposed formulation to numerical simulation results satisfies the mass balance budget, provided the bed concentration criterion is sufficiently high. A low iso-concentration criterion can generate significant numerical errors in the estimation of the bed flux term ④, while a coarse mesh (notably in the horizontal plane ($\Delta x, \Delta y > 14d_{50}$)) does not allow the mass balance to be closed. If these two criteria are met — namely, a large bed concentration and a fine enough grid — the mass balance diagnostic can be used to calculate the contribution of each term to erosion and deposition processes: bed evolution (term ①), storage evolution (term ②), lateral flux (term ③), bed flux (term ④) and top flux (term ⑤). In general, term ① is key in the Exner equation, representing the vertical bed interface velocity. In the sedimentation case, this term is largely driven by the storage evolution term ②, and by the lateral flux term ③ for scour upstream of a vertical square pile. Scour downstream of an apron depicts a combination of lateral flux (term ③) and bed flux (term ④). Depending on the control volume, different results are, of course, expected. Ultimately, the morphological evolution predicted by the two-phase flow model can be explained at the mesoscale using the proposed diagnosis, and the numerical results can be used to perform an upscaling from two-phase to single-phase flow approaches.
<figure><img src="image_20.png" /><figcaption>(a) Time evolution of each term of the mass balance equation (25)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Time evolution of the different components of the bed flux term ④ and the lateral flux term ③</figcaption></figure>

<caption>Fig. 17. Analysis of the two three dimensional configuration with mesh M4. referred to the web version of this article.)</caption>

This study provides a solid foundation for analyzing complex two-phase flow numerical simulations of morphologically active systems, offering a wide range of perspectives. For instance, the formulation can be applied to configurations subject to significant bed evolution, such as scouring and dune migration, or to any processes involving deposition and erosion. The ultimate aim is to perform an upscaling from grain-scale processes to system-scale dynamics. The proposed diagnosis and post-processing tools will improve our understanding of sediment flux parameterization in the Exner equation. This should enable the proposal of new formulas that are better suited to the instantaneous and local evolution of the bed than the classical empirical formulas derived from uniform and steady flow configurations (e.g. Meyer-Peter and Müller, 1948). Last but not least, all the Python post-processing codes used to compute the diagnosis and generate the figures in this article are made available on github (see Data availability statement section below).

### CRediT authorship contribution statement

**Alban Gilletta:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Cyrille Bonamy:** Writing – review & editing, Validation, Supervision, Software, Investigation, Formal analysis, Data curation, Conceptualization.

**Julien Chauchat:** Writing – review & editing, Validation, Supervision, Software, Resources, Project administration, Funding acquisition, Formal analysis, Conceptualization.

### Funding

This work was supported by the ANR through the MODULLES project (ANR-10-IEED-0006-34).

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors gratefully acknowledge the support of the French project MOdelling of marine DUnes: Local and Large-scale EvolutionS in an OWF context (MODULLES) led by the Service Hydrographique et Océanique de la Marine (SHOM), the Zentrum für Marine Umweltwissenschaften (MARUM) and France Énergies Marines (FEM). Most of the computations presented in this study were performed using the GENCI infrastructure projects number A0160107567 and A0180107567. We are thankful to the developers involved in OpenFOAM.

### Appendix A. Demonstration of the bed interface velocity

Let $f: (x, y, t) \mapsto z_b$ be the function of the bed interface. The differential form of $f$ is:

$$d f = \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y} dy + \frac{\partial f}{\partial t} dt \quad (A.1)$$

The bed interface velocity at point M of coordinates $(x, y, f(x, y) = z_b)$ is $u^b = \frac{dOM}{dt} = \left(\frac{dx}{dt}, \frac{dy}{dt}, \frac{df}{dt}\right)$.

At a given time $t_0$, let $F(x, y, z) = f(x, y, t_0) - z$ be the implicit function of the bed interface. On a point $(x_0, y_0, z_0)$ belonging to the bed interface, the normal vector is oriented by $\nabla F(x_0, y_0, z_0)$. Consequently, defining the norm of a vector $v = (v_1, v_2, v_3)$ as follows: $||v|| = \sqrt{v_1^2 + v_2^2 + v_3^2}$, the bottom normal unit vector is defined as:

$$n_b = \left( \frac{\frac{\partial F}{\partial x}}{\sqrt{\left(\frac{\partial F}{\partial x}\right)^2 + \left(\frac{\partial F}{\partial y}\right)^2 + \left(\frac{\partial F}{\partial z}\right)^2}} \right) \left( \frac{\frac{\partial f}{\partial x}}{\sqrt{\left(\frac{\partial f}{\partial x}\right)^2 + \left(\frac{\partial f}{\partial y}\right)^2 + \left(\frac{\partial f}{\partial z}\right)^2}} \right) \quad (A.2)$$

The scalar product of vectors $u^b$ and $n_b$ gives:

$$u_b \cdot n_b = \left( \frac{\frac{\partial f}{\partial x}}{\sqrt{\left(\frac{\partial f}{\partial x}\right)^2 + \left(\frac{\partial f}{\partial y}\right)^2 + \left(\frac{\partial f}{\partial z}\right)^2}} \right) \left( \frac{\frac{\partial f}{\partial y}}{\sqrt{\left(\frac{\partial f}{\partial x}\right)^2 + \left(\frac{\partial f}{\partial y}\right)^2 + \left(\frac{\partial f}{\partial z}\right)^2}} \right) \quad (A.3)$$

If $dS$ bottom surface with the normal unit vector $n_b$, we have the following equality:

$$\cos(\beta) = \frac{dS_b}{dS} = -n_b \cdot z = \frac{1}{\|\nabla F\|} \quad (A.4)$$
<figure><img src="image_22.png" /><figcaption>Fig. A.18. Sketch representing the projected surface of the bed.</figcaption></figure>

with z the vertical unit vector, S<sub>b</sub> the projected surface of S along the vertical axis and β the slope (see sketch on Fig. A.18).

Finally:

$$u_b \cdot n_b dS = - \frac{\partial z_b}{\partial t} dS_b \qquad (A.5)$$

## Appendix B. Sketch of the discrete formulation of Eq. (24) for the two dimensional case

Fig. B.19 represents the different indices of the discrete formulation of mass conservation equation. The red dashed line represents the control volume considered, the black dashed line represent the bed interface everywhere else of the control volume.

<figure><img src="image_23.png" /><figcaption>Fig. B.19. Example of indices representation for the formulation of the Exner mass balance equation for the two dimensional case.</figcaption></figure>

## Data availability

The data that support the findings of this study are openly available in GitHub repository at https://github.com/SedFoam/exnerDiagnosis.git (Bonamy and A., 2025).

## References

* Amoudry, L., Liu, P.-F., 2009. Two-dimensional, two-phase granular sediment transport model with applications to scouring downstream of an apron. Coast. Eng. 56 (7), 693-702.

* Balachandar, S., 2009. A scaling analysis for point-particle approaches to turbulent multiphase flows. Int. J. Multiph. Flow 35 (9), 801-810. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2009.02.013, Special Issue: Point-Particle Model for Disperse Turbulent Flows. URL https://www.sciencedirect.com/science/article/pii/S0301932209000251.

* Bonamy, C., A., G., 2025. SedFoam/exnerDiagnosis: Release 1.0. Zenodo. https://doi.org/10.5281/zenodo.17201939.

* Bonamy, C., Chauchat, J., Clemencot, Q., Puig-Montella, E., P., H., Mathieu, A., Bernard, A., Gonçalves, G., Skorlic, M., I., N., Chassagne, R., 2025. Fluiddyn/fluidfoam: v0.2.9. Zenodo http://dx.doi.org/10.5281/zenodo.6453089.

* Bonamy, C., Chauchat, J., Hsu, T.-J., Puig-Montella, E., Chassagne, R., Mathieu, A., Cheng, Z., Nagel, T., Salimi-Tarazouj, A., Tsai, B., Higuera, P., 2024. SedFoam/sedFoam: Release 2312. Zenodo http://dx.doi.org/10.5281/zenodo.836642.

* Borland, W.M., 1960. Stream Channel Stability. United States Bureau of Reclamation.

* Celik, I., Rodi, W., 1988. Modeling suspended sediment transport in nonequilibrium situations. J. Hydraul. Eng. 114 (10), 1157-1191.

* Chassagne, R., Bonamy, C., Chauchat, J., 2023. A frictional–collisional model for bedload transport based on kinetic theory of granular flows: Discrete and continuum approaches. J. Fluid Mech. 964, A27. http://dx.doi.org/10.1017/jfm.2023.335.

* Chauchat, J., 2018. A comprehensive two-phase flow model for unidirectional sheet-flows. J. Hydraul. Res. 56 (1), 15-28. http://dx.doi.org/10.1080/00221686.2017.1289260, arXiv:https://doi.org/10.1080/00221686.2017.1289260.

* Chauchat, J., Cheng, Z., Nagel, T., Bonamy, C., Hsu, T.-J., 2017. SedFoam-2.0: A 3-D two-phase flow numerical model for sediment transport. Geosci. Model. Dev. 10 (12), 4367-4392. http://dx.doi.org/10.5194/gmd-10-4367-2017.

* Chauchat, J., Guillou, S., 2008. On turbulence closures for two-phase sediment-laden flows models. J. Geophys. Res.-Oceans 113, 20.

* Chauchat, J., Lacaze, L., 2024. Mathematical modelling of particulate flows. In: Physics of Granular Suspensions. In: CISM International Centre for Mechanical Sciences, vol. 612, Springer Nature Switzerland, pp. 13-43. http://dx.doi.org/10.1007/978-3-031-55509-1_2, URL https://hal.science/hal-04852521.

* Chauchat, J., Mazzuoli, M., 2024. Sediment transport (fluid driven flows). In: Physics of Granular Suspensions. In: CISM International Centre for Mechanical Sciences, vol. 612, Springer Nature Switzerland, pp. 101-153. http://dx.doi.org/10.1007/978-3-031-55509-1_4, URL https://hal.science/hal-04852506.

* Chen, X., Li, Y., Niu, X., Li, M., Chen, D., Yu, X., 2011. A general two-phase turbulent flow model applied to the study of sediment transport in open channels. Int. J. Multiph. Flow 37 (9), 1099-1108. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2011.05.013, URL https://www.sciencedirect.com/science/article/pii/S0301932211001200.

* Cheng, Z., Hsu, T.-J., Calantoni, J., 2017. SedFoam: A multi-dimensional Eulerian two-phase model for sediment transport and its application to momentary bed failure. Coast. Eng. 119, 32-50. http://dx.doi.org/10.1016/j.coastaleng.2016.08.007.

* Chou, Y.-J., Wu, F.-C., Shih, W.-R., 2014. Toward numerical modeling of fine particle suspension using a two-way coupled Euler-Euler model. Part 1: Theoretical formulation and implications. Int. J. Multiph. Flow 64, 35-43. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2013.12.008, URL https://www.sciencedirect.com/science/article/pii/S0301932214000032.

* Coleman, S.E., Nikora, V.I., 2009. Exner equation: A continuum approximation of a discrete granular system. Water Resour. Res. 45 (9).

* Engelund, F., Fredsoe, J., 1976. A sediment transport model for straight alluvial channels. Nord. Hydrol. 7 (5), 293-306.

* Exner, F.M., 1925. Über die wechselwirkung zwischen wasser und geschiebe in flüssen (in German). Sitz. Acad. Wiss. Wien Math. Natur- Wiss. Abt. 2a 134, 165-203.

* Finn, J.R., Li, M., 2016. Regimes of sediment-turbulence interaction and guidelines for simulating the multiphase bottom boundary layer. Int. J. Multiph. Flow 85, 278-283. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2016.06.007, URL https://www.sciencedirect.com/science/article/pii/S0301932215301282.

* Garcia, M.H. (Ed.), 2008. Sediment Transport and Morphodynamics. Sedimentation Engineering.

* Gilletta, A., Bonamy, C., Robert, M., Chauchat, J., 2025. Assessment of hybrid RANS/LES models for the prediction of the flow and scour around a wall-mounted cylinder. Adv. Water Resour. 199, 104934. http://dx.doi.org/10.1016/j.adwatres.2025.104934, URL https://www.sciencedirect.com/science/article/pii/S030917082500048X.

* Guan, L., Salinas, J., Zgheib, N., Balachandar, S., 2021. Force and torque model sensitivity and coarse graining for bedload-dominated sediment transport. Eur. J. Mech. B Fluids 90, 137-151. http://dx.doi.org/10.1016/j.euromechflu.2021.09.005, URL https://www.sciencedirect.com/science/article/pii/S0997754621001400.

* Guevara, O., Guan, L., Salinas, J.S., Zgheib, N., Balachandar, S., 2024. Sediment transport on rippled beds. Phys. Fluids 36 (11), 113339. http://dx.doi.org/10.1063/5.0236116, arXiv:https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/5.0236116/20246923/113339_1_5.0236116.pdf.

* Hsu, T., Jenkins, J.T., Liu, L.F., 2003. On two-phase sediment transport: Dilute flow. J. Geophys. Res. 108, 14.

* Hsu, T.-J., Jenkins, J.T., Liu, P.L.-F., 2004. On two-phase sediment transport: sheet flow of massive particles. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 460 (2048), 2223-2250. http://dx.doi.org/10.1098/rspa.2003.1273, arXiv: https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2003.1273. URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2003.1273.

* Jenkins, J.T., Hanes, D.M., 1998. Collisional sheet flows of sediment driven by a turbulent fluid. J. Fluid Mech. 370 (1-2), 29-52, doi: null.

* Kidanemariam, A.G., Uhlmann, M., 2014. Interface-resolved direct numerical simulation of the erosion of a sediment bed sheared by laminar channel flow. Int. J. Multiph. Flow 67, 174-188. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2014.08.008, URL https://www.sciencedirect.com/science/article/pii/S0301932214001505.
Lane, E.W., 1955. The importance of fluvial morphology in hydraulic engineering. In: Proceedings of the American Society of Civil Engineers, Journal of the Hydraulics Division. Vol. 81, pp. 45–761.

Lee, C.-H., Low, Y.M., Chiew, Y.-M., 2016. Multi-dimensional rheology-based two-phase model for sediment transport and applications to sheet flow and pipeline scour. Phys. Fluids 28 (5), 053305. http://dx.doi.org/10.1063/1.4948987, arXiv: http://dx.doi.org/10.1063/1.4948987.

Mathieu, A., Cheng, Z., Chauchat, J., Bonamy, C., Hsu, T.-J., 2022. Numerical investigation of unsteady effects in oscillatory sheet flows. J. Fluid Mech. 943, A7. http://dx.doi.org/10.1017/jfm.2022.405.

Mazzuoli, M., Balachandar, S., Chauchat, J., 2024a. Dispersed turbulent multiphase flows. In: Physics of Granular Suspensions. In: CISM International Centre for Mechanical Sciences, vol. 612, Springer Nature Switzerland, pp. 45–97. http://dx.doi.org/10.1007/978-3-031-55509-1_3, URL https://hal.science/hal-04852515.

Mazzuoli, M., Blondeaux, P., Vittori, G., 2024b. Particle-resolved direct numerical simulation of the oscillatory flow and sediment motion over a rippled bed. Int. J. Multiph. Flow 172, 104707. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2023.104707, URL https://www.sciencedirect.com/science/article/pii/S0301932223003270.

Meyer-Peter, E., Müller, R., 1948. Formulas for bed-load transport. In: Proc. 2nd Meeting. IAHR, pp. 39–64.

Montellà, E., Chauchat, J., Bonamy, C., Weij, D., Keetels, G., Hsu, T., 2023. Numerical investigation of mode failures in submerged granular columns. Flow 3, E28. http://dx.doi.org/10.1017/flo.2023.23.

Montellà, E., Chauchat, J., Chareyre, B., Bonamy, C., Hsu, T., 2021. A two-fluid model for immersed granular avalanches with dilatancy effects. J. Fluid Mech. 925, A13. http://dx.doi.org/10.1017/jfm.2021.666.

Nagel, T., 2018. Numerical Study of Multi-Scale Flow-Sediment Structure Interactions Using a Multiphase Approach (Doctoral dissertation). University of Grenoble Alpes.

Nagel, T., Chauchat, J., Bonamy, C., Liu, X., Cheng, Z., Hsu, T.-J., 2020. Three-dimensional scour simulations with a two-phase flow model. Adv. Water Resour. 138, 103544. http://dx.doi.org/10.1016/j.advwatres.2020.103544.

Ozel, A., Fede, P., Simonin, O., 2013. Development of filtered Euler–Euler two-phase model for circulating fluidised bed: High resolution simulation, formulation and a priori analyses. Int. J. Multiph. Flow 55, 43–63. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2013.04.002, URL https://www.sciencedirect.com/science/article/pii/S030193221300058X.

Paola, C., Voller, V.R., 2005. A generalized exner equation for sediment mass balance. J. Geophys. Res. Earth Surf. 110 (F4).

Pham Van Bang, D., Lefrançois, E., Sergent, P., Bertrand, F., 2008. MRI experimental and finite elements modelling of the sedimentation-consolidation of mud. Houille Blanche 3, 39–44.

Renaud, M., Bertrand, O., Bonamy, C., Chauchat, J., 2025. SedExnerFoam 2412: A 3D exner-based sediment transport and morphodynamics model. Geosci. Model. Dev. submitted for publication.

Revil-Baudard, T., Chauchat, J., 2013. A two-phase model for sheet flow regime based on dense granular flow rheology. J. Geophys. Res. Ocean. 118 (2), 619–634.

Salimi-Tarazouj, A., Hsu, T.-J., Traykovski, P., Chauchat, J., 2021a. Eulerian two-phase model reveals the importance of wave period in ripple evolution and equilibrium geometry. J. Geophys. Res. Earth Surf. 126 (7), e2021JF006132. http://dx.doi.org/10.1029/2021JF006132, e2021JF006132 2021JF006132. arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2021JF006132. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021JF006132.

Salimi-Tarazouj, A., Hsu, T.-J., Traykovski, P., Chauchat, J., 2024. Investigating wave shape effect on sediment transport over migrating ripples using an Eulerian two-phase model. Coast. Eng. 189, 104470. http://dx.doi.org/10.1016/j.coastaleng.2024.104470, URL https://www.sciencedirect.com/science/article/pii/S0378383924000188.

Salimi-Tarazouj, A., Hsu, T.-J., Traykovski, P., Cheng, Z., Chauchat, J., 2021b. A numerical study of onshore ripple migration using a Eulerian two-phase model. J. Geophys. Res. Ocean. 126 (2), e2020JC016773. http://dx.doi.org/10.1029/2020JC016773, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020JC016773.

Tsai, B., Mathieu, A., Montellà, E.P., Hsu, T.-J., Chauchat, J., 2022. An Eulerian two-phase flow model investigation on scour onset and backfill of a 2D pipeline. Eur. J. Mech. B Fluids 91, 10–26.

Van Rijn, L.C., 1984a. Sediment pickup functions. J. Hydraul. Eng. 110, 1494–1502.

Van Rijn, L.C., 1984b. Sediment transport, part II: Suspended load transport. J. Hydraul. Eng. 110, 1613–1641.

Warner, J.C., Sherwood, C.R., Signell, R.P., Harris, C.K., Arango, H.G., 2008. Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model. Comput. Geosci. 34 (10), 1284–1306.

Wilson, K.C., 1987. Analysis of bed-load motion at high shear stress. J. Hydraul. Eng. 113, 97–103.

Wong, M., Parker, G., 2006. Reanalysis and correction of bed-load relation of Meyer-Peter and Müller using their own database. J. Hydraul. Eng. 132 (11), 1159–1168.

Wu, W., Rodi, W., Wenka, T., 2000. 3D numerical modeling of flow and sediment transport in open channels. J. Hydraul. Eng. 126, 4–15.

Zhang, H., Mizutani, H., Nakagawa, H., Kawaike, K., 2015. Euler–Lagrange model for local scour and grain size variation around a spur dyke. Int. J. Multiph. Flow 68, 59–70. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2014.10.003, URL https://www.sciencedirect.com/science/article/pii/S0301932214001827.