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

# sedInterFoam 1.0: a three-phase numerical model for sediment transport applications with free surfaces

Antoine Mathieu¹, Yeulwoo Kim², Tian-Jian Hsu¹, Cyrille Bonamy³, and Julien Chauchat⁴

¹Center for Applied Coastal Research, University of Delaware, Newark, DE 19716, USA

²Department of Sustainable Engineering, Pukyong National University, Busan 48513, Republic of Korea

³LEGI, University of Grenoble Alpes, G-INP, CNRS, 38000 Grenoble, France

⁴Univ. Grenoble Alpes, INRAE, CNRS, Grenoble INP, IRD, IGE, 38000 Grenoble, France

**Correspondence:** Antoine Mathieu (amathieu@udel.edu)

Received: 28 January 2024 – Discussion started: 10 June 2024

Revised: 18 October 2024 – Accepted: 24 December 2024 – Published: 11 March 2025

**Abstract.** In this paper, an Eulerian two-phase flow model, sedFoam, is extended to include an air phase together with its water and sediment phases. The numerical model called sedInterFoam is implemented using the open-source library OpenFOAM. sedInterFoam includes the previous features of sedFoam for sediment transport modeling and also solves the air-water interface using the volume-of-fluid method coupled with the waves2Foam toolbox for free-surface wave generation and absorption. Using sedInterFoam, four test cases are successfully reproduced to validate the free-surface evolution algorithm’s implementation, mass conservation of sediment and fluid phases, and predictive capabilities and to demonstrate its potential in modeling a broader range of coastal applications with sediment transport dominated by surface waves.

## 1 Introduction

Sediment transport is the main driver of morphological changes in coastal and fluvial environments (Sherwood et al., 2022). Understanding and modeling the physical processes involved in the movement of sediment particles forced by free-surface flow is a major issue. Accounting for and integrating the complex mechanisms involved in sediment transport in numerical models, e.g., interactions between solid particles and fluid turbulence, intergranular interactions, and surface-wave-driven processes, is a critical step in improving and enhancing existing model capabilities.

In the conventional single-phase modeling approach, Navier–Stokes equations are only solved for the fluid phase, sediment transport is related to the bottom shear stress using bed-load and suspended-load formulations, and the resulting morphological evolution is calculated using the Exner equation (Jacobsen and Fredsøe, 2014; Baykal et al., 2017). This modeling methodology has several limitations, although it is widely used because of its simplicity and computational efficiency. One of the major limitations is related to the empirical relation between horizontal and vertical sediment fluxes and the bottom shear stress. Indeed, most of these relations are valid for steady or mild transient flows and relatively coarse sediments in a controlled laboratory environment. As a result, a large discrepancy is observed between predicted and measured sediment fluxes for a given wave time series (Yu et al., 2012; van der A et al., 2013), suggesting that sediment flux solely parameterized by bottom shear stress is incomplete. Another major limitation of a conventional single-phase model is its inability to simulate the onset of scour due to piping (Sumer et al., 2001), which involves gap formation and tunneling due to seepage flow driven by an upstream–downstream pressure gradient.

To overcome these limitations, the Eulerian two-phase modeling approach which encompasses most of the physical mechanisms involved in the coupling between water and particles has been actively developed in the last 2 decades (Dong and Zhang, 2002; Hsu et al., 2004; Amoudry, 2014; Lee et al., 2016; Chauchat and Guillou, 2008; Chauchat, 2018; Mathieu et al., 2021). In the Eulerian two-phase flow approach, both the carrier phase and the dispersed sediment
phase composed of the particles are seen as interpenetrating continua. Coupling between the two phases is modeled using spatially averaged interaction forces between the fluid and the particles (e.g., buoyancy, drag, or added mass), and particle-particle interactions can be included in the solid-phase stress closure models. As demonstrated in Tsai et al. (2022), an Eulerian two-phase model is able to simulate scour onset underneath a 2D pipeline due to seepage flow driven by an upstream-downstream pressure difference.

To make the scientific community benefit from the Eulerian two-phase model for sediment transport applications, an open-source solver called sedFoam (Chauchat et al., 2017; Cheng et al., 2017) implemented using the CFD library OpenFOAM (Jasak and Uroić, 2020) has been developed. The solver includes the most recent developments for intergranular stress modeling, e.g., the kinetic theory for granular flows (Chassagne et al., 2023), the $\mu(I)$ rheology (Boyer et al., 2011), turbulence closure models such as mixing length (Revil-Baudard and Chauchat, 2013), $k-\epsilon$ (Hsu et al., 2004), $k-\omega$ (Amoudry, 2014; Nagel et al., 2020), and large-eddy simulation models (Mathieu et al., 2021). sedFoam has been validated using many benchmarks (Chauchat et al., 2017) and successfully applied to various practical configurations, e.g., scour applications (Mathieu et al., 2019; Nagel et al., 2020; Tsai et al., 2022), ripple migration and geometry evolution (Salimi-Tarazouj et al., 2021a, b), sheet flow in an oscillatory boundary layer (Delisle et al., 2022; Mathieu et al., 2022), and immersed granular avalanches (Montellà et al., 2021, 2023).

Although two-phase flow models allowed us to significantly improve our knowledge of sediment transport processes in coastal and fluvial environments, this kind of modeling methodology is only applicable for configurations where free-surface effects are negligible or under simplifying assumptions, such as sediment transport in an oscillating water tunnel rather than under surface waves. Indeed, it is critical to include the ability to solve the propagation and breaking of surface waves in order to reproduce more realistic cross-shore sediment transport and morphodynamics that include the swash zone. To do so, resolving a third separated gaseous phase and implementing a numerical algorithm to resolve the evolution of the free surface is necessary. Kim et al. (2018) were the first to tackle this problem by combining an early version of sedFoam, the free-surface solver interFoam (Klostermann et al., 2013), and a library to generate free-surface waves (waves2Foam) (Jacobsen et al., 2012). With their solver, they successfully predicted sheet flow under monochromatic non-breaking waves by comparing numerical results with experimental data from Dohmen-Janssen and Hanes (2002) and studying sheet flow near a sandbar under near-breaking surface waves (Kim et al., 2019, 2021). However, the numerical implementation of the model prevented the sediment from being present in the air phase and, as a consequence, the sediment bed had to be located far away from the free surface. Despite the significant advances made

possible by adding a third gaseous phase to their model, this limitation prevented simulation of configurations for which there were strong interactions between air, water, and sediment, such as in the swash zone. More recently, Lee et al. (2019) proposed a three-phase model similar to the one presented by Kim et al. (2018) using an equivalent but more general set of equations that allows sediment in both the air and water phases.

