Impact Statement
Flow slides and breaching failures represent a major risk for buildings, roads and other infrastructures. Additionally, in the dredging industry, the breaching process is used to extract sand and it may become unstable which may result in a loss of land. Although these phenomena are broadly observed in nature, the lack of understanding of the underlying physical processes is partly due to the absence of accurate experimental data and partly due to the complexity of the models. Laboratory experiments are challenging due to the impossibility of seeing through a dense suspension of sand particles. Numerical models, thus, may be a promising alternative to gain insight into these complex failure mechanisms and help engineers to better assess the risk of flow slides and breaching. This contribution is one of the first attempts to develop a physically consistent two-phase flow model to predict the dynamics of a wide range of failure modes.
1. Introduction
Although granular flows are ubiquitous in nature and industrial applications, researchers still struggle to completely understand the underlying physics of such flows. Difficulties arise when the granular material interacts with a viscous fluid: the deformation of the soil skeleton induces changes on the fluid pressure field which, subsequently, affects the topology of the soil skeleton. Indeed, the complex inner nature of such flows is the main reason for the absence of simplified models that describe immersed granular flows.
In this work, we explore the collapse of granular columns immersed in a viscous fluid. Special attention is given to Reference Rondon, Pouliquen and AussillousRondon, Pouliquen, and Aussillous (2011), an experimental investigation of the role of the initial solid volume fraction in the dynamics of the granular collapse. Reference Rondon, Pouliquen and AussillousRondon et al. (2011) observed that in initially loose packings, the excess of pore pressure built up and enhanced a rapid mobility of the granular column, whereas in initially dense granular packings, negative pore pressure developed within the porous medium increasing the shear resistance and, therefore, delaying the granular collapse. This mechanism, due to the contracting/dilating behaviour of granular material, is commonly known as pore pressure feedback and was first reported by Reference Iverson, Reid and LahusenIverson, Reid, and Lahusen (1997) and Reference IversonIverson (2005) and largely studied experimentally (Reference Bougouin and LacazeBougouin & Lacaze, 2018; Reference Iverson, Reid, Iverson, Lahusen, Logan, Mann and BrienIverson et al., 2000; Reference Pailha, Nicolas and PouliquenPailha, Nicolas, & Pouliquen, 2008; Reference Rondon, Pouliquen and AussillousRondon et al., 2011) and numerically (Reference Bouchut, Fernández-nieto, Mangeney and Narbona-reinaBouchut, Fernández-nieto, Mangeney, & Narbona-reina, 2016; Reference LeeLee, 2021; Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà, Chauchat, Chareyre, Bonamy, & Hsu, 2021; Reference Wang, Wang, Peng and MengWang, Wang, Peng, & Meng, 2017a) using multiple configurations where the dense/loose granular material is sheared. Overall, on the collapse of immersed granular columns, literature (Reference LeeLee, 2021; Reference Rondon, Pouliquen and AussillousRondon et al., 2011) distinguishes two main processes largely linked to the initial volume fraction: on the one hand, initially loose packings lead to shear failure forming a sliding wedge. On the other hand, dense granular packings are more prone to collapse through the breaching mechanism, which is the process of front particles progressively released as a turbidity current while the fluid penetrates the granular column to enhance the dilation of the medium. Breaching failure carries potential danger in submerged dense sandy soils. Indeed, Reference Beinssen, Neil and MastbergenBeinssen, Neil, and Mastbergen (2014) have reported multiple large-scale failures due to the breach face slowly receding from the original position. Even though this paper focuses on the breaching and shear plane failures, it is worth mentioning that other failure mechanisms, such as soil liquefaction, are also influenced by dilatancy effects (Reference PrevostPrevost, 1985; Reference YoudYoud, 1973). Liquefaction occurs in loose soils where a rapid particle rearrangement leads to a pore pressure build-up vanishing the effective stress. Structures on liquefiable soils may have terrible consequences under earthquakes or other shear-induced situations (Reference Koutsourelakis, Prévost and DeodatisKoutsourelakis, Prévost, & Deodatis, 2002).
Several numerical studies of immersed granular collapses have been reported in the literature. Reference Kumar, Delenne and SogaKumar, Delenne, and Soga (2017) explored the effect of initial volume fraction on the dynamics of two-dimensional (2-D) granular collapses by means of the discrete element method (DEM) coupled with the lattice Boltzmann method (LBM). Reference Izard, Lacaze, Bonometti and PedronoIzard, Lacaze, Bonometti, and Pedrono (2018) were able to reproduce three-dimensional granular collapses with an immersed boundary method coupled with DEM. Similarly, Reference Xu, Dong and DingXu, Dong, and Ding (2019) and Reference Yang, Jing, Kwok and SobralYang, Jing, Kwok, and Sobral (2019) relied on smoothed particle hydrodynamics–DEM and LBM–DEM approaches, respectively, to study the process of submerged granular collapse. Most previous work is computationally expensive and time-consuming. Thus, affordable simulations are typically restricted to a low number of particles, which limits the range of applications, especially if one is interested in large-scale applications such as coastal breaching or landslides. Alternatively, continuum approaches are more suited for large-scale problems but their accuracy highly depends on the closure models and an adequate coupling between the fluid and the solid phase. On the one hand, the mixture model proposed by Reference Savage, Babaei and DabrosSavage, Babaei, and Dabros (2014) was capable of predicting the initial dynamics of loose granular collapse but less satisfactory results were found for dense granular collapse. This approach neglected the excess of pore pressure despite the fact that the pore pressure feedback mechanism plays a key role in the collapse. On the other hand, Reference Bouchut, Fernández-Nieto, Koné, Mangeney and Narbona-ReinaBouchut, Fernández-Nieto, Koné, Mangeney, and Narbona-Reina (2017) proposed a depth-averaged approach to model the dilatancy effects and pore pressure feedback mechanism of different submerged granular collapses. Based on the Eulerian–Eulerian framework, the first attempts by Reference Lee, Huang and ChiewLee, Huang, and Chiew (2015) and Reference Lee and HuangLee and Huang (2018) to model the intricate dynamics of granular column collapses led to promising results. However, such models adopted simple elastic relationships to determine the solid pressure ignoring the shear-induced volume changes which led to an insufficient description of the pore pressure feedback. Reference Shi, Dong, Yu and ZhouShi, Dong, Yu, and Zhou (2021) and Reference LeeLee (2021) went one step further and proposed modifications of the elastic equation that captured the sheared-induced volume changes and reproduced the granular collapse dynamics with great success. Previous studies are mainly focused on the morphology of the deposit and the pore pressure feedback mechanism of small-scale granular collapses. Yet, experimental studies (Reference Eke, Viparelli and ParkerEke, Viparelli, & Parker, 2011; Reference Vandenberg, Vangelder and MastbergenVandenberg, Vangelder, & Mastbergen, 2002; Reference van Rhee and Bezuijenvan Rhee & Bezuijen, 1999) show that the breaching process forming turbidity currents is a crucial mechanism to reproduce practical applications such as dredging engineering or protection measures against coastal erosion (Reference Beinssen, Neil and MastbergenBeinssen et al., 2014; Reference Mastbergen, Beinssen and NédélecMastbergen, Beinssen, & Nédélec, 2019; Reference ShipwayShipway, 2015; Reference Vandenberg, Vangelder and MastbergenVandenberg et al., 2002). Accordingly, it seems reasonable to upscale the problem and examine whether the collapse is driven by the same dynamics as observed in the experiments of Reference Rondon, Pouliquen and AussillousRondon et al. (2011) or, conversely, whether other physics apply. Although studies of breaching are scarce (Reference Alhaddad, Labeur and UijttewaalAlhaddad, Labeur, & Uijttewaal, 2020; Reference Alhaddad, Weij, Vanrhee and KeetelsAlhaddad, Weij, Vanrhee, & Keetels, 2023; Reference van Rhee and Bezuijenvan Rhee & Bezuijen, 1999; Reference WeijWeij, 2020; Reference You, Flemings and MohrigYou, Flemings, & Mohrig, 2012, Reference You, Flemings and Mohrig2014), most laboratory experiments reported that the granular columns manifest a dual-mode slope failure. Instead of observing a progressive erosion induced by the turbidity currents of the breaching process, the evolution of the granular column is governed by a combination of breaching and occasional sliding wedges. Reference You, Flemings and MohrigYou et al. (2014) suggest that slide failures take place when the negative pore pressure developed in the porous medium is not enough to keep the shear resistance against slide failure. Reference You, Flemings and MohrigYou et al. (2014) remarked that slides are associated with a drop in excess pore pressure, and the magnitude of the jump is related to the size of the slide. Additionally, experiments (Reference Lee and ChenLee & Chen, 2022; Reference Lee and KuanLee & Kuan, 2021) have reported that dense granular packings exhibit shear failure for coarse particles and breaching for fine particle. It is worth noting that experiments investigating the breaching process are limited in size because large-scale experiments are not affordable. Thus, numerical simulations arise as a potential alternative to study the physics of the breaching process in large-scale applications. Reference Breusers, Nicollet and ShenBreusers, Nicollet, and Shen (1977) introduced the concept of active wall velocity defined as the horizontal velocity at which steep slopes move due to the breaching process. Since then, several breaching erosion models have been proposed based on different closures for the wall velocity. A quasi-static one-dimensional depth-averaged approach developed by Reference Mastbergen and VandenbergMastbergen and Vandenberg (2003) was used to investigate the breaching process showing that turbidity currents can be strong enough to periodically flush large deposits of sediments from canyons. Reference Eke, Viparelli and ParkerEke et al. (2011) considered a similar model to study another flushing event in submarine canyons. More recently, Reference van Rheevan Rhee (2015) proposed a 2-D drift-flux model based on the Reynolds-averaged Navier–Stokes equations to study the stability of the breaching process. Previous studies, nonetheless, are subject to some limitations. Firstly, current models do not account for slide failures and, secondly, they neglect the evolution of soil properties. Therefore, in order to further investigate the breaching process, numerical models should be able to capture not only the turbidity currents, but also the transition from static to yielding soil and the effects of the pore pressure feedback.
In this work, we first validate the elastoplastic model presented in Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al. (2021) using the experiments of Reference Rondon, Pouliquen and AussillousRondon et al. (2011) as reference. Then, we report on the application of the numerical model to study the effects of the breaching process of a granular column collapse 
$\sim$35 times larger (Reference WeijWeij, 2020) than that from the experiments of Rondon et al.
2. Mathematical formulation
The governing equations for the Eulerian–Eulerian two-phase formulation are shown below along with the closure forms for drag force, turbulence model and stresses for the fluid and particle phases.
2.1 Two-phase flow governing equations
The mass continuity equations for the solid and fluid phase are written as follows:
$$\begin{gather} \frac{\partial \phi }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{u^s} \phi ) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial (1-\phi) }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u^f} (1-\phi) ) = 0. \end{gather}$$ Here, 
$\phi$, 
$\boldsymbol {u^s}$ and 
$\boldsymbol {u^f}$ are the solid volume faction, the particle phase velocity and the fluid phase velocity, respectively.
The momentum conservation equations for the solid phase and fluid phase are written as
\begin{align} \rho^s \phi \left[ \frac{\partial \boldsymbol{u^s} }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\boldsymbol{u^s}\otimes \boldsymbol{u^s} \right)\right] & = \phi(\rho^s-\rho^f) \boldsymbol{g} + \dfrac{(1-\phi) \rho^f \nu^f}{K} (\boldsymbol{u^f} - \boldsymbol{u^s}) \nonumber\\ &\quad -\boldsymbol{\nabla} p^s + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau^s} - \phi \boldsymbol{\nabla} p^f, \end{align}
\begin{align} \rho^f (1-\phi) \left[\frac{\partial \boldsymbol{u^f} }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u^f}\otimes \boldsymbol{u^f})\right] &= \frac{(1-\phi) \rho^f \nu^f}{K} (\boldsymbol{u^s} - \boldsymbol{u^f}) + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau^f} - (1-\phi) \boldsymbol{\nabla} p^f, \end{align}
where 
$\rho ^s$ is the solid density, 
$\rho ^f$ is the fluid density, 
$\nu ^f$ stands for the fluid kinematic viscosity, 
$\otimes$ is the outer product of two vectors, 
$p^f$ is the excess of pore pressure defined as the difference between the pore pressure and the hydrostatic pressure, 
$p^s$ is the solid pressure, 
$\boldsymbol {\tau ^s}$ is the granular shear stress, 
$\boldsymbol {\tau ^f}$ is the fluid shear stress and 
$K$ is the permeability of the porous medium.
This work uses the approach of Reference EngelundEngelund (1953) to compute the permeability based on the pressure drop in steady porous flow:
\begin{equation} K = \frac{{\rm d}^2 (1-\phi)}{\alpha_E \phi^2 } + \frac{{\rm d} \nu^f (1-\phi)^2}{\beta_E \left\| \boldsymbol{u^f} - \boldsymbol{u^s} \right\| }, \end{equation}
where 
$d$ is the mean particle diameter. According to Reference Burcharth and ChristensenBurcharth and Christensen (1991) and Reference Burcharth and AndersenBurcharth and Andersen (1995), 
$\alpha _E$ ranges from 780 (for uniform and spherical particles) to 1500 or more (for irregular and angular grains), while 
$\beta _E$ ranges from 1.8 (for uniform and spherical particles) to 3.6 or more (for irregular and angular grains). The choice of the Reference EngelundEngelund (1953) model over other approaches employed for high packing fractions (such as Reference ErgunErgun, 1952) lies in the advantage of adjusting the factors 
$\alpha _E$ and 
$\beta _E$ to accurately represent the permeability of soils made up of irregular-shaped particles.
2.1.1 Solid and fluid stress
The fluid phase shear stress is expressed as the sum of the Reynolds stresses due to turbulent fluctuations (
$\boldsymbol {R^{t}}$) and the viscous stresses (
$\boldsymbol {r^f}$):
\begin{equation} \boldsymbol{\tau^f} = \boldsymbol{R^{t}} + \boldsymbol{r^{f}}, \end{equation}
where the Reynolds stress tensor 
$\boldsymbol {R^{t}}$ is modelled as
\begin{equation} \boldsymbol{R^{t}} = 2 \rho^f (1-\phi) \left [ \nu^{t} \boldsymbol{S^f} -\tfrac{1}{3}k \boldsymbol{I} \right ] \end{equation}and the viscous stress is written as
\begin{equation} \boldsymbol{r^{f}} = 2 \rho^f (1-\phi) \nu^{mix} \boldsymbol{S^f}, \end{equation}where
\begin{equation} \boldsymbol{S^f} = \tfrac{1}{2} (\boldsymbol{\nabla} \boldsymbol{u^f} + (\boldsymbol{\nabla} \boldsymbol{u^f})^{\rm T}) - \tfrac{1}{3} {\rm tr}(\boldsymbol{\nabla} \boldsymbol{u^f}) \end{equation}
is the deviatoric and symmetric part of the velocity gradient for the fluid phase, 
$k$ is the turbulent kinetic energy, 
$\nu ^{t}$ stands for the eddy viscosity calculated with a turbulence closure model (see § 2.1.2) and 
$\nu ^{mix}$ is the mixture viscosity, which, according to Reference Boyer, Guazzelli and PouliquenBoyer, Guazzelli, and Pouliquen (2011), can be computed with the following phenomenological expression:
\begin{equation} \nu^{mix} = \nu^{f} \left [1 + 2.5 \phi \left ( 1- \frac{\phi}{\phi_{max}}\right )^{{-}1} \right ], \end{equation}
where 
$\phi _{max}$ is the maximum solid volume fraction set to 
$\phi _{max}=0.625$ for spheres.
Following Reference Cheng, Hsu and CalantoniCheng, Hsu, and Calantoni (2017) and Reference Chauchat, Cheng, Nagel, Bonamy and HsuChauchat, Cheng, Nagel, Bonamy, and Hsu (2017), the solid phase pressure 
$p^s$ is defined as the sum of a viscous shear-rate-dependent pressure 
$p_s^s$ and the contribution of enduring contacts 
$p_c^s$:
\begin{equation} p^s = p^s_s + p^s_c, \end{equation}
where 
$p_c^s$ is proportional to the difference between the solid volume fraction 
$\phi$ and the reference solid fraction 
$\phi _{pl}$ where dilatancy effects are embedded:
\begin{equation} p^s_{c} = \begin{cases} 0 & \phi < \phi_{pl}\\ E \dfrac{(\phi-\phi_{pl})^3}{(\phi_{rcp}-\phi)^5} & \phi \geq \phi_{pl}, \end{cases} \end{equation}
where 
$\phi _{rcp}$ is the random close packing volume fraction. We adopt the value for sphere packings (
$\phi _{rcp}=0.625$). It is important to note that (2.12) is based on the work of Reference Johnson and JacksonJohnson and Jackson (1987), which assumes a constant value called the random loose packing fraction (
$\phi _{rlp}$ instead of 
$\phi _{pl}$). However, to accurately account for the effects of dilatancy, we need to consider initial and transient packing fractions that are different from the random loose packing fraction. The initial volume fraction (
$\phi _{o}$) is calculated as the average volume fraction throughout the height of the bed, which is given by
\begin{equation} \phi_{o} = \frac{1}{h_o}\int_{o}^{h_o} \phi(y,t = 0)\,{{\rm d}y}. \end{equation}
Here, 
$h_o$ represents the lowest position above which 
$\phi \le \phi _{top}=0.53$. By adjusting the initial plastic volume fraction (
$\phi _{pl,t=0}$), which remains constant during the process of gravitational deposition for preparing the sample, we can achieve different initial volume fractions. Higher values (in this study, 
$\phi _{pl,t=0} = 0.609$ to match Reference Rondon, Pouliquen and AussillousRondon et al. (2011)) result in initially dense granular beds (
$\phi _o \approx 0.61$), while lower values (in this study, 
$\phi _{pl,t=0} = 0.54$ to reproduce Reference Rondon, Pouliquen and AussillousRondon et al. (2011)) yield initially loose granular packings (
$\phi _o \approx 0.55$).
Once the system reaches an equilibrium state and the numerical sedimentation is complete, the granular collapse begins and 
$\phi _{pl}$ evolves to account for the plastic effects. Following Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al. (2021), the plastic effects that arise from local rearrangements during shearing deformations are captured as an increase/reduction of 
$\phi _{pl}$, which further changes the solid pressure (see (2.12)). More specifically, the expression that governs the evolution of 
$\phi _{pl}$ is
\begin{equation} \frac{\partial\phi_{pl}}{\partial t} + \boldsymbol{u^s} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi_{pl} = {-} \phi_{pl} \delta \| \boldsymbol{S^s}\|, \end{equation}
where 
$\boldsymbol {S^s}$ is the deviatoric strain rate of the solid phase computed as in (2.9) but replacing the fluid velocity with the solid velocity and 
$\delta$ is the dilatancy coefficient defined as
\begin{equation} \delta = K_1 (\phi-\phi_{\infty}), \end{equation}
where 
$K_1$ is a calibration parameter and 
$\phi _{\infty }$ stands for the equilibrium volume fraction. As reported by Reference Roux and RadjaiRoux and Radjai (1998) and Reference Pailha and PouliquenPailha and Pouliquen (2009), the linear variation of the dilatancy angle with the volume fraction as written in (2.15) is derived from the linerization of the dilation rate:
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u^s} = \frac{1}{\phi} \frac{{\rm d} \phi}{{\rm d} t} = \delta \| \boldsymbol{S^s}\| . \end{equation}
Parameter 
$\phi _{\infty }$ is modelled as a function of the particle pressure and the shear rate through the viscous number:
\begin{equation} \phi_{\infty} = \frac{\phi_{c}}{1 + I_v^{1/2}}, \end{equation}
where the viscous number (
$I_v$) is defined as
\begin{equation} I_v = \frac{\rho^f \nu^f \| \boldsymbol{S^s}\| }{p^s} \end{equation}
and 
$\phi _{c}$ is the critical volume fraction in quasti-static shear (
$I_v \rightarrow 0$).
In this work 
$\delta$ is limited to a range of 
$-0.4 \leq \delta \leq 0.4$ which falls into a range of physical values proposed by previous works (Reference Alshibli and CilAlshibli & Cil, 2018; Reference Iverson and GeorgeIverson & George, 2014; Reference Pouliquen and RenautPouliquen & Renaut, 1996). The influence of the dilatancy prefactor 
$K_1$ is be studied in § 3.1.3.
Equation (2.17) suggests that different values of the equilibrium solid volume fraction are expected in transient conditions, such as the onset of granular collapse, the fully developed flow and the final arrest. In the context of a granular column collapse, steady flows are unlikely: viscous forces are expected to slow down the granular flow before it is fully developed; therefore, dilatancy effects will arise provided that 
$\phi \neq \phi _\infty$ and 
$\| \boldsymbol {S^s}\| > 0$. It is noteworthy that despite the fact that 
$\phi _{pl}$ is just a numerical parameter, the consequences of changes in 
$\phi _{pl}$ are completely physical. Indeed, this model not only extends the critical state soil mechanics to a rate-dependent critical state (
$\phi _{\infty }$) but also leads to an increase or a decrease of the granular pressure depending on the initial packing and, consequently, to pressure-driven expansion or compaction of the solid phase under shear conditions. Furthermore, the value of 
$\phi _{pl}$ remains unbounded because during granular flow it cannot decrease or increase indefinitely. This is because for loose materials (low 
$\phi _{pl}$ value), the dilatancy coefficient is negative, and according to (2.14), 
$\phi _{pl}$ must increase. Similarly, for dense packings, it is also not possible for 
$\phi _{pl}$ (initially large) to continuously increase because positive dilatancy coefficients lead to a reduction in 
$\phi _{pl}$.
The expression for the shear-rate-dependent pressure induced by collisional interactions was derived by Reference Chauchat, Cheng, Nagel, Bonamy and HsuChauchat et al. (2017) inverting (2.17) to give the rate-dependent normal stress 
$p^s_{\infty }$:
\begin{equation} p^s_{\infty} = \rho^f \nu^f \| \boldsymbol{S^s}\| \left ( 1- \frac{\phi_{c}}{\phi} \right )^{{-}2}. \end{equation} However, as suggested by Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al. (2021), 
$p^s_{\infty }$ is consistently defined to be the stationary shear-induced pressure whereas the actual pressure is supposed to converge asymptotically to that value with accumulated strain. Therefore, the following equation governs the progressive mobilization of 
$p_s^s$:
\begin{equation} \frac{\partial p^s_s}{\partial t} + \boldsymbol{u^s} \boldsymbol{\cdot} \boldsymbol{\nabla} p^s_s = {-} K_2 ( p^s_s-p^s_{\infty}) \| \boldsymbol{S^s}\| . \end{equation} In short, dilatancy is a complex physical phenomenon and our approach consists of decomposing its effect into the enduring contact pressure 
$p_c^s$ and the shear-rate-dependent pressure 
$p_s^s$. The former pressure (
$p_c^s$) is closely related to the microstructure. For instance, in initially dense packings, the particles are interlocking and are not able to move freely. Therefore, the grains need to be rearranged to allow shearing deformations. During grain reorganization, the contacts become stronger which comes with an increase of the enduring contact pressure 
$p_c^s$. The shear-rate-dependent pressure 
$p_s^s$, on the contrary, is derived from the 
$\mu (I_v)$ rheology. Reference Boyer, Guazzelli and PouliquenBoyer et al. (2011) showed the granular medium dilates when increasing the shear rate (
$I_v \uparrow$) which is accompanied by an increase of the solid pressure that scales with the fluid viscosity and the shear rate.
Reference Jop, Forterre and PouliquenJop, Forterre, and Pouliquen (2006) proposed that the ratio between shear stress and pressure can be scaled by the inertial number, 
$I={{\rm d} \| \boldsymbol {S^s}\| }/{\sqrt {p^s / \rho ^s}}$, defined as the ratio between the macro and micro time scales of granular flow. Reference Cassar, Nicolas and PouliquenCassar, Nicolas, and Pouliquen (2005) followed a similar approach using the viscous number 
$I_v$ to model granular flows immersed in viscous fluids. Even though the present work utilizes the viscous number 
$I_v$ because the analysed granular collapses are found to be in the viscous regime (the Stokes number 
$St={\rho ^s d^2 \| \boldsymbol {S^s}\|}/{\rho ^f \nu ^f}$ ranges from 
$St=0.005$ to 
$St=0.3$ depending on the granular collapse), different scenarios or real large-scale events, where the rheology belongs to an inertial regime rather than a viscous regime, may be studied with the present numerical model by simply adapting the constitutive laws for the inertial number 
$I$ as reported by Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al. (2021).
In this work, the solid shear stress 
$\boldsymbol {\tau ^s}$ is proportional to the solid pressure following a frictional law depending on the viscous number 
$I_v$:
\begin{equation} \boldsymbol{\tau^s} = \mu(I_v) p^s \frac{\boldsymbol{S^s}}{ \| \boldsymbol{S^s}\|}, \end{equation}
where 
$\mu (I_v)$ is the friction coefficient for a certain shear state described in Reference Boyer, Guazzelli and PouliquenBoyer et al. (2011) as
\begin{equation} \mu(I_v) = \mu_s + \frac{ \Delta \mu }{ I_o/I_v + 1}, \end{equation}
where the empirical material constants correspond to the static friction coefficient 
$\mu _s$, the dynamic friction coefficient 
$\Delta \mu$ and the reference viscous number 
$I_o$.
In order to have an expression for 
$\tau ^s$ resembling the definition for the fluid shear stress, the shear stress due to frictional contacts can be rewritten as
\begin{equation} \boldsymbol{\tau^s} = \rho^s \nu^s \boldsymbol{S^s}, \end{equation}
where 
$\nu ^s$ is the frictional shear viscosity:
\begin{equation} \nu^s = \frac{\mu(I_v) p^s}{ \rho^s \left ( \| \boldsymbol{S^s}\|^2 + \lambda^2_r \right ) ^{1/2}}, \end{equation}
where 
$\lambda _r$ is a regularization parameter from Reference Chauchat and MédaleChauchat and Médale (2014) taken to be 
$\lambda _r=10^{-6}$ s
$^{-1}$. Moreover, 
$\nu ^s$ is limited to be smaller than 
$\nu ^s_{max}=10^{5}\,{\rm m}^2\,{\rm s}^{-1}$ to avoid numerical issues.
2.1.2 Turbulence model
To model the turbulent eddy viscosity 
$\nu ^t$, the 
$k$–
$\epsilon$ model (Reference Cheng, Hsu and CalantoniCheng et al., 2017; Reference Hsu and LiuHsu & Liu, 2004) is used in this study. Parameter 
$\nu ^t$ is computed as
\begin{equation} \nu^t = C_{\mu} \frac{k^2}{\epsilon}, \end{equation}
where 
$C_{\mu }=0.09$ is an empirical coefficient, 
$k$ is the turbulent kinetic energy and 
$\epsilon$ is the dissipation rate of turbulent kinetic energy. The turbulent kinetic energy 
$k$ is determined with the following transport equation:
\begin{equation} \frac{\partial k}{\partial t} + u^f_j \frac{\partial k}{\partial x_j} = \frac{R^{t}_{ij}}{\rho^f} \frac{\partial u^f_i}{ \partial x_j} + \frac{\partial}{\partial x_j} \left[\left (\nu^f + \frac{\nu^t}{\sigma_k} \right ) \frac{\partial k}{\partial x_j} \right ] -\epsilon - \frac{2 K (1-t_{mf}) \phi k }{\rho^f} , \end{equation}
where 
$R^t_{ij}$ is the Reynolds stress tensor, 
$\sigma _k$ is an empirical coefficient and 
$t_{mf}$ is a parameter that characterizes the degree of correlation between particles and fluid velocity fluctuations modelled as 
$t_{mf} = \exp ({-B \cdot St})$, in which 
$B$ is an empirical coefficient chosen as 
$B = 0.25$ and 
$St$ is the Stokes number defined as 
$St={t_p}/{t_l}$, where 
$t_p=\rho ^s K \phi / ((1-\phi ) \rho ^f \nu ^f)$ is the particle response and 
$t_l=k/(6 \epsilon )$ is the characteristic time scale of energetic eddies.
The balance equation for the dissipation rate of turbulent kinetic energy is written as
\begin{equation} \frac{\partial \epsilon}{\partial t} + u^f_j \frac{\partial \epsilon}{\partial x_j} = C_{1 \epsilon} \frac{\epsilon}{k} \frac{R^{t}_{ij}}{\rho^f} \frac{\partial u^f_i}{ \partial x_j} + \frac{\partial}{\partial x_j} \left [ \left (\nu^f + \frac{\nu^t}{\sigma_{\epsilon}} \right ) \frac{\partial \epsilon}{\partial x_j} \right ] -C_{2 \epsilon}\frac{\epsilon^2}{k} -C_{3 \epsilon} \frac{\epsilon}{k} \frac{2 \nu^f (1-t_{mf}) k }{K}. \end{equation}For this work, the values of the turbulent empirical coefficients are found in table 1.
Empirical coefficients for the 
$k$–
$\epsilon$ turbulence model taken from Reference Chauchat, Cheng, Nagel, Bonamy and HsuChauchat et al. (2017).
2.2 Numerical implementation
Simulations are conducted with the open-source software SedFoam, a two-phase flow solver used for sediment transport applications (Reference Chassagne, Bonamy and ChauchatChassagne, Bonamy, & Chauchat, 2023; Reference Chauchat, Cheng, Nagel, Bonamy and HsuChauchat et al., 2017; Reference Mathieu, Chauchat, Bonamy and NagelMathieu, Chauchat, Bonamy, & Nagel, 2019; Reference Mathieu, Cheng, Chauchat, Bonamy and HsuMathieu, Cheng, Chauchat, Bonamy, & Hsu, 2022; Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al., 2021; Reference Tsai, Mathieu, Montellà, Hsu and ChauchatTsai, Mathieu, Montellà, Hsu, & Chauchat, 2022) based on the open-source finite volume library OpenFOAM (Reference Jasak and UroićJasak & Uroić, 2020) (v2212 release from ESI). The solver is available for download from GitHub (https://github.com/SedFoam/SedFoam). Several interpolation/discretization techniques can be used to evaluate the face fluxes. Table 2 shows schemes used for temporal and spatial discretization. It is worth noting that the Gauss upwind scheme is only adopted for the divergence discretization of 
$\phi _{pl}$ and 
$p^s_s$ fields. In order to resolve the velocity–pressure coupling, we rely on the pressure-implicit split-operator (PISO) algorithm.
Numerical schemes for the interpolation of the convective fluxes.
3. Results
In this section, the two-phase flow model is used to reproduce initially loose and dense granular column collapses. The first part of this section consists of reproducing the granular collapses of Reference Rondon, Pouliquen and AussillousRondon et al. (2011). Although, a broad number of works (Reference Izard, Lacaze, Bonometti and PedronoIzard et al., 2018; Reference Jing, Yang, Kwok and SobralJing, Yang, Kwok, & Sobral, 2019; Reference Meng, Liao, Yu, Li and AnMeng, Liao, Yu, Li, & An, 2021; Reference Polanía, Cabrera, Renouf and AzémaPolanía, Cabrera, Renouf, & Azéma, 2022; Reference Riffard and RisRiffard & Ris, 2022; Reference Sun, Zhang, Wang and LiuSun, Zhang, Wang, & Liu, 2020; Reference Wang, Wang, Peng and MengWang et al., 2017a) have predicted the main features of granular column collapses immersed in a viscous fluid, only a few (Reference Wang, Wang, Peng and MengWang, Wang, Peng, & Meng, 2017b; Reference Yang, Jing, Kwok and SobralYang, Jing, Kwok, & Sobral, 2020) have successfully captured the dynamics of Rondon's experiments for both initially loose and dense granular columns, and even fewer have done it with a continuum approach (Reference Baumgarten and KamrinBaumgarten & Kamrin, 2019; Reference Bouchut, Fernández-Nieto, Koné, Mangeney and Narbona-ReinaBouchut et al., 2017; Reference LeeLee, 2021; Reference Phan, Bui and NguyenPhan, Bui, & Nguyen, 2022; Reference RauterRauter, 2021; Reference Shi, Dong, Yu and ZhouShi et al., 2021; Reference Si, Shi and YuSi, Shi, & Yu, 2018). The dilatancy model presented herein was able to capture the pore pressure feedback in one-dimensional and 2-D granular avalanches (Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al., 2021) with reasonably good agreement. The granular column collapse is, thus, used to extend the model to a more complex and realistic configuration. Additionally, a sensitivity analysis is summarized in this section to underline the parameters that govern the dynamics of the granular collapse complementing previous work on this topic (Reference LeeLee, 2021; Reference RauterRauter, 2021). Once the model is validated, the second part of this section is devoted to gain insight into the breaching process. In order to achieve this goal, the laboratory experiments of Reference WeijWeij (2020) and Reference Alhaddad, Weij, Vanrhee and KeetelsAlhaddad et al. (2023) are reproduced numerically using a 2-D approach.
3.1 Collapse of a granular column: Rondon's experiments
The experiments performed by Reference Rondon, Pouliquen and AussillousRondon et al. (2011) investigated the collapse of a granular column in a viscous liquid. Initially dense columns resulted in negative pore pressures that slowed down the collapse, while in loose granular packings, the collapsing process was triggered instantaneously with positive pore pressures that enhanced a rapid flow. Although several volume fractions and aspect ratios were analysed, only two representative cases (an initially loose column with 
$\phi _o \approx 0.55$ and an initially dense packing with 
$\phi _o \approx 0.61$) are considered in this work. The experimental set-up consisted of a tank with a length of 0.7 m and a width and height of 0.15 m. A removable gate was placed vertically and glass beads were poured behind the gate. The rest of the physical and geometrical parameters have been taken from Reference Rondon, Pouliquen and AussillousRondon et al. (2011) (see table 3). During the experiments, the excess of pore pressure was measured at the bottom at 
$2$ cm from the left-hand side of the tank. The numerical set-up is presented in figure 1. The numerical domain is decomposed into square cells of 
$0.41\,{\rm mm} \times 0.41$ mm.
Physical and geometrical variables used in the numerical simulations.
Numerical set-up to predict the granular collapse of Reference Rondon, Pouliquen and AussillousRondon et al. (2011).
According to Reference Rondon, Pouliquen and AussillousRondon et al. (2011), the internal friction angle is around 
$20^{\circ }$ and the critical volume fraction is 
$\phi _c=0.58$. In this work, however, we adopt the rheological parameters as in Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al. (2021) with a critical volume fraction of 
$\phi _c=0.57$. The rest of the rheological coefficients and calibration parameters are summarized in table 4. The permeability in (2.4) is modelled according to Engelund's model presented before with the coefficients 
$\alpha _E=780$ and 
$\beta _E=1.8$ corresponding to smooth spherical particles. This set of parameters led to the best fit to the experimental results. The influence of these parameters is further discussed in § 3.1.3. Moreover, one may argue that 
$K_2$ values should be taken the same for both loose and dense scenarios. Ideally, 
$K_2$ should be 
$1$ so the relaxation time is dominated solely by the shear rate (see (2.20)). However, numerical instabilities in the dense case forced us to set an additional relaxation. The choice of 
$K_2=0.01$, nonetheless, has a minor influence on the results because the inherent slow flow of dense granular collapse is driven mainly by changes in the contact pressure (see § 3.1.3). Finally, numerical simulations of Rondon's experiments belong to the viscous regime for the fluid phase; consequently, the turbulent viscosity is set to zero (
$\nu ^t=0$) for simplicity.
Rheological and numerical parameters used to reproduce Reference Rondon, Pouliquen and AussillousRondon et al. (2011).
3.1.1 Morphology
Figure 2 shows the evolution of the deposit shapes during the granular column collapses. As reported by Reference Rondon, Pouliquen and AussillousRondon et al. (2011), the dynamics of the granular column collapse is very different depending on the initial volume fraction: initial dense granular packings are mobilized slowly and show short run-out distances (see figure 2a). On the contrary, initially loose granular packings are characterized by a rapid flow and elongated fronts (see figure 2b).
Evolution of the morphology and solid volume fraction during the collapse of an initially (a) dense and (b) loose column. A grey line is included to illustrate the evolution of the isoline with the initial volume fraction (
$\phi _o = 0.55$ for the initially loose column and 
$\phi _o = 0.61$ for the initially dense column).
In addition to the dynamics of the granular flow, dilatancy effects are also illustrated in figure 2 as changes in the solid volume fraction. The sheared region in figure 2(a) is expanding progressively from 
$\phi _0 \approx 0.61$ to 
$\phi \approx 0.57$. Conversely, figure 2(b) shows a contraction of the sheared region from 
$\phi _0 \approx 0.55$ to 
$\phi \approx 0.565$. It is worth noting the abrupt change of concentration along a straight line displayed in figure 2(b) for the loose granular collapse. This line splits the moving and the static regions, and it is commonly referred to as the failure line. The collapse of the upper right part of the column is triggered right after the gate is removed. The sliding region is contracting more rapidly close to the bottom and the failure line because the shear rate is higher near the non-moving regions (the variable 
$\phi _{pl}$ embedding the dilatancy effects is proportional to the shear rate – see (2.14)). The failure line and the fluid velocity field are also visualized in figure 3. As mentioned before, only the sliding wedge is moving significantly, hence contracting. This leads to an expansion along the failure line to ensure the conservation of mass. Figure 3(c) illustrates the dilation rate as defined by Reference Iverson and GeorgeIverson and George (2014), i.e. the divergence of the solid phase velocity (
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u^s}$) so we can distinguish the contractancy and dilatancy regions. After the gate removal, 
$\phi _{pl}$ increases along the failure line. According to (2.12), the contact pressure is reduced; therefore, a reduction of the shear strength is expected enhancing a rapid flow slide. As the collapse carries on, the contractancy behaviour of the sheared region entails an increase of the volume fraction, which according to (2.12), is accompanied by an increase of the contact pressure. Figures 2(b) and 3 also show a mild expansion within the non-moving region. One may wonder why the so-called non-moving region is deforming if we defined it as static region. It is worth mentioning that the expansion of the static region is mainly caused by the reduction of the column height and its subsequent decompression. As a matter of fact, the expansion rate of this region is significantly low (
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u^s}$ values in figure 3c are 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u^s} \approx 0.004\,{\rm s}^{-1}$ whereas 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u^s} \approx -0.020\,{\rm s}^{-1}$ inside the contracting band). The gentle expansion of the static region is accompanied by an inward flux (see figure 3c) that occupies the growing pore space.
(a) Solid volume fraction and zoom-in view along the failure surface with (b) detailed volume fraction and (c) divergence of the solid phase velocity and fluid flow field for the initially loose granular column. It must be noted that the arrows displayed in (c) represent the magnitude of the fluid velocity by their colour and not their size.
Figure 2 exhibits qualitatively good agreement with the experimental data (Reference Rondon, Pouliquen and AussillousRondon et al., 2011) showing comparable time scales and remarkably similar run-out distances. Although there is an undeniable resemblance between the numerical and the experimental final deposits, figure 2 suggests that the numerical solution for the dense collapse leads to a milder final slope compared with the experimental data and the loose case is slightly slower during the first moments after gate removal.
3.1.2 Excess of pore pressure
Figure 4(a) shows that negative pore pressure develops, stabilizing the dense granular material. The initial vertical front has a much steeper slope than the angle of repose. Because of this unstable configuration, shearing deformation is triggered. Under such circumstances, the granular material dilates, pore bodies are enlarged and the fluid phase flows inwards to the porous medium to accommodate the expansion. Consequently, negative pore pressure is generated increasing the effective strength, and, overall, stabilizing the granular material. The negative pore pressure is, therefore, responsible for the characteristic creeping flow observed in figures 2(a) and 4(a).
Evolution of the morphology and excess of pore pressure (
$p^f$) during the collapse of an initially (a) dense and (b) loose column. A grey line is included to illustrate the zero pressure isoline. (c) Evolution of basal pore pressure (
$p^f$) measured at 2 cm (dark continuous line) and 3 cm (light dashed line). Shaded areas correspond to the region between the two probe results.
The loose scenario portrayed in figure 4(b) shows a dual positive/negative pore pressure map revealing the complex dynamics of the collapse. The contracting behaviour of the sheared region comes along with an expulsion of the pore fluid. Therefore, positive pore pressure develops within the moving area reducing the shear resistance and enhancing a rapid granular flow. As mentioned before, dilation also occurs in the loose collapse along the failure line (see figure 3c). Subsequently, the fluid is sucked into the failure line and partially into the non-moving region due to its decompression leading to negative pore pressures. This phenomenon was already reported by Reference LeeLee (2021): the solid volume fraction along the failure line decreases inducing a reduction of the solid pressure (see (2.12)), thus a lower shear strength. As a result of the combination of positive pore pressure within the sliding zone and a significant decrease of contact forces along the failure line, the shearing region is partially fluidized and the lower effective stress is accompanied by a rapid slide failure. The pore pressure map depicted in figure 4(b) is consistent with the pore pressure feedback mechanism and the sliding failure reported in the experiments. However, pore pressure values at the bottom of the loose column (see red point located at 
$2$ cm from the end of the tank in figure 4b) differ from the experimental measurements. Figure 4(c) shows that the experimental data and the numerical solution have a similar trend for the loose and dense granular collapse: a pore pressure jump is registered after gate removal that gradually dissipates. The negative pore pressure jump simulated for the dense granular collapse is properly predicted; however, the dissipation dynamics observed in figure 4(c) is slightly different. Even though the failure mode and the pore pressure feedback mechanism are generally well reproduced, the positive pore pressure values are underpredicted for the loose case. It is worth noting that the model fails to replicate the exact position of the failure line, and therefore the pore pressure probe (at 2 cm from the end of the tank) falls into the non-moving area instead of the sliding region where pore pressures present two distinct behaviours: negative pore pressure within the static zone and positive pore pressure in the sliding region. For this reason, figure 4(c) includes the pore pressure numerical measurements between 2 and 3 cm.
3.1.3 Sensitivity study
This section summarizes the role of different parameters on the dynamics of the granular collapse. Figures for the sensitivity analysis are provided as supplementary material available at https://doi.org/10.1017/flo.2023.23. The results of the sensitivity study were used to find the optimal set of parameters to reproduce Rondon's experiments numerically.
The 
$K_1$ sensitivity. The dilatancy prefactor 
$K_1$ introduced in (2.15) is responsible for controlling the plastic effects that arise from particle rearrangements during shear deformations. From the results shown in the supplementary material for the dense case, we observe that large 
$K_1$ values result in slightly slower creeping flow while the final deposit shape is unaltered. Dilatancy is relevant right after gate removal; then, dilatancy effects fade out and the granular collapse starts flowing after 
$t \approx 6$ s as it would flow at the equilibrium state (
$\phi =\phi _{\infty }$), eventually reaching the same deposit shape regardless of the 
$K_1$ value. Numerical results show that the dense granular collapse follows very similar dynamics for both 
$K_1=4$ (value proposed in Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al., 2021) and 
$K_1=40$ (reference value in this article). The 
$K_1=4$ scenario provides a slightly faster collapse but a better pore pressure dissipation curve matching the experimental points with striking accuracy.
Concerning the loose case, the supplementary material illustrates that increasing 
$K_1$ has a strong effect on the dynamics of the spreading deposit. According to (2.14) and (2.15), the increase of 
$K_1$ is responsible for 
$\phi _{pl}$ increasing more rapidly. Microscopically, it means contact forces are reduced more abruptly leading to a lower solid pressure and shear strength. Consequently, large 
$K_1$ values enhance a rapid flow with longer run-out distances. Adopting 
$K_1=4$ for the loose case is not enough to trigger compacting effects and the corresponding pore pressure feedback; meanwhile, the largest value (
$K_1=100$) reveals a high positive pore pressure jump. However, the morphology of the deposit remains barely affected. The lack of difference between 
$K_1=40$ and 
$K_1=100$ in terms of deposit shape for both loose and dense collapses is a consequence of the dilatancy coefficient (
$\delta$) limits imposed to keep 
$\delta$ bounded to the physical values reported in Reference Pouliquen and RenautPouliquen and Renaut (1996), Reference Iverson and GeorgeIverson and George (2014) and Reference Alshibli and CilAlshibli and Cil (2018). Values of 
$\delta$ are significantly important after gate removal. At this point, both 
$K_1=40$ and 
$K_1=100$ scenarios reach the limit 
$| \delta | =0.4$ in some regions of the granular column; thus, no relevant difference is observed between 
$K_1=40$ and 
$K_1=100$. As we approach the equilibrium state, 
$\delta$ values decay suppressing the dilatancy effect.
The results presented in this subsection indicate that, overall, good agreement is found between the experimental data and the numerical simulations. However, it is worth mentioning that the choice of 
$K_1$ is strongly influenced by the reference critical volume fraction, which is discussed in the next subsection; therefore, one may find other optimal 
$K_1$ values after increasing/decreasing the critical volume fraction 
$\phi _c$. This dense–loose asymmetric trend in terms of 
$K_1$ sensitivity suggests that the dilatancy model still has room for improvement. Indeed, the dilatancy model is governed by the evolution of the dilation angle through changes in the plastic volume fraction using the first term of the Taylor expansion of the dilatancy expression given by (2.16). Additionally, the elasto-plastic expression given by (2.12) is a crude simplification to model the stress field neglecting the anisotropy and non-local effects of real soils. Thus, the current formulation is a simplified approach to model dilatancy effects based on the amount of plastic volumetric strain. Further research could explore the use of nonlinear dilatancy laws and perform DEM simulations to improve or introduce new phenomenological expressions to predict the effects of dilatancy with even greater accuracy.
The 
$\phi _c$ sensitivity. The influence of the critical volume fraction (
$\phi _c$) is closely linked to the dilatancy effects. According to (2.14) and (2.17), dilatancy effects are proportional to 
$\phi -\phi _c/(1+I_v^{1/2})$. Therefore, assuming lower 
$\phi _c$ values result in lower contractancy effects for initially loose cases and stronger dilatancy effects for initially dense packings (see figure in the supplementary material). Likewise, larger 
$\phi _c$ values are associated with weaker dilatancy in dense cases and enhanced contractancy for loose granular columns. Changes in 
$\phi _c$ are markedly more important in the loose packings. Contractancy effects partially fluidize the granular column, which leads to rapid collapse with a final deposit of a very gentle slope. Correspondingly, the pore pressure dynamics is significantly higher for the loose case. In particular, the 
$\phi _c=0.58$ case reproduces the magnitude of the positive pore pressure jump reported in the experiments. On the contrary, minor differences in terms of pore pressure are observed for the scenarios with initially dense columns.
The 
$K_2$ sensitivity. In this section, we examine the impact of the parameter 
$K_2$ in (2.20). Parameter 
$K_2$ affects how quickly the shear-induced pressure reaches its equilibrium state. Figures in the supplementary material demonstrate that there is barely no difference in terms of excess pore pressure and deposit spreading response. It is important to note that the numerical simulations in this study are conducted under a dense viscous granular flow regime. Therefore, it is not surprising that the contact pressure has a much greater impact compared with the shear-induced pressure. In the dense scenario 
$K_2=1$ poses numerical stability issues that require a smaller time step. Thus, 
$K_2=0.01$ is preferred for the dense granular collapse to circumvent numerical instabilities.
Elastic modulus. This section explores the influence of the elastic modulus (
$E$) of (2.12) on the granular collapse dynamics. Before getting into the discussion, it is pertinent to note that 
$E$ values used in the present numerical model are considerably far from those of real materials such as glass beads (
$E \approx 70$ GPa) or sand (
$E \in [5\unicode{x2013}80]$ MPa); however, these values would induce numerical issues with the nonlinear approach of (2.12). Instead, 
$E$ values used in the present sensitivity analysis remain in a lower range (
$E \in [0.1\unicode{x2013}100]$ Pa) where numerical instabilities are not detected. Limited differences are observed for the dense granular collapse in terms of pore pressure and deposit shape. However, results are more sensitive in the loose scenario. As detailed in the supplementary material, differences arise as a consequence of changes on the shape of the vertical concentration profile. The nature of the contact pressure expression (2.12) plays a key role in the distribution of the initial volume fraction along the vertical by increasing the concentration vertical gradient with soft elastic modulus. A low elastic modulus (i.e. 
$E=0.1$ Pa) leads to a vertical concentration curve that ranges from 
$\phi _{top}=0.525$ to 
$\phi _{bottom}=0.565$ (
$\Delta \phi = \phi _{bottom} - \phi _{top} \approx 0.04$) a whereas stiffer modulus (i.e. 
$E=100$ Pa) has a narrower range (
$\Delta \phi = \phi _{bottom} - \phi _{top} \approx 0.01$). Choosing a low elastic modulus leads to a dual behaviour within the granular column; the material close to the bottom shows the classic features of a dense soil while the material located close to the surface presents a very loose-like behaviour. It is, therefore, recommended to use a high elastic modulus, provided that the numerical model remains stable, in order to have a realistic soil behaviour.
Permeability coefficients. Figures in the supplementary material show lower permeabilities lead to a slow mobilization. In turn, the pressure dissipation takes longer as expected, specially for the initially dense column. Figures in the supplementary material evidence that different Engelund coefficients have a minor impact on the results for the loose granular column regarding the morphology and the pore pressure curve.
Frictional coefficients. The friction coefficient (see (2.22)) has a certain effect on the shape of the deposit. Large friction coefficients delay the collapse and the deposit ends up with a steeper slope. The lower mobilization entails a weaker pore pressure feedback: the soil is more difficult to shear; thus, pore volume changes take longer. Conversely, low fiction angles promote a rapid failure with abrupt pore volume changes; therefore, higher pore pressure jumps are observed.
Discussion of the sensitivity study. In this work we have optimized the set of parameters based on the shape of the deposit. The sensitivity analysis of the present section reveals that a different set (
$\phi _c \uparrow$, 
$K_1 \uparrow$ and 
$\mu _s \downarrow$) would definitely provide a better pore pressure prediction for the loose case to the detriment of the prediction of the deposit morphology. Nevertheless, none of the combinations (except for 
$\phi _c=0.58$) presented in the sensitivity analysis reaches the positive pore pressure developed within the loose granular column. Several factors could explain the underprediction of the positive pore pressure curve. (1) Results are very sensitive to the initial volume fraction. Indeed, Reference LeeLee (2021) reported different 
$p^f$–time curves and deposit shapes for a narrow range of initial concentrations: 
$\phi _o = 0.553$ and 
$\phi _o = 0.550$. (2) Measurement imprecisions and external factors inherent to the experiment may not be completely modelled in the numerical simulations. For instance, free surface perturbations could be induced at the moment of gate removal and/or rapid collapse of the loose column. Such perturbations may increase the pore pressure and partially fluidize the granular column. Reference Rondon, Pouliquen and AussillousRondon et al. (2011) claimed that wall effects are negligible; however, the fluid flowing through the porous medium may have three-dimensional effects that cannot be captured by the 2-D numerical approach. (3) Limitations of the numerical model to fully reproduce the pore pressure feedback mechanism for a wide range of concentrations using a single formulation and set of numerical parameters.
3.2 Breaching process
3.2.1 Experimental set-up
Unlike loose shear failures, where the collapse is almost instantaneous, initially dense packings are distinguished by their slow collapse and the presence of negative pore pressures that stabilize the soil. At larger scales, it is also possible to observe the breaching process. As described in § 3.1, dilatancy induces negative pore pressure that reduces the shearing preventing a shear failure. However, near the front face, particles are released due to the expansion of the granular material causing the breach front to slowly regress.
In this section, the laboratory experimental data provided by Reference WeijWeij (2020) and Reference Alhaddad, Weij, Vanrhee and KeetelsAlhaddad et al. (2023) are used to evaluate the accuracy of the present numerical model to predict the breaching process. The experimental set-up of Reference WeijWeij (2020) and Reference Alhaddad, Weij, Vanrhee and KeetelsAlhaddad et al. (2023) consists of a tank with a height of 2 m, a length of 5.1 m and a width of 0.5 m (see figure 5). An impermeable removable gate, similar to the one used in Rondon's experiments, is placed at a distance of 2.5 m to divide the tank in two sections. A pump is located at the bottom right corner as shown in figure 5. The location of the pump prevents the reflection of the turbidity current. This pump takes out the mixture to a basin which is 4.5 m long, 1.25 m wide and 1.25 m high. In addition to this pump, a second pump is placed behind a 1 m high divider to maintain a constant water level by pumping clean water back into the reservoir. The left part of the tank is filled with sand prepared using the following procedure: the sand is placed horizontally layer by layer. In order to achieve a dense granular packing, layers are compacted using a vibrating needle. This process goes on until the height of the sand column is 1.47 m.
(a) Experimental set-up to study the breaching process. Image taken from Reference WeijWeij (2020). (b) Numerical set-up. In SedFoam, the ‘InletOutlet’ condition is written as ‘pressureInletOutletVelocity’ so the velocity is set to have a zero gradient condition when the flow leaves the domain, whereas the velocity assigned when the flow goes into the domain is based on the flux in the patch-normal direction.
Two different types of sand are considered in the laboratory experiments: (1) GEBA sand with a median diameter of 
$D_{50}=120\,\mathrm {\mu }$m, initial volume fraction of 
$\phi _o=0.585$ and an internal friction angle of 
$\mu _s=35.8^{\circ }$ and (2) D9 sand with a median diameter of 
$D_{50}=330\,\mathrm {\mu }$m, initial volume fraction of 
$\phi _o=0.570$ and an internal friction angle of 
$\mu _s=40.1^{\circ }$. For the numerical simulations, we observed that the size of the mesh elements should be at least 
$\Delta x =2.5$ mm (a uniform grid is considered). Larger grid sizes accelerate the breaching and produce more rounded shapes at the corner of the top right column. Although the mesh convergence is not completely reached, the computational cost becomes too expensive for finer meshes with little effect on the results in terms of morphology and wall velocity (a mesh convergence study has been carried out, the results of which are available in the supplementary material). Thus, the choice of 
$\Delta x =2.5$ mm seems reasonable to study the problem without compromising the accuracy of the solution. As shown in § 3.1.3 and suggested by Reference WeijWeij (2020), the critical volume fraction 
$\phi _c$ controls the amount of dilation, hence, the velocity of the receding wall. Accordingly to Reference WeijWeij (2020), the critical volume fraction is chosen be the same as the sand concentration during breaching, just before it is released from the breach face. In this case, 
$\phi _c=0.545$ and 
$\phi _c=0.56$ are chosen for the GEBA and D9 sands, respectively. Finally, the Engelund coefficients (linked to the permeability of the soil) are chosen as 
$\alpha _E=1200$ and 
$\beta _E=3.6$ for the GEBA sand and 
$\alpha _E=900$ and 
$\beta _E=1.8$ for the D9 sand. These values have been calibrated to match the experimental wall velocity (the horizontal velocity at which the steep slope moves due to the breaching process). Rheological and numerical parameters are summarized in table 5. Numerical results on the Rondon experiment (see § 3.1) evidenced small differences between 
$K_1=4$ and 
$K_1=40$. In this section we take 
$K_1=4$ because most of the works in the literature adopt dilatancy factors with the same order of magnitude.
Geometric, rheological and numerical parameters used to reproduce Reference WeijWeij (2020).
At this point, its is worth noting that (2.3) and (2.4) neglect the contribution of the turbulent suspension term. One may argue that at the breach face a turbidity current may be triggered, in which case, turbulent suspension/dispersion could significantly contribute to the erosion and this would need to be modelled through a turbulent suspension term in the momentum equation. However, this process takes place at a length scale which is of the order of the grain scale, smaller than the grid size of our numerical simulations; therefore, it cannot be resolved by the two-phase flow model in the experimental breaching configuration. In addition to the smaller grid sizes at the interface required to capture the transition from regions dominated by the granular rheology to the areas where the turbulent suspension term is dominant, other missing multiphase forces in our model, such as the lift force near the interface, may change the turbulent kinetic energy and affect the transition towards the turbulence suspension. Numerical tests with the turbulence term switched on have been performed showing an overestimated wall velocity as well as a complete erosion of the bottom deposit. The rapid erosion is probably due to the inability of the model to capture the correct turbulent suspension term. This term is proportional to the gradient of the solid volume fraction which becomes significantly high at the interface, where the volume fraction changes from 
${\sim }0$ to 
${\sim }0.6$ over a very short distance. As the capability of the two-phase model to simulate turbulent erosion by turbidity currents has not been checked even in simpler idealized situations, this would deserve further investigation that is beyond the scope of the present paper.
3.2.2 Dynamics of breaching
The removal of the gate triggers an initial creeping phase. The low mobility of the column is associated with the negative pore pressure measured with the sensors. Depending on the type of sand, different dynamics are observed.
In figure 6 we observe the evolution of the column morphology for the GEBA sand. The velocity of the receding front is well captured by the numerical model; however, slight discrepancies are found near the right top corner, where the numerical model predicts a rounded shape rather than a sharp corner observed in the experiments. The final shape consists of a deposit with a slope around 
$20^{\circ }$ modelled by SedFoam with remarkable agreement.
Comparison of the morphology between the experiments and the numerical simulations for the GEBA sand.
A key feature of the breaching process is the turbidity current formed near the front face as the particles are slowly released and pulled down by gravity. Such currents are illustrated in figure 7(a) and become less intense with time. Some particles released from the breaching front settle and contribute to form the deposit at the bottom. Although the main physics of turbidity currents are well captured, the lack of the turbulent suspension term may introduce inaccuracies in our numerical results. The presence of the deposit reduces the height of the breaching face and, by the end of the experiments and the numerical simulation, the deposit adopts a triangular shape with a slope milder than the friction angle in which the material is at rest with no turbidity current.
(a) Solid phase velocity field and (b) pore pressure field extracted in the numerical simulations for the GEBA sand. A grey line is included to illustrate the zero pressure isoline. (c) Comparison of the excess of pore pressure (
$p^f$) evolution within the granular column between the experiments and the numerical simulations for the GEBA sand.
The negative pore pressure plays the role of stabilizing the column and contributes to delay the granular flow. Figure 7(b) shows a pore pressure map very similar to the Rondon granular collapse studied in § 3.1. The negative pore pressure is significantly higher in the top right region of the granular column, where the material is more likely to fail along a shear plane. On the one hand, the negative pore pressure is responsible for a higher shear strength that keeps the granular column as a whole without shear failures. On the other hand, the dense material at the top right area is expanding which enhances the turbidity currents previously discussed. As the pore pressure dissipates, the shear resistance is reduced and the shear stress may reach the yield point forming a shear failure plane. This situation was reported in a few experiments (Reference WeijWeij, 2020) where occasional shear slide failures were observed during the breaching process. When slides occur, small drops in the excess of pore pressure are measured by the pressure sensors near the shear plane. In particular, the pressure drops located at 
$t \approx 80$ s and 
$t \approx 125$ s could have caused two successive minor slides in the experiments as shown in figure 7(c). These slides erode the front face much faster than the breaching process. Nonetheless, no sliding failures were observed in the numerical situations for the GEBA sand where the receding front face was mainly conducted by the breaching process. In the numerical simulations, figure 7 shows a pore pressure drop right after gate removal followed by a continuous dissipation with a time scale and dissipation rate comparable with the experimental ones.
The experiment and numerical simulations using the D9 sand show a similar trend, yet slightly different dynamics is observed in the experiments. The larger permeability of the D9 sand (grains roughly three times larger than those of the GEBA sand and lower initial volume fraction) leads to a faster breaching process. In figure 8(a) the evolution of the deposit morphology is well predicted by the numerical model in terms of shape and time scale. Despite the good agreement observed in figure 8(a), in the experiments, the initial breaching process is followed by a slow sliding failure. This feature is not captured by the numerical model. The dual-mode failure, including the breaching process and punctual shear failure, is evidenced in figure 8(b). Initially, the negative pore pressure builds up within the granular material stabilizing the whole granular column. Then, the progressive dilation of the medium is accompanied by a reduction of the pore pressure. At 
$t \approx 18$ s, the negative pore pressure is no longer capable of counteracting the shear forces; therefore, a wedge close to the breaching face starts sliding down. Simultaneously, the sliding region increases the shear rate significantly, especially near the failure line, which is accompanied by a negative pore pressure build-up (second minimum point in figure 8b). In contrast to the experimental measurements, the numerical approach is unable to reproduce this particular dual-mode failure.
(a) Comparison of the morphology and (b) the pore pressure (
$p^f$) evolution within the granular column between the experiments and the numerical simulations for the D9 sand.
4. Conclusion
The present article studied the sliding and breaching failures for different granular columns. In a first series of numerical tests, the model was able to reproduce the experimental data of Reference Rondon, Pouliquen and AussillousRondon et al. (2011) with a reasonable agreement. We showed that dilatancy effects are crucial to reproduce the rapid collapse of initially loose columns and the low mobility of initially dense columns. Additionally, the plastic effects embedded in the present dilatancy model are a key feature to predict the consequent pore pressure feedback mechanism signified by negative pore pressure for the dense case and positive pore pressure in the loose columns. The results presented in the sensitivity analysis suggest that dilatancy effects intensify with large 
$K_1$ values or with the choice of 
$\phi _c$ values far from the actual volume fraction. Concerning the basal pore pressure, little sensitivity to the dilatancy effects (
$K_1$ and 
$\phi _c$) is found for the initially dense column. However, the influence of dilatancy effects is much stronger for the loose case. Indeed, large 
$K_1$ or 
$\phi _c$ values lead to a partially fluidized bed. In such cases, some regions of the column exhibit high positive pore pressures that counterbalance the gravitational forces and the mixture flows easily until the pore pressure has dissipated so the shear strength builds up again. Although the model provided fairly good agreement with the experiments capturing the dilatancy effects at the grain scale through the dilatancy model proposed by Reference Montellà, Chauchat, Chareyre, Bonamy and HsuMontellà et al. (2021) and the expansion/contraction of the granular medium governed by 
$\phi (I)$ (see (2.17)), improvements of the dilatancy model could reduce and explain the current discrepancies, in particular the underpredicted positive pore pressures observed during the loose collapse. On the one hand, other dilatancy laws beyond the linear expression given by (2.15) should be explored. On the other hand, non-local dynamics should be introduced in the model to take into account the system size, which, in turn, has an important effect on the level of stresses as reported by Reference Athani, Metzger, Forterre and MariAthani, Metzger, Forterre, and Mari (2022).
For the second experimental configuration, the breaching process of Reference WeijWeij (2020) was numerically reproduced with very good agreement. The present numerical model was able to reproduce the breaching process with great success in terms of morphology and a reasonably good prediction of the pore pressure measured within the granular medium. The experiments of Reference WeijWeij (2020) evidenced that in some cases, sliding failures occur in addition to the breaching process. Although Reference WeijWeij (2020) did not report sliding failures in most of the experiments, they suggested that dual-mode failure is rather common in nature. Our numerical simulation mainly reproduces breaching but a different set of parameters may predict features observed in slide failures. It seems, therefore, reasonable to include the transition of slope failures from breaching to slides in future research in order to make an accurate prediction of the slope failure modes. As a main conclusion, it has been shown that a two-phase flow numerical model including dilatancy effects is able to reproduce the various failure modes of underwater particles collapse, breaching and slide failures, and their sensitivities to the initial dense volume fractions as well as their dynamics and morphology of the final deposits. The results of the present model could be extended in the future to investigate the influence of several geometrical and physical parameters on the dynamics of granular collapse without the need of experiments. More particularly, the dredging industry could benefit from the present two-phase model to predict flow slides during excavation processes, thereby further action may be considered to mitigate the erosion process. Finally, the present model could shed some light on the transition between breaching and liquefaction flow slides.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/flo.2023.23. The numerical simulations corresponding to the Reference Rondon, Pouliquen and AussillousRondon et al. (2011) configuration are available on GitHub (https://github.com/SedFoam/sedfoam/tree/master/tutorials/laminar/2DCollapse). Results, post-process scripts and numerical set-ups for the breaching cases (Reference WeijWeij, 2020) are available for download on Zenodo (https://zenodo.org/record/8116463).
Acknowledgements
The authors also would like to acknowledge the Strategic Environmental Research and Development Program, USA (MR20-1478). We would also like to acknowledge the support from the US National Science Foundation (CMMI-2050854). Most of the computations presented in this paper were performed using the GENCI infrastructure under Allocations A0060107567 and A0080107567 and the GRICAD infrastructure. Finally, we are also grateful to the developers involved in OpenFOAM.
Funding statement
This work received funds from the Strategic Environmental Research and Development Program, USA, through grant number MR20-1478.
Declaration of interests
The authors declare no conflict of interest.
Open access