To provide an open-source model for nearshore sediment transport applications, the purpose of this work is to present and validate the open-source solver sedInterFoam, an extension of the two-phase model sedFoam for simulating sediment transport with a free surface and for intense interactions in the swash zone. Benefiting from the earlier work of Kim et al. (2018) in combining sedFoam, interFoam, and the wave generation library waves2Foam with the set of equations proposed by Lee et al. (2019), sedInterFoam is able to simulate configurations for which sediment is present in both air and water. The model is first presented in Sect. 2. Then, the model implementation and the solution algorithm are detailed in Sect. 3. In Sect. 4, the model is validated using benchmarks and applied to configurations related to cross-shore sediment transport and beach profile evolution.

## 2 Mathematical model

In the three-phase flow solver sedInterFoam, mass and momentum conservation equations are solved for the dispersed phase representing solid particles and the fluid phase constituted of the separated gas (air) and liquid (water) phases. The mass conservation equations for the solid, liquid, and gas phases are given by

$$ \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_i} [\phi u_i^s] = 0, \qquad (1) $$

$$ \frac{\partial}{\partial t} [\gamma(1-\phi)] + \frac{\partial}{\partial x_i} [\gamma(1-\phi)u_i^l] = 0, \qquad (2) $$

$$ \frac{\partial}{\partial t} [(1-\gamma)(1-\phi)] + \frac{\partial}{\partial x_i} [(1-\gamma)(1-\phi)u_i^g] = 0, \qquad (3) $$

with $\phi$ the sediment volume concentration; $u_i^s$, $u_i^l$, and $u_i^g$ the solid-, liquid-, and gas-phase velocities, respectively; $x_i$ the position vector with $i = 1, 2, 3$ representing the three spatial components; and $\gamma$ the fluid indicator function defined as the ratio between the volume occupied by water and the total volume occupied by air and water in a cell ($\gamma = 1$ in water, $\gamma = 0$ in air, and $0 < \gamma < 1$ at the interface). $\gamma$ is also commonly referred to as the volume of the fluid (Hirt and Nichols, 1981).

Instead of solving the three mass conservation equations, mass conservation can be fully accounted for by solving the conservation equations for $\phi$ and $\gamma$. Combining Eqs. (1), (2), and (3) and assuming $u_i^l = u_i^g$ at the interface because of the no-slip boundary condition for velocities between air and
water, we obtain the conservation equation for the indicator function $\gamma$ following

$$ \frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} \left[ \gamma u_i^f \right] - \gamma \frac{\partial u_i^f}{\partial x_i} = 0, \quad (4) $$

with $u_i^f = \gamma u_i^l + (1 - \gamma) u_i^g$ representing the velocity of the fluid phase constituted of the air and water. More details about the derivation of Eq. (4) can be found in Appendix A.

Similarly, instead of solving momentum conservation equations for the three phases, the air and water mixture is considered a single fluid with varying densities and viscosities across the interface. The momentum conservation equations for the solid and fluid phases are given by

$$ \frac{\partial}{\partial t} [\rho^s \phi u_i^s] + \frac{\partial}{\partial x_i} [\rho^s \phi u_i^s u_j^s] = \frac{\partial}{\partial x_i} \Sigma_{ij}^s + \phi \rho^s g_i + \phi f_i + M_i, \quad (5) $$

$$ \frac{\partial}{\partial t} [\rho^f (1 - \phi) u_i^f] + \frac{\partial}{\partial x_i} [\rho^f (1 - \phi) u_i^f u_j^f] = \frac{\partial}{\partial x_i} \Sigma_{ij}^f + (1 - \phi) \rho^f g_i + (1 - \phi) f_i - M_i + \sigma \kappa \frac{\partial \gamma}{\partial x_i}, \quad (6) $$

with $\rho^s$ and $\rho^f = \gamma \rho^l + (1 - \gamma) \rho^g$ representing the solid and fluid densities, $\rho^l$ and $\rho^g$ the water- and air-phase densities, $\Sigma_{ij}^s$ and $\Sigma_{ij}^f$ the solid- and fluid-phase effective stress tensors defined later in Sect. 2.1, $g_i$ the gravitational acceleration, $f_i$ the momentum forcing term used to drive the flow, $M_i$ the momentum exchange term between the fluid and solid phases to be described in Sect. 2.2, $\sigma$ the surface tension coefficient, and $\kappa$ the local air-water interface curvature (more details about this term are available in Klostermann et al., 2013).

Compared with the two-phase model sedFoam, only one additional unknown field $\gamma$ is added to the model, and therefore only one additional equation for $\gamma$ (Eq. 4) needs to be solved. Also, the last term on the right-hand side of the fluid-phase momentum conservation equation (Eq. 6) needs to be included to model the surface tension between the air and water.

## 2.1 Effective stress tensors

The effective stress tensors of the fluid and solid phases are decomposed into normal and shear stresses following $\Sigma_{ij}^f = -P^f \delta_{ij} + T_{ij}^f$ and $\Sigma_{ij}^s = -P^s \delta_{ij} + T_{ij}^s$. Here, $P^f$ and $P^s$ denote the fluid and solid pressures, $\delta_{ij}$ is the Kronecker symbol, and $T_{ij}^f$ and $T_{ij}^s$ represent the fluid and solid shear stress tensors. These stress tensors are expressed in terms of the flow variables, which are defined as

$$ T_{ij}^f = \rho^f (1 - \phi) \times \left[ (\nu^{\text{mix}} + v_t^f) \left( \frac{\partial u_i^f}{\partial x_j} + \frac{\partial u_j^f}{\partial x_i} - \frac{2}{3} \frac{\partial u_k^f}{\partial x_k} \delta_{ij} \right) - \frac{2}{3} \lambda^f \delta_{ij} \right] \quad (7) $$

for the fluid phase and

$$ T_{ij}^s = \rho^s \phi \left[ (\nu^s + \nu_t^s) \left( \frac{\partial u_i^s}{\partial x_j} + \frac{\partial u_j^s}{\partial x_i} - \frac{2}{3} \frac{\partial u_k^s}{\partial x_k} \delta_{ij} \right) - \frac{2}{3} \lambda^s \delta_{ij} \right] \quad (8) $$

for the solid phase.

In these equations, $\nu^s$ denotes the solid-phase viscosity, $\nu_t^f$ and $\nu_i^s$ are the fluid- and solid-phase eddy (turbulent) viscosities, $\lambda^f$ and $\lambda^s$ are the fluid- and solid-phase bulk viscosities, and $\nu^{\text{mix}}$ is the mixture viscosity equal by default to the fluid-phase viscosity $\nu^f = \gamma \nu^l + (1 - \gamma) \nu^g$, with $\nu^l$ and $\nu^g$ the liquid- and gas-phase viscosities or functions of the sediment concentration when using the $\mu(I)$ rheology.

While air and water are considered Newtonian fluids with constant viscosities, the solid-phase viscosity $\nu^s$ together with the solid-phase bulk viscosity $\lambda^s$ and pressure $P^s$ are modeled using either the kinetic theory for granular flows or the $\mu(I)$ rheology. sedInterFoam is a direct extension of sedFoam documented in Chauchat et al. (2017). As a consequence, the solid-phase stress modeling in the two solvers is exactly the same. The reader is referred to Chauchat et al. (2017) for more details.

However, compared with the sedFoam version presented in Chauchat et al. (2017), the sedInterFoam user can choose between using several Reynolds-averaged Navier-Stokes (RANS) turbulence models, e.g., $k-\epsilon$ or $k-\omega$ models (Chauchat et al., 2017; Nagel et al., 2020), or using a dynamic Lagrangian large-eddy simulation (LES) model introduced by Mathieu et al. (2021) to model the subgrid eddy viscosities.

## 2.2 Momentum exchange between the phases

The momentum exchange term $M_i$ between the two phases is composed of buoyancy, drag, lift, added mass forces, and unresolved fluid-particle interaction forces $B_i$, $D_i$, $L_i$, $A_i$, and $I_i$, respectively. They are written in the following expressions:

$$ M_i = B_i + D_i + L_i + A_i + I_i \quad \text{with} \quad \left\{ \begin{array}{l} B_i = -\phi \frac{\partial P^f}{\partial x_i}, \\ D_i = \phi (1 - \phi) K (u_i^f - u_i^s), \\ L_i = \phi (1 - \phi) C_l \rho^m \| u_i^f - u_i^s \| \epsilon_{ijk} \frac{\partial u_k^m}{\partial x_j}, \\ A_i = \phi (1 - \phi) C_a \rho^f \times \left[ \frac{\partial u_i^f}{\partial t} + \frac{\partial u_i^f u_j^f}{\partial x_j} - \left( \frac{\partial u_i^s}{\partial t} + \frac{\partial u_j^s u_i^s}{\partial x_j} \right) \right], \\ I_i = -S_{US} (1 - \phi) K v_t^f \frac{\partial \phi}{\partial x_i} \end{array} \right. \quad (9) $$

where $C_l = 0.5$ and $C_a = 0.5$ are the lift and added mass coefficients, $u_i^m = \phi u_i^s + (1 - \phi) u_i^f$ and $\rho^m = \phi \rho^s + (1 - \phi) \rho^f$
are the mixture velocity and density, $\epsilon_{ijk}$ is the Levi-Civita symbol, $S_{US} = 1/Sc$ is the inverse of the Schmidt number, and $K = \rho^s/t_s(1-\phi)$ is the drag parameter with $t_s$ the particle response time modeled using a drag law. Several drag laws available in sedFoam are also present in sedInterFoam. As an example, the drag law proposed by Ding and Gidaspow (1990) combines the model of Ergun (1952) for high concentrations ($\phi > 0.2$) and the model of Wen and Yu (1966) for low concentrations ($\phi < 0.2$) following

$$t_s = \begin{cases} \rho^s \left( \frac{150\phi v^f \rho^f}{(1-\phi)d_p^2} + \frac{1.75\rho^f \|u_i^f - u_i^s\|}{d_p} \right)^{-1}, & \phi > 0.2, \\ \frac{4}{3} \frac{\rho^s}{\rho^f} \frac{d_p}{C_D \|u_i^f - u_i^s\|} (1-\phi)^{1.65}, & \phi < 0.2, \end{cases} \quad (10)$$

with $d_p$ the particle diameter and $C_D$ the drag coefficient given by

$$C_D = \frac{24}{Re_p} (1 + 0.15 Re_p^{0.687}). \quad (11)$$

The particle Reynolds number $Re_p$ is expressed as

$$Re_p = \frac{(1 - \phi)d_p \|u_i^f - u_i^s\|}{v^f}. \quad (12)$$

The model proposed for the unresolved fluid-particle interaction term $I_i$ should only be included for RANS simulations. Considering that no subgrid interaction model has been validated rigorously in sedInterFoam for LES, we recommend setting $S_{US}$ to zero. This term also plays a role in LES but is only important for particles with a diameter much smaller than the grid size. According to Ozel et al. (2013), it can be neglected when the grid size is on the order of the particle size.

### 3 Model implementation and solution algorithm

sedInterFoam is an extension of the Eulerian two-phase model for sediment transport applications, sedFoam (https://github.com/sedFoam/sedFoam, last access: 1 March 2025) (Chauchat et al., 2017; Mathieu et al., 2021), implemented using the open-source CFD library OpenFOAM (Jasak and Uroić, 2020). Mass and momentum equations are solved using a finite-volume method, and a pressure-implicit splitting of operators (PISO) algorithm is utilized for velocity-pressure coupling (Rusche, 2003).

As a comparison with sedFoam, the novelty lies in resolving the spatial and temporal evolution of free surfaces (air-water interface) by solving Eq. (4) using the volume-of-fluids (VOF) method and the Multidimensional Universal Limiter for Explicit Solution (MULES) algorithm (Rusche, 2003; Klostermann et al., 2013). The idea behind the VOF

method is to maintain $\gamma$ as a step function across the air-water interface and to ensure boundedness ($0 < \gamma < 1$). In other words, the interface has to be artificially compressed to balance numerical diffusion. To do so, we define the artificial relative velocity between air and water $u_i^f = u_i^g - u_i^l$ and rewrite Eq. (4) as

$$\frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} [\gamma u_i^f] - \gamma \frac{\partial u_i^f}{\partial x_i} + \frac{\partial}{\partial x_i} [\gamma (1 - \gamma) u_i^l] = 0. \quad (13)$$

The last term of Eq. (13) only acts in the region $0 < \gamma < 1$, and the relative velocity $u_i^f$ is explicitly estimated to compress the interface. More details about the numerical treatment of Eq. (13) are available in Klostermann et al. (2013).

After solving for the free-surface evolution, fluid-phase viscosity and density are updated and momentum conservation equations for the fluid and solid phases are solved as a two-phase system using the discretization procedure presented in Chauchat et al. (2017). For free-surface wave applications, a popular toolbox called waves2Foam (Jacobsen et al., 2012) is included for wave generation and absorption. The solution procedure is outlined as follows:

1. Solve for the indicator function $\gamma$ using the interface compression method (Eq. 13).

2. Update the curvature $\kappa$, viscosity $v^f$, and density $\rho^f$.

3. Solve for the sediment concentration $\phi$ (Eq. 1).

4. Call the waves2Foam library to update the wave generation (optional).

5. Update the momentum exchange term (Eq. 9).

6. Solve for solid-phase stress using the kinetic theory for granular flows or $\mu(I)$ rheology.

7. Solve for velocity-pressure coupling using a PISO loop.

8. Solve for turbulence closure (nothing for laminar, eddy viscosity for RANS, and subgrid closure for LES).

9. Go to the next time step.

The present model and implementation have been used with OpenFOAM-v2106 and OpenFOAM-v2112.

### 4 Model benchmarks and applications

In order to validate sedInterFoam and its implementation and to demonstrate its capability for challenging applications, four benchmarks available as tutorials along with the source code are presented in this section. A dam-break tutorial case included in the official OpenFOAM distribution is selected as the first benchmark to demonstrate the VOF capability of sedInterFoam (without sediment) compared with
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Sketch of the dam-break benchmark.</figcaption></figure>

interFoam. The second benchmark is the sedimentation tutorial distributed with sedFoam to verify the coupling between the solid and fluid phases. The third benchmark is a sheet flow configuration under monochromatic non-breaking waves measured in a large-wave flume by Dohmen-Janssen and Hanes (2002) in order to validate the coupling with waves2Foam for wave generation and the generation of sheet flows. Lastly, sedInterFoam is applied to simulate a solitary wave plunging on an erodible sloping beach, similar to the laboratory experiment of Sumer et al. (2011), in order to validate the model for predicting the beach profile evolution in the swash zone.

## 4.1 Dam break

A 2D dam-break tutorial distributed with OpenFOAM is reproduced with sedInterFoam to validate the free-surface algorithm without sediment.

A 0.15 m × 0.29 m water column initially at rest on the left-hand side of a 0.584 m × 0.584 m tank otherwise filled with air is released and impacts an obstacle at the center of the tank (Fig. 1). No turbulence model is used following the original tutorial (i.e., $v_t^f = 0$). The liquid-phase viscosity and density are specified as $v^l = 1 \times 10^{-6} \text{ m}^2\text{s}^{-1}$ and $\rho^l = 1000 \text{ kg m}^{-3}$, while the gas-phase viscosity and density are $v^g = 1.48 \times 10^{-5} \text{ m}^2\text{s}^{-1}$ and $\rho^g = 1 \text{ kg m}^{-3}$. The surface tension coefficient is $\sigma = 0.07$.

The same numerical parameters are used for both the interFoam tutorial case and the present test case for a fair comparison. No-slip boundary conditions are applied on the left, right, and bottom walls, and the top boundary is exposed to the atmosphere and permits both inflow and outflow. The mesh is composed of 9072 cells with a grid size on the order of 4 cm. An implicit first-order scheme (Eulerian) is used for time integration, and second-order upwind schemes are used to discretize the convection terms. The time step is adaptive and is calculated to ensure a Courant-Friedrichs-Lewy (CFL) number that is lower than 1. This case is run sequentially.

The comparisons between sedInterFoam and interFoam at 0, 0.3, and 0.6 s are presented in Fig. 2. The time evolution

of the air-water interface (defined by $\gamma = 0.5$) calculated by sedInterFoam is in very good agreement with interFoam. The shape and dynamics of the liquid phase impacting the obstacle are recovered, and generation of air pockets downstream of the obstacle is reproduced. The implementation of the free-surface evolution algorithm in sedInterFoam is therefore equivalent to that in interFoam. Discrepancies only appear at 0.6 s, when the flow becomes highly chaotic and initially small differences can lead to totally different results. Considering the numerical treatment of the three-phase flow equations, this behavior can be expected.

## 4.2 Sedimentation

The sedimentation benchmark taken from Chauchat et al. (2017) is reproduced with sedInterFoam but with an air-water interface specified in the top portion of the domain for comparison with experimental data obtained by Pham-Van-Bang et al. (2008) and validation of the solid and liquid coupling using the pressure-velocity algorithm, momentum exchange term, and mass conservation.

The configuration consists of a suspension of monodispersed polystyrene particles of diameter $d_p = 290 \text{ µm}$ and density $\rho^s = 1050 \text{ kg m}^{-3}$ in Rhodorsil silicone oil with a density of $\rho^l = 950 \text{ kg m}^{-3}$ and viscosity $\nu^l = 1.02 \times 10^{-5} \text{ m}^2\text{s}^{-1}$. The mixture is initially well-mixed in the bottom portion of the tank, with a particle concentration of $\phi = 0.5$ before sedimentation. The time evolution of the concentration is monitored in the experiment using a proton MRI device.

The numerical domain is 1D vertically decomposed into 240 cells over a depth of 12 cm. Compared with the sedFoam tutorial, an additional air layer with the same physical parameters as the dam-break benchmark presented in Sect. 4.1 is added at the top of the domain in order to be able to validate the mass conservation of the three phases throughout the simulation (Fig. 3). The time step is fixed at $\Delta t = 0.01 \text{ s}$. A first-order implicit time-integration scheme (Eulerian) and a second-order total variation diminish (TVD) scheme for convection terms are used. A wall boundary condition is applied at the bottom, and the pressure is fixed at the top to serve as a reference. The gradients of all the other quantities are set to zero at the top boundary. This case is run sequentially.

While settling, the sediment concentration profile shows two distinct interfaces. The upper interface corresponds to the transition between $\phi = 0$ and $\phi = 0.5$, while the lower interface corresponds to the transition between $\phi = 0.5$ and the maximum concentration in the settled bed. A comparison of the temporal locations of these two interfaces in the experiment and in the simulation is presented in Fig. 4a, while the concentration profiles at 232, 652, 1072, and 1492 s are shown in Fig. 4b. Very good agreements are observed with the measured data. For this benchmark, the performance of sedInterFoam is very similar to that of sedFoam reported in Chauchat et al. (2017). The spikes in the volume fraction for
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison between the $\gamma$ field from sedInterFoam represented by the color map and interFoam with the air-water interface represented by a dashed line (a–c). Differences between the $\gamma$ fields from sedInterFoam and interFoam (d–f) at 0 s (a, d), 0.3 s (b, e), and 0.6 s (c, f) for the dam-break configuration.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Sketch of the 1D sedimentation benchmark.</figcaption></figure>

the numerical results come from numerical oscillations at the interface because of the sharp concentration gradient. This behavior has also been observed with the sedFoam solver and results from the difficulty in handling the propagation of a shock by the mass conservation equations.

The temporal evolution of the total volume fraction difference compared with the initial values for the three phases is presented in Fig. 5. While the total volume fraction of the particles is conserved throughout the simulation, some variations are observed for the air and oil phases at the beginning. This is attributed to the time needed for the model to adapt from the initial conditions. After a few seconds, their total volume fractions evolve very slowly. Overall, the total volume fraction differences for the air and water phases remain lower than 0.5 %. To conclude, the model and its implementation

allow us to reproduce the pure sedimentation, agreeing with the measured data, and the mass conservation is preserved well.

## 4.3 Sheet flow under monochromatic waves

To validate the coupling with the waves2Foam library for surface wave generation and to show sedInterFoam capabilities on unstructured grids, an experimental configuration of the sheet flow driven by monochromatic non-breaking waves from Dohmen-Janssen and Hanes (2002) is reproduced numerically with sedInterFoam. Cnoidal waves with a wave period of $T = 6.5$ s and a wave height of $H = 1.55$ m are generated and propagate over a bed composed of sand particles with a median diameter $d_p = 240$ µm and a density $\rho^s = 2650$ kg m$^{-3}$ in a flume 3.5 m deep in the measurement section. Conductivity concentration meters buried in the sand measured the time series of the sediment concentration and allow for comparison with the model results.

The numerical configuration used for this benchmark is the same as the configuration investigated by Kim et al. (2018). In a Reynolds average field, the flow is homogeneous in the spanwise direction and a 2D configuration is established. The numerical flume is 151.56 m in length, and a sediment pit 4 m long and 0.1 m deep is located in the middle of the flume (Fig. 6). The mesh is generated with the OpenFOAM utilities blockMesh and snappyHexMesh. It is composed of around 2.85 million grid points with refined cells at the air-water and water-sediment interfaces. Details of the mesh close to
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Temporal evolution of the upper (blue) and lower (red) interfaces in the pure sedimentation configuration (panel a, with symbols for the experimental results and the solid line for the simulation results) and concentration profiles at *t* = 232, 652, 1072, and 1492 s (panel b, with the experimental results represented by the dashed blue line and the simulation results represented by the solid red line).</figcaption></figure>

the sand pit can be observed in Fig. 7. Wall boundary conditions are applied at the bottom, while the top boundary is exposed to the atmosphere, waves are generated using a 10th-order streamfunction on the left-hand side of the flume in the relaxation zone, and a sponge layer on the right-hand side of the flume allows us to absorb incoming waves and prevent reflection. The time step is adaptive to ensure a maximum CFL number lower than 0.4. A second-order backward time integrator is used, and a second-order TVD scheme is used for the convection terms. Turbulence is modeled using the two-phase $k-\epsilon$ model, and the kinetic theory for granular flows is used to model the shear-dependent solid-phase stresses. This case was run in parallel on 180 cores for 72 h.

To validate the wave generation implementation through coupling with waves2Foam and the model's capacity to predict sheet flow processes under surface waves, time series of the free stream velocity above the sediment bed and concentration profiles taken in the middle of the sediment pit at the wave crest, flow reversal, and wave trough are compared with measured data in Fig. 8. The good agreement in the time series of the free stream velocity highlights the ability of sedInterFoam to simulate nonlinear surface wave processes and validates a proper implementation of waves2Foam.

For the concentration profiles, numerical results are in good agreement with the measured data. Discrepancies can be observed at flow reversal, for which suspended sediments are slightly overpredicted. Turbulence closure parameters controlling the equilibrium between settling and upward turbulent diffusion of sediment may be responsible for the overpredicted sediment suspension. Careful tuning of the model parameters would allow more quantitative agreement between numerical and experimental results but is beyond the scope of this study. Results may also be sensitive to grid resolution and other turbulence models. The same conclusions apply for the comparison with the previous model, sedWaveFoam. The results are almost identical, except at flow reversal. Differences in the implementation of the equations and turbulence closure parameters can affect the vertical turbulent diffusion of sediment.

## 4.4 Plunging solitary wave

In order to demonstrate the novel capabilities of sedInterFoam for a more complex coastal application, an experimental configuration reported by Sumer et al. (2011) for a plunging solitary wave over an erodible bed is reproduced numerically. The solitary wave propagates in a flume with a 0.4 m water depth and breaks on top of a beach with a 1 : 14 slope composed of sand with a median diameter of $d_p = 180$ µm. Information about the solitary wave characteristics is given in Sumer et al. (2011). In the experiment, Sumer et al. (2011) recorded the beach profile evolution after the impact of four identical solitary waves. In the simulation, successive solitary waves are generated after intervals of 20 s. This has been shown to be enough time for the free surface to become still again. This is a critical test of sedInterFoam for simulating the challenging process of beach erosion in the swash zone.

Assuming spanwise homogeneity of the Reynolds average flow, the numerical configuration is 2D. The left part of the flume is 12 m long with a constant water depth of 0.4 m. The right portion of the flume consists of a 1 : 14 sloping beach loaded with sediment (see Fig. 9). The sediment bed is 0.15 m deep at the toe and 0.05 m deep at the top of the beach. The mesh is composed of 1.2 million grid points. The time step is adaptive to ensure a maximum CFL number lower than 0.3. The same numerical schemes as for the sheet flow in the monochromatic-wave configuration presented in Sect. 4.3 were used. The $k-\epsilon$ turbulence model is used, and
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Time evolution of the total volume fraction difference for the solid, liquid, and gas phases in the 1D sedimentation configuration.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Sketch of the sheet flow configuration under monochromatic waves, similar to the large-wave flume experiment of Dohmen-Janssen and Hanes (2002).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Snapshot of the sheet flow with a monochromatic-wave numerical configuration showing air in white, water in blue, and sediment in red (a), together with the details of the mesh (b). For visibility, the mesh is downsampled and the vertical scale is stretched by a factor of 7. The dashed lines in the top panel correspond to the relaxation zones for wave generation on the left and the sponge layer on the right.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Time series of the free stream velocity above the sediment pit (a) and the sediment concentration profiles at the wave crest (b), flow reversal (c), and wave trough (d).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Sketch of the plunging solitary wave configuration, similar to the laboratory experiment reported by Sumer et al. (2011).</figcaption></figure>

solid-phase stress is modeled using the $\mu(I)$ rheology. This case was run in parallel on six cores for 30 h (four successive solitary waves).

Snapshots of the simulation showing the impact of the first solitary wave are presented in Fig. 10. Please note that the maximum value of the color bar is larger than the maximum concentration in order to avoid saturation. At $t = 8.5$ s (top-left panel of Fig. 10), the incoming solitary wave starts to break and erodes sediment at $x = 5.5$ m (the origin of the $x$ axis corresponds to the toe of the sloping beach). The morphological changes resulting from the passage of breaking waves at around $x = 5.6$ to $5.7$ m are already visible at $t = 11$ s (top-right panel of Fig. 10), which corresponds to the end of the uprush phase. During backwash at $t = 13$ s (bottom-left panel of Fig. 10), a hydraulic jump reported by Sumer et al. (2011) is also observed between 4.5 and 5.2 m, under which sediment accumulates. Eventually, at $t = 20$ s (bottom-right panel of Fig. 10), the free surface returns to rest, and the accretive pattern between 4.5 and 5 m and the erosive pattern further landward become more evident.

Surface elevation numerically monitored at the toe of the beach ($x = 0$ m) and in the swash zone ($x = 4.63$ and $5.35$ m) is compared to the rigid-bed experimental data in Fig. 11. The shape of the incoming solitary wave is captured well by

sedInterFoam at the different locations across the domain. At $x = 4.63$ m (second panel of Fig. 11), a significant decrease in the water surface between 8 and 12 s is underestimated by the model. This time interval corresponds to the rundown phase when the newly formed hydraulic jump is migrating rapidly seaward, while strong shear suspends sediment close to the free surface. As a result, the differences observed in the free-surface elevation may be a consequence of interactions between sediment and water that are not present in the rigid-bed experiments.

To compare the predictive capabilities of sedInterFoam, the evolution of the bed profile after four solitary waves is compared to the experiments in Fig. 12. Following the laboratory experiment, each successive solitary wave is sent after the previous solitary wave impact to the flow field, and the bathymetry diminishes. The distinctive beach profile evolution features that can be observed from the measured data are accretions between 4 and 5 m corresponding to the location where a hydraulic jump is observed during backwash and erosion between 5 and 6 m, corresponding to the intermittently wet and dry areas of the beach (i.e., the upper swash). The position and amplitude of the morphological changes in the erosion and accretion zone between 3 and 7 m are predicted fairly well by sedInterFoam, with an error of less than 10% in the region between 2 and 7 m with respect to the distance between the highest deposition and deepest scouring points. More experimental points would be required to draw conclusions about the accuracy of the model in the swash zone between 6 and 7 m. However, the model does not reproduce the accretion zone located at around 1 m that corresponds to the formation of a sandbar with an error of around 20%. The processes involved in the generation of this morphological feature most certainly include interaction between
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Snapshots of the sediment concentration $\phi$ at $t = 8.5, 11, 13,$ and $20$ s from the simulation of the first solitary wave impacting the sloping beach. The local coordinate $x = 0$ m is defined at the toe of the sloping beach.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Free-surface elevation at $x = 0$ m (toe of the beach), $x = 4.63$ m, and $x = 5.35$ m from simulations with sedInterFoam (black solid lines) compared with rigid-bed experimental data from Sumer et al. (2011) (red dotted lines).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Morphological evolution of the sediment bed after four consecutive solitary waves from sedInterFoam simulation compared with the experimental results from Sumer et al. (2011).</figcaption></figure>
turbulence and sand particles in suspension. Such processes would require better parameterization of RANS turbulence models or use of turbulence-resolving simulations.

Overall, the presented simulation is a proof of concept where the morphodynamics in the swash zone can be simulated numerically using sedInterFoam. This new open-source model will allow us to further investigate fine-scale hydro-morphodynamic processes in order to answer the many open questions in this field.

## 5 Conclusions

In this paper, sedInterFoam, an extension of the two-phase flow solver sedFoam (Chauchat et al., 2017) implemented using the open-source library OpenFOAM (Jasak and Uroic, 2020), was presented. A third phase representing the air is included to model sediment transport applications driven by surface waves. The air-water interface is solved using the VOF method, similar to interFoam (Rusche, 2003; Klostermann et al., 2013) but adapted for the present miscible liquid-solid phase. Coupling with the waves2Foam library (Jacobsen et al., 2012) is implemented to generate and absorb free-surface waves. sedInterFoam includes all the previous features present in sedFoam, such as turbulence models (RANS or LES) and solid-phase stress models (kinetic theory for granular flows and $\mu(I)$ rheology).

The model has been applied successfully to a 2D dam-break configuration to benchmark the numerical results with those produced by the existing OpenFOAM solver interFoam in order to verify the implementation of the free-surface evolution algorithm. A comparison between sedInterFoam and sedFoam has been performed by simulating a 1D sedimentation configuration in order to confirm that the implementation of a third phase did not affect the mass conservation of each phase. A sheet flow configuration under monochromatic waves from Kim et al. (2018) was reproduced successfully, indicating that the implementation of waves2Foam for wave generation and absorption is appropriate. Eventually, an experimental configuration of a solitary wave plunging on a sandy beach with a slope of 1:14 (Sumer et al., 2011) was reproduced numerically to highlight the capabilities of sedInterFoam in simulating complex wave breaking, swash dynamics, and beach profile evolution. Quantitative agreement between the measured and simulated results is obtained, particularly regarding the erosion and deposition processes in the swash zone.

The development of the open-source three-phase flow model sedInterFoam provides a modeling tool for investigating coastal sediment transport applications dominated by surface waves. While relatively computationally expensive, the physics-based model sedInterFoam can be a useful tool for gaining insight into the complex physical processes associated with breaking waves, sediment transport, and morphodynamics and providing simulation data to improve em-

pirical parameterizations in regional-scale morphodynamic models.

## Appendix A: Derivation of the indicator function transport equation

In this section, the steps to derive the $\gamma$ transport equation are presented. Summing the mass conservation equations for the solid, liquid, and gas phases (Eqs. 1, 2, and 3) and assuming $u_i^l = u_i^g = u_i^f$ because of the non-slip boundary condition at the air-water interface allows us to write the mixture conservation equation as

$$ \frac{\partial}{\partial x_i} [\phi u_i^s + (1-\phi) u_i^f] = 0. \qquad (A1) $$

Expanding Eq. (A1) and rearranging the terms gives

$$ \frac{\partial}{\partial x_i} [\phi u_i^f] = \frac{\partial}{\partial x_i} [\phi u_i^s] + \frac{\partial u_i^f}{\partial x_i}. \qquad (A2) $$

Starting from the liquid-phase mass conservation Eq. (2), expanding it, and using the chain rule for the derivation of the products gives

$$ \frac{\partial}{\partial t} [\gamma(1-\phi)] + \frac{\partial}{\partial x_i} [\gamma(1-\phi)u_i^f] = 0, \qquad (A3) $$

$$ \Rightarrow \frac{\partial \gamma}{\partial t} - \frac{\partial \phi \gamma}{\partial t} + \frac{\partial}{\partial x_i} [\gamma u_i^f] - \frac{\partial}{\partial x_i} [\gamma \phi u_i^f] = 0, \qquad (A4) $$

$$ \Rightarrow \frac{\partial \gamma}{\partial t} - \phi \frac{\partial \gamma}{\partial t} - \gamma \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_i} [\gamma u_i^f] \qquad (A5) $$

$$ -\phi \frac{\partial}{\partial x_i} [\gamma u_i^f] - \gamma u_i^f \frac{\partial \phi}{\partial x_i} = 0. $$

Rearranging the terms of Eq. (A5) allows us to write

$$ (1-\phi) \left[ \frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} [\gamma u_i^f] \right] - \gamma \left[ \frac{\partial \phi}{\partial t} + u_i^f \frac{\partial \phi}{\partial x_i} \right] = 0, \qquad (A6) $$

$$ \Rightarrow (1-\phi) \left[ \frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} [\gamma u_i^f] \right] \qquad (A7) $$

$$ -\gamma \left[ \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_i} [\phi u_i^f] - \phi \frac{\partial u_i^f}{\partial x_i} \right] = 0. $$

Replacing the second-to-last term of Eq. (A7) with Eq. (A2) and remembering the solid-phase mass conservation Eq. (1) reads as

$$ (1-\phi) \left[ \frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} [\gamma u_i^f] \right] - \gamma \underbrace{\left[ \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_i} [\phi u_i^s] + \frac{\partial u_i^f}{\partial x_i} - \phi \frac{\partial u_i^f}{\partial x_i} \right]}_{=0} = 0. \qquad (A8) $$

Rearranging the terms allows us to obtain the final form of the $\gamma$ transport equation following
$$ \Rightarrow (1 - \phi) \left[ \frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} \left[ \gamma u_i^f \right] \right] - \gamma \frac{\partial u_i^f}{\partial x_i} + \gamma \phi \frac{\partial u_i^f}{\partial x_i} = 0, \quad (A9) $$

$$ \Rightarrow (1 - \phi) \left[ \frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} \left[ \gamma u_i^f \right] - \gamma \frac{\partial u_i^f}{\partial x_i} \right] = 0, \quad (A10) $$

$$ \Rightarrow \frac{\partial \gamma}{\partial t} + \frac{\partial}{\partial x_i} \left[ \gamma u_i^f \right] - \gamma \frac{\partial u_i^f}{\partial x_i} = 0. \quad (A11) $$

**Code availability.** sedInterFoam is distributed under the GNU General Public License v2.0 (GNU GPL v2.0) and is available from Zenodo at the following DOI: https://doi.org/10.5281/zenodo.10577879 (Mathieu et al., 2024).

**Data availability.** The benchmarks from this article are available in the code repository on Zenodo for reproducibility (https://doi.org/10.5281/zenodo.10577879, Mathieu et al., 2024).

**Author contributions.** Model development: AM, YK, CB, JC. Validation case development: AM, YK, CB, JC, TJH. Project management: YK, JC, TJH. Manuscript preparation: AM with contributions from all the co-authors.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

**Acknowledgements.** The numerical simulations presented in this study benefited from computing resources provided by the DARWIN and Caviness clusters at the University of Delaware and the EXPANSE cluster at the San Diego Supercomputer Center via ACCESS.

**Financial support.** This research has been supported by the National Research Foundation of Korea (grant no. NRF-2021R1F1A1062223), the Korea Institute of Ocean Science and Technology (grant no. PEA-0131), the Directorate for Engineering (grant no. CMMI-2050854), the Directorate for Geosciences (grant no. OCE-2242113), the Office of Naval Research (grant no. N00014-22-1-2412), and the Strategic Environmental Research and Development Program (grant no. MR-1478).

**Review statement.** This paper was edited by Andy Wickert and reviewed by two anonymous referees.

## References

* Amoudry, L. O.: Extension of $k-\omega$ turbulence closure to two-phase sediment transport modelling: Application to oscillatory sheet flows, Adv. Water Resour., 72, 110–121, 2014.

* Baykal, C., Sumer, B., Fuhrman, D., Jacobsen, N., and Fredsøe, J.: Numerical simulation of scour and backfilling processes around a circular pile in waves, Coast. Eng., 122, 87–107, 2017.

* Boyer, F., Guazzelli, E., and Pouliquen, O.: Unifying Suspension and Granular Rheology, Phys. Rev. Lett., 107, 188301, https://doi.org/10.1103/PhysRevLett.107.188301, 2011.

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

* Chauchat, J.: A comprehensive two-phase flow model for unidirectional sheet-flows, J. Hydraul. Res., 56, 15–28, 2018.

* Chauchat, J. and Guillou, S.: On turbulence closures for two-phase sediment-laden flow models, J. Geophys. Res.-Oceans, 113, C11017, https://doi.org/10.1029/2007JC004708, 2008.

* Chauchat, J., Cheng, Z., Nagel, T., Bonamy, C., and Hsu, T.-J.: SedFoam-2.0: a 3-D two-phase flow numerical model for sediment transport, Geosci. Model Dev., 10, 4367–4392, https://doi.org/10.5194/gmd-10-4367-2017, 2017.

* Cheng, Z., Hsu, T.-J., and Calantoni, J.: SedFoam: A multi-dimensional Eulerian two-phase model for sediment transport and its application to momentary bed failure, Coast. Eng., 119, 32–50, 2017.

* Delisle, M.-P. C., Kim, Y., Mieras, R. S., and Gallien, T. W.: Numerical investigation of sheet flow driven by a near-breaking transient wave using SedFoam, Eur. J. Mech. B-Fluids, 96, 51–64, 2022.

* Ding, J. and Gidaspow, D.: A bubbling fluidization model using kinetic theory of granular flow, AIChE J., 36, 523–538, 1990.

* Dohmen-Janssen, C. M. and Hanes, D. M.: Sheet flow dynamics under monochromatic nonbreaking waves, J. Geophys. Res.-Oceans, 107, 13-1–13-21, 2002.

* Dong, P. and Zhang, K.: Intense near-bed sediment motions in waves and currents, Coast. Eng., 45, 75–87, 2002.

* Ergun, S.: Fluid Flow through Packed Columns, Chem. Eng. Prog., 48, 89–94, 1952.

* Hirt, C. and Nichols, B.: Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys., 39, 201–225, 1981.

* Hsu, T.-J., Jenkins, J. T., and Liu, P. L.-F.: On two-phase sediment transport: sheet flow of massive particles, P. Roy. Soc. Lond. A Mat., 460, 2223–2250, 2004.

* Jacobsen, N. G. and Fredsøe, J.: Formation and development of a breaker bar under regular waves. Part 2: Sediment transport and morphology, Coast. Eng., 88, 55–68, 2014.

* Jacobsen, N. G., Fuhrman, D. R., and Fredsøe, J.: A wave generation toolbox for the open-source CFD library: OpenFoam®, Int. J. Numer. Meth. Fl., 70, 1073–1088, 2012.

* Jasak, H. and Uroic, T.: Practical Computational Fluid Dynamics with the Finite Volume Method, Springer, 103–161, https://doi.org/10.1007/978-3-030-37518-8_4, 2020.
Kim, Y., Cheng, Z., Hsu, T.-J., and Chauchat, J.: A Numerical Study of Sheet Flow Under Monochromatic Nonbreaking Waves Using a Free Surface Resolving Eulerian Two-Phase Flow Model, J. Geophys. Res.-Oceans, 123, 4693–4719, 2018.

Kim, Y., Mieras, R. S., Cheng, Z., Anderson, D., Hsu, T.-J., Puleo, J. A., and Cox, D.: A numerical study of sheet flow driven by velocity and acceleration skewed near-breaking waves on a sandbar using SedWaveFoam, Coast. Eng., 152, 103526, https://doi.org/10.1016/j.coastaleng.2019.103526, 2019.

Kim, Y., Mieras, R. S., Anderson, D., and Gallien, T.: A Numerical Study of Sheet Flow Driven by Skewed-Asymmetric Shoaling Waves Using SedWaveFoam, J. Mar. Sci. Eng., 9, 936, https://doi.org/10.3390/jmse9090936, 2021.

Klostermann, J., Schaaake, K., and Schwarze, R.: Numerical simulation of a single rising bubble by VOF with surface compression, Int. J. Numer. Meth. Fl., 71, 960–982, 2013.

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

Lee, C.-H., Xu, C., and Huang, Z.: A three-phase flow simulation of local scour caused by a submerged wall jet with a water-air interface, Adv. Water Resour., 129, 373–384, 2019.

Mathieu, A., Chauchat, J., Bonamy, C., and Nagel, T.: Two-Phase Flow Simulation of Tunnel and Lee-Wake Erosion of Scour below a Submarine Pipeline, Water, 11, 1727, https://doi.org/10.3390/w11081727, 2019.

Mathieu, A., Chauchat, J., Bonamy, C., Balarac, G., and Hsu, T.-J.: A finite-size correction model for two-fluid large-eddy simulation of particle-laden boundary layer flow, J. Fluid Mech., 913, A26, https://doi.org/10.1017/jfm.2021.4, 2021.

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

Mathieu, A., Kim, Y., Hsu, T.-J., Bonamy, C., and Chauchat, J.: sedInterFoam, Zenodo [data set and code], https://doi.org/10.5281/zenodo.10577879, 2024.

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

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

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

Ozel, A., Fede, P., and Simonin, O.: Development of filtered Euler-Euler two-phase model for circulating fluidised bed: High resolution simulation, formulation and a priori analyses, Int. J. Multiphas. Flow, 55, 43–63, 2013.

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

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

Rusche, H.: Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions, PhD thesis, Imperial College London, https://www.researchgate.net/publication/271830940_Computational_Fluid_Dynamics_of_Dispersed Two-Phase_Flows_at_High_Phase_Fractions (last access: 1 March 2025), 2003.

Salimi-Tarazouj, A., Hsu, T.-J., Traykovski, P., and Chauchat, J.: Eulerian Two-Phase Model Reveals the Importance of Wave Period in Ripple Evolution and Equilibrium Geometry, J. Geophys. Res.-Earth, 126, e2021JF006132, 2021a.

Salimi-Tarazouj, A., Hsu, T.-J., Traykovski, P., Cheng, Z., and Chauchat, J.: A Numerical Study of Onshore Ripple Migration Using a Eulerian Two-phase Model, J. Geophys. Res.-Oceans, 126, e2020JC016773, https://doi.org/10.1029/2020JC016773, 2021b.

Sherwood, C. R., van Dongeren, A., Doyle, J., Hegermiller, C. A., Hsu, T.-J., Kalra, T. S., Olabarrieta, M., Penko, A. M., Rafati, Y., Roelvink, D., van der Lugt, M., Veeramony, J., and Warner, J. C.: Modeling the Morphodynamics of Coastal Responses to Extreme Events: What Shape Are We In?, Annu. Rev. Mar. Sci., 14, 457–492, 2022.

Sumer, B., Truelsen, C., Sichmann, T., and Fredsøe, J.: Onset of scour below pipelines and self-burial, Coast. Eng., 42, 313–335, 2001.

Sumer, B. M., Sen, M. B., Karagali, I., Ceren, B., Fredsøe, J., Søttila, M., Zilioli, L., and Fuhrman, D. R.: Flow and sediment transport induced by a plunging solitary wave, J. Geophys. Res.-Oceans, 116, C01008, https://doi.org/10.1029/2010JC006435, 2011.

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

van der A, D. A., Ribberink, J. S., van der Werf, J. J., O’Donoghue, T., Buijsrogge, R. H., and Kranenburg, W. M.: Practical sand transport formula for non-breaking waves and currents, Coast. Eng., 76, 26–42, 2013.

Wen, C. Y. and Yu, Y. H.: A generalized method for predicting the minimum fluidization velocity, AIChE J., 12, 610–612, 1966.

Yu, X., Hsu, T.-J., Jenkins, J. T., and Liu, P. L.-F.: Predictions of vertical sediment flux in oscillatory flows using a two-phase, sheet-flow model, Adv. Water Resour., 48, 2–17, 2012.