The compressible granular collapse in a fluid as a continuum: validity of a Navier-Stokes model with $\mu(J)$-$\phi(J)$-rheology

The incompressible $\mu(I)$-rheology has been used to study subaerial granular flows with remarkable success. For subaquatic granular flows, drag between grains and the pore fluid is substantially higher and the physical behaviour is more complex. High drag forces constrain the rearrangement of grains and dilatancy, leading to a considerable build-up of pore pressure. Its transient and dynamic description is the key to modelling subaquatic granular flows but out of the scope of incompressible models. In this work, we advance from the incompressible $\mu(I)$-rheology to the compressible $\mu(J)$-$\phi(J)$-rheology to account for pore pressure, dilatancy, and the scaling laws under subaquatic conditions. The model is supplemented with critical state theory to yield the correct properties in the quasi-static limit. The pore fluid is described by an additional set of conservation equations and the interaction with grains is described by a drag model. This new implementation enables us to include most of the physical processes relevant for submerged granular flows in a highly transparent manner. Both, the incompressible and compressible rheologies are implemented into OpenFOAM and various simulations at low and high Stokes numbers are conducted with both frameworks. We found a good agreement of the $\mu(J)$-$\phi(J)$-rheology with low Stokes number experiments, that incompressible models fail to describe. The combination of granular rheology, pore pressure, and drag model leads to complex phenomena such as apparent cohesion, remoulding, hydroplaning, and turbidity currents. The simulations give remarkable insights into these phenomena and increase our understanding of subaquatic mass transports.

to simulate these phenomena with Navier-Stokes type models. The vast amount of studies relies on the ( )-rheology and its derivatives. The core of the ( )-rheology is the Drucker-Prager yield criterion (Drucker & Prager 1952;Rauter et al. 2020) and the recognition that the friction coefficient is solely a function of the inertial number (GDR MiDi 2004;Jop et al. 2006). Further studies found a similar correlation between the inertial number and the packing density (Forterre & Pouliquen 2008).
A similar scaling was found in granular flows with low Stokes numbers (see Eq. (2.31)). The Stokes number is related to the ratio between inertia and drag force on a particle and thus describes the influence of ambient fluid on the granular flow dynamics (e.g. Finlay 2001). Small Stokes numbers indicate a strong influence of the pore fluid on the particles, and hence also on the landslide dynamics. In this regime, the viscous number replaces the inertial number as a control parameter for the friction coefficient and the packing density , forming the so-called ( )-( )-rheology (Boyer et al. 2011). Furthermore, excess pore pressure can be remarkably high under these conditions and it is imperative to explicitly consider it in numerical simulations. High drag forces and respectively small Stokes numbers are usually related to small particles. They are virtually omnipresent in geophysical flows: submarine landslides (Kim et al. 2019), turbidity currents (Heerema et al. 2020), powder snow avalanches (Sovilla et al. 2015), and pyroclastic flows (Druitt 1998) can be dominated by fine grained components. It follows that a large portion of gravitational mass flows occurs at low Stokes numbers and a deeper understanding of the respective processes is relevant for many researchers.
Incompressible granular flow models have been applied in different forms to various problems in the last decade. Lagrée et al. (2011) were the first to conduct numerical simulations of subaerial granular collapses with the ( )-rheology and the finite volume method. Staron et al. (2012) used the same method to simulate silo outflows, and Domnik et al. Most of the mentioned applications rely on standard methods from computational fluid dynamics (CFD). This is reasonable, considering the similarity between the hydrodynamic (Navier-Stokes) equations and the granular flow equations. However, the pressure dependent and shear thinning viscosity associated with granular flows introduces considerable conceptual and numerical problems. The unconditional ill-posedness of an incompressible granular flow model with constant friction coefficient was described by Schaeffer (1987) and the partial ill-posedness of the ( )-rheology by Barker et al. (2015). By carefully tuning the respective relations, Barker & Gray (2017) were able to regularize the ( )-rheology for all but very high inertia numbers. Barker et al. (2017) described a well-posed compressible rheology, incorporating the ( )-rheology as a special case.
Another pitfall of granular rheologies is the concept of effective pressure. When pore pressure is considerably high (i.e. at low Stokes numbers), it is imperative to distinguish between effective pressure and total pressure (first described by Terzaghi 1925). Effective pressure represents normal forces in the grain skeleton that have a stabilizing effect, in contrast to pore pressure which has no stabilizing effect. This has shown to be a major issue, as pore pressure and consequently the effective pressure, react very sensitively to the packing density and dilatancy (Rondon et al. 2011).
The volume-of-fluid method, which is also used in this work, allows to track the slide as a single component but also as a mixture of multiple phases (grains and pore fluid). Components are defined in here as objects (e.g. the landslide) that completely cover a bounded region in space without mixing with other components (e.g. the ambient fluid), see Fig. 3. The tracking becomes a purely geometric problem (see e.g. Roenby et al. 2016, for a geometric interpretation). In contrast, phases (e.g. grains) are dispersed and mixed with other phases (e.g. pore fluid) to represent the dynamic bulk of the landslide, see Fig. 1.
The component-wise tracking is used in various landslide models (e.g. Lagrée et al. 2011;Domnik et al. 2013;Barker & Gray 2017). Components, i.e. the slide and the surrounding fluid, are immiscible and separated by a sharp interface. Usually, this also implies that the model is incompressible. The phase-wise tracking is commonly applied in chemical engineering (Gidaspow 1994;van Wachem 2000;Passalacqua & Fox 2011) and has lately been introduced to environmental engineering (e.g. Cheng et al. 2017;Chauchat et al. 2017;Si et al. 2018a). This approach allows to describe a variable mixture of grains and pore fluid that merges smoothly into the ambient fluid. The description of the pore fluid as an individual phase enables the model to decouple effective pressure from pore pressure, which is imperative in many flow configurations, e.g. for low Stokes numbers.
In this work, a two-component and a two-phase Navier-Stokes type model are applied to granular flows. Both models are implemented into the open-source toolkit OpenFOAM (Weller et al. 1998;Rusche 2002;OpenCFD Ltd. 2004), using the volume-of-fluid method for component-and phase-wise tracking (see section 2). Subaerial (Balmforth & Kerswell 2005) and subaquatic granular collapses (Rondon et al. 2011) are simulated with both models and results are compared to the respective experiments and with each other.
We apply the ( )-( )-rheology to subaerial cases ( 1) and the ( )-( )-rheology to subaquatic cases ( 1). The two-component model applies simplified rheologies in form of the incompressible ( )and ( )-rheologies. The ( )and ( )-curves are merged into the particle pressure relation of Johnson & Jackson (1987) to achieve the correct quasi-static limits (Vescovi et al. 2013). This yields reasonable values for the packing density at rest which is imperative for granular collapses with static regions. In contrast to many previous works (e.g. Savage et al. 2014;von Boetticher et al. 2017;Si et al. 2018a), we renounce additional contributions to shear strength (e.g. cohesion) because we do not see any physical justification (e.g. electrostatic forces, capillary forces, cementing) in the investigated cases. We apply a very transparent and simple model, focusing on the relevant physical processes and achieve a remarkable accuracy, especially in comparison to more complex models (e.g. Si et al. 2018a;Baumgarten & Kamrin 2019). Further, it is shown that various experimental setups with different initial packing densities can be simulated with the same constitutive parameters, whereas many previous attempts required individual parameters for different cases (e.g. Savage et al. 2014;Wang et al. 2017;Si et al. 2018a).
The paper is organised as follows: The multi-phase (section 2.1) and multi-component (section 2.2) models are introduced in section 2, including models for granular viscosity (section 2.3), granular particle pressure (sections 2.4 and 2.5) and drag (section 2.6). Results are shown and discussed in section 3 for a subaerial case and in section 4 for two subaquatic cases. A conclusion is drawn in section 5 and a summary is given in section 6. Furthermore, a thorough sensitivity analysis is provided in the appendix.

Methods
2.1. Two-phase landslide-model The two-phase model is based on the phase momentum and mass conservation equations (see e.g. Rusche 2002). The governing equations for the continuous fluid phase are given as and for the grains as Phase-fraction fields g and c , i.e. the phase volume over the total volume = , (2.5) describe the composition of the grain-fluid mixture, see Fig. 1 (the index indicates either c or g). The granular phase-fraction is identical with the packing density = g . Phase-fractions take values between zero and one and the sum of all phase-fractions yields one. The pore fluid is assumed to match the surrounding fluid and the respective phase-fraction c is therefore one outside the slide. This way, phase-fraction fields provide not only a mechanism to track the packing density of the slide, but also its geometry. Every phase moves with a unique velocity field , which is not divergence-free. This allows the mixture to change, yielding a variable packing density and thus bulk-compressibility, although phase densities g and c are constant. The volume weighted average velocity is divergence free, ∇ · = ∇ · g g + c c = 0, (2.6) which allows to use numerical methods for incompressible flow. The pore pressure (or shared pressure) is acting on all phases equally, while the grain phase experiences additional pressure due to force chains between particles, the so called effective pressure (or particle pressure) s , see Fig. 2. The effective pressure is a function of the packing density in this model and the balance between effective pressure and external pressure (e.g. overburden pressure) ensures realistic packing densities. The total pressure can be assembled as tot = + s .
(2.7) The deviatoric phase stress tensors are expressed as with phase viscosity , phase density and deviatoric phase strain rate tensor (2.9) Focus on Fluids articles must not exceed this page length The effective pressure s (red arrows) represents normal forces in the grain skeleton (black arrows). The pore-pressure (blue arrows) represents pressure that is equally shared by pore fluid and grains.
The viscosity of the pore fluid c is usually constant and the granular viscosity g is following from constitutive models like the ( )-rheology (see section 2.3). The total deviatoric stress tensor can be calculated as (2.10) The last terms in Eqs. (2.2) and (2.4) represent drag forces between phases and gc is the drag coefficient of the grains in the pore fluid. Lift and virtual mass forces are neglected in this work, because they play a minor role (Si et al. 2018a).
The granular viscosity g , the effective pressure s , and the drag coefficient gc represent interfaces to exchangeable sub-models, presented in sections 2.3-2.6.

Two-component landslide-model
Many two-phase systems can be substantially simplified by assuming that phases move together, i.e. that phase velocities are equal, (2.11) This fits very well to completely separated phases that are divided by a sharp interface (e.g. surface waves in water, Rauter et al. 2021) but also systems of mixed phases (e.g. grains and fluid) can be handled to some extent (e.g. Lagrée et al. 2011). The phase momentum conservation equations (2.2) and (2.4) can be combined into a single momentum conservation equation and the system takes the form of the ordinary Navier-Stokes Equations with variable fluid properties (see e.g. Rusche 2002), Values between zero and one are not intended by this method and only appear due to numerical reasons, i.e. the discretisation of the discontinuous field (see section 2.7). In here, two component indicator functions are used, one for the ambient fluid component, c , and one for the slide component, s (see Fig. 3). Evolution equations for component indicator functions can be derived from mass conservation equations as The definition of components is straight forward for completely separated phases, where components can be matched with phases, e.g. water and air. The definition of the slide component, on the other hand, is not unambiguous, as it consists of a variable mixture of grains and pore fluid. A boundary of the slide component can, for example, be found by defining a limit for the packing density (e.g. 50% of the average packing density). Further, a constant reference packing density has to be determined, which is assigned to the whole slide component. The density of the slide component follows as   with the component viscosity and the deviatoric shear rate tensor S. Note that the deviatoric shear rate tensor S matches the shear rate tensor D, because the volume weighted averaged velocity field is divergence free, (2.20) The viscosity of the ambient fluid c is usually constant and the viscosity of the slide region s is following from granular rheology, see section 2.3.
2.3. Rheology 2.3.1. Unifying rheologies Most granular rheologies (e.g. the ( )-rheology) are defined in terms of the total deviatoric stress tensor in the slide component T s . This has to be accounted for and corrected in the two-phase model if the same viscosity model is used in both models. Similar to Eq. (2.16), component viscosities can be related to phase viscosities as The contribution of the granular phase to stresses is assumed to be much higher than the contribution of the pore fluid, g g S g ≫ (1 − ) c c S c . Further, by neglecting the mass of the pore fluid, s ≈ g , it follows that kinematic viscosities have to be similar in both models, s ≈ g . (2.23) Alternatively, one can match the dynamic viscosities s s and g g if the factor g is removed from the viscous term in Eq. (2.4). Note, that this assumptions are fairly accurate for subaerial granular flows but questionable for subaquatic granular flows. However, multiphase and multi-component models differ substantially under subaquatic conditions and a unification is not possible.

Drucker-Prager plasticity model
An important characteristic of granular materials is the pressure dependent shear stress, described by the Drucker-Prager yield criterion (Drucker & Prager 1952). Schaeffer (1987) was the first to include granular friction in the Navier-Stokes equations by expressing the Drucker-Prager yield criterion in terms of the shear rate tensor and the pressure, where the norm of a tensor A is defined as The friction coefficient is constant and a material parameter in the first model by Schaeffer (1987 where is the characteristic height of the investigated case, was found to give a good estimate for a reasonable viscosity cut-off. Notably, the Drucker-Prager yield surface leads to an ill-posed model (Schaeffer 1987) and the truncation of the viscosity is not sufficient for a regularization.
Schaeffer (1987) did not distinguish between effective and total pressure in Eq. (2.26), limiting the applications of his model substantially. We will explicitly consider effective pressure in Eqs. (2.26) and (2.27) using Eq. (2.34) or (2.36) in the two-component model and Eq. (2.37), (2.40), or (2.43) in the two-phase model to avoid such limitations.

( )-rheology
The ( )-rheology (GDR MiDi 2004;Jop et al. 2006;Forterre & Pouliquen 2008) states that the friction coefficient is not constant in dense, dry, granular flows but rather a function of the inertial number . The inertial number is defined as the ratio between the typical time scale for microscopic rearrangements of grains with diameter , micro = g / s , and the macroscopic time scale of the deformation, macro = 1/2 S −1 , In the two-phase model, the shear rate S is replaced with the deviatoric shear rate of grains S g . Various approaches have been proposed for the ( )-curve, in here we apply the classic relation, given as The two-component model is limited in considering pore pressure and dilatancy effects because the packing density is not described by this model. The effective pressure can only be reconstructed from total pressure tot and various assumptions. The simplest model assumes that the pore pressure is negligibly small, leading to s ≈ tot . (2.34) This assumption is reasonable for subaerial granular flows and has been applied to such by e.g. Lagrée et al. (2011)  Here, 0 is the position of the free water surface, where the total pressure is supposed to be zero. For a variable and non-horizontal free water surface, common in e.g. landslidetsunamis, this concept is complicated substantially, and to the authors knowledge, not applied. Furthermore, excess pore pressure, which is common in low Stokes number flows, is out of the scope for this model.

2.5.
Effective pressure in the two-phase model 2.5.1. Critical state theory The structure of the two-phase model allows us to include the packing density in the effective pressure equation. Critical state theory (Roscoe et al. 1958;Roscoe 1970;Schofield & Wroth 1968) was the first model to describe the relationship between the effective pressure and the packing density. The critical state is defined as a state of constant packing density and constant shear stress, which is reached after a certain amount of shearing of an initially dense or loose sample. The packing density in this state, called critical packing density crit , is a function of the effective pressure s . This function can be inverted to get the effective pressure as a function of the critical packing density. It is further assumed that the flow is in its critical state g = crit to get a model that is compatible with the governing equations. This assumption is reasonable for avalanches, slides, and other granular flows but questionable for the initial release and deposition. At small deformations, the packing density might be lower (underconsolidated) or higher (overconsolidated) than the critical packing density and the effective pressure model will over-or underestimate the effective pressure.
A popular relation for the effective pressure (the so-called critical state line) has been described by Johnson & Jackson (1987); Johnson et al. (1990) as where rlp is the random loose packing density in critical state, rcp the random close packing density in critical state and a scaling parameter. The scaling parameter can be interpreted as the effective pressure at the packing density 1 2 rcp + rlp . Note that we apply a simplified version of the original relation, similar to Vescovi et al. (2013). Packing densities above rcp are not valid and avoided by the asymptote of the effective pressure at rcp . If packing densities higher or equal rcp appear in simulations, they should be terminated and restarted with refined numerical parameters (e.g. time step duration).

( )-relation
Equation ( where max and Δ are material parameters. This relation can be transformed into a model for the effective pressure by introducing the inertial number , This relation has two substantial problems: For S g = 0 it yields s = 0 and for g = 0 it yields s ≠ 0, which causes numerical problems and unrealistic results. The first problem is addressed by superposing Eq. (2.39) with the quasi-static relation (2.37), similar to Vescovi et al. (2013). The second problem is solved by multiplying Eq. (2.39) with the normalized packing density g / , which ensures that the pressure vanishes for g = 0. The normalization with the reference packing density ensures that parameters ( max , Δ ) will be similar to the original equation. Further, to reduce the number of material parameters, we set the maximum packing density in the ( )-relation equal to the random close packing density rcp . The final relation reads . Notably, the inertial number is a function of only the packing density and the shear rate, = g , S g , because the effective pressure is calculated as function of the packing density. The same follows for the friction coefficient = g , S g and the deviatoric stress tensor T g = g , S g . This highlights that the two-phase model implements a density-dependent rheology, rather than a pressure-dependent rheology.
It should be noted that there are various possibilities to combine critical state theory and Notably, Boyer et al. (2011) emphasised that m is not matching the random close packing density rcp ≈ 0.63 but rather a value close to 0.585. This leads to substantial problems for large values of g as the relation is only valid for g < m = 0.585 or S g = 0. In other words, shearing is only possible for g < m . We solve this issue by allowing a creeping shear rate of 0 at packing densities above m . Further and as before, we superpose the relation with the quasi-static relation of Johnson & Jackson (1987) to yield the correct asymptotic values for S g → 0 known from critical state theory. The final relation reads (2.44) The respective relation is shown in Fig with the grain diameter as the only parameter. This relation is supposed to be valid for small relative velocities and densely packed granular material. It has been modified to account for higher relative velocities (Ergun 1952) and lower packing densities (Gidaspow 1994), however, which is not relevant for the investigated configurations (see Si et al. (2018a) for an application of the extended relation). This relation is visualized in Fig. 6a for various diameters and packing densities. ( 2.47) The permeability is visualized in Fig. 6b. In a similar manner, the drag coefficient can be calculated as if the intrinsic permeability of the granular material is known.

Two-component landslide-model
The two-component model is based on the solver multiphaseInterFoam, using the PISOalgorithm (Issa 1986) and interpolations following Rhie & Chow (1983) to solve the coupled system of pressure and velocity. First, an updated velocity field is calculated without the contribution of pressure. The predicted velocity field is later corrected to be divergence-free and the pressure follows from the required correction. Finally, all other fields, e.g. the phase indicator functions, are updated. This procedure is repeated in each time step. Components (slide and ambient air or water) are divided by an interface which is supposed to be sharp. However, the interface is often smeared by numerical diffusion. To keep the interface between components sharp, the relative velocity between phases ij , which was previously eliminated from the system, is reintroduced in Eq. (2.15), Eq. (2.49) is finally solved using the MULES algorithm (Multidimensional Universal Limiter with Explicit Solution) (Weller 2008). This scheme limits the interface compression term (i.e. the term containing ) to avoid over-( > 1) and undershoots ( < 0) of the component indicator fields.
There is no conservation equation for the relative velocity in the two-component model and it has to be reconstructed from assumptions. Two methods are known to construct the relative velocity for granular flows. Barker et al. (2021) suggest to construct the relative velocity for granular flows from physical effects such as segregation and settling. The relative velocity follows as the terminal velocity of spheres in the surrounding fluid under the influence of gravity. Alternatively, one can construct a velocity field that is normal to the interface and of the same magnitude as the average velocity , This method has a maximum sharpening effect (Weller 2008) and is thus also applied in this work.

Two-phase landslide-model
The two-phase model is based on the solver multiphaseEulerFoam. The system of pressure and average velocity is solved with the same concept as in the two-component solver. The velocity fields for all phases are first predicted without contributions from pore pressure , but including effective pressure s . The average velocity is then corrected to be divergencefree and the pore pressure follows from the required correction. In a further step, the velocity correction is applied to phase velocities. The solution procedure is described in depth by Rusche (2002). The interface compression term is not required in this model because settling and segregation is directly simulated and counteracting numerical diffusion. The implementation of the effective pressure term is taken from SedFoam 2.0 (Chauchat et al. 2017).

Time stepping
The numerical solution of transport equations is subject to limitations that pose restrictions on the solution method. One of these limitations is known as the Courant-Friedrichs-Lewy (CFL) condition and enforced by limiting the CFL number. In convection dominated problems, the CFL number is defined as the ratio of the time step duration Δ and the cell convection time Δ / , i.e. the time required for a particle to pass a cell with size Δ , For the stability of e.g. the forward Euler method, it is required, that the convection time is smaller than the time step duration, and similar limits exist for other explicit methods. This limitation has to be enforced by choosing the time step duration Δ according to mesh size and flow velocity. However, Eq. (2.51) is only valid for convection dominated problems. In the case of granular flows, the viscosity term is dominating over all other terms. Therefore, the viscosity has to be considered in the calculation of the CFL number and the time step duration. The respective definition, ignoring the contribution of convection follows as This relation is imperative for stability of explicit and semi-implicit Navier-Stokes solvers when viscous forces are dominating. The squared cell size in the denominator and the high viscosity introduce very strict limitations on the time step, making computations very expensive. Note that simplified relations for the one-dimensional case are given in here. The full multi-dimensional conditions for arbitrary finite volume cells can be found in Rauter et al. .

Subaerial granular collapse (Balmforth & Kerswell 2005)
As a first test of the numerical models, we simulate the granular collapse experiments of Balmforth & Kerswell (2005) under subaerial conditions. A sketch of the experiment is shown in Fig. 7. The experiment was conducted between two parallel, smooth walls and the setup is approximated as a 2D granular collapse. Balmforth & Kerswell (2005) conducted multiple experiments with different geometries, in here we focus on the experiments with an aspect ratio of / = 1/2, but similar results have been obtained for other aspect ratios. In theory, both, the two-component and the two-phase model should be equally capable of simulating this case because pore pressure plays a minor role. Most parameters, such as density, quasi-static friction coefficient, and particle diameter are reported by Balmforth & Kerswell (2005). The missing parameters are completed with data from the literature. Notably, the experiments are conducted on a smooth surface, which was incorporated in simulations by switching to a constant friction coefficient wall at smooth surfaces. This modification is simple in the finite volume method because stresses are calculated on cell faces before their divergence is calculated as a sum over faces. The Stokes number is estimated to be of order 10 3 (with S = 10 s −1 ) for this experiments and the ( )-( )-rheology is chosen to describe friction and effective pressure. Parameters . These parameters are therefore optimized to the subaquatic case (section 4), where extended measurements are available, and applied to this case without further modification. The average packing density is assumed to be = 0.6 following the critical state line at this pressure level. The applied pressure equation is visualized in Fig. 4. From the height = 0.1 m the required viscosity threshold max can be estimated following Eq. (2.28) to be of order 1 m 2 s −1 . This estimation was validated with a sensitivity analysis (see appendix B.1). The final set of parameters is given in Tab. 1. Regular meshes of square cells are used to cover a simulation domain of 0.5 × 0.2 m, which was sufficient to have no artificial influences from boundaries. Standard boundary conditions are applied at walls (zero velocity, zero pressure gradient) and the permeable top (zero velocity gradient, zero pressure). Multiple mesh resolutions were applied to investigate the influence of the grid resolution on the results (see appendix B.2). The time stepping was investigated with a similar approach, modifying the limit for CFL diff max between 1 and 1000 (depending on model and solver mode, see appendix B.3). In the following, the CFL-number is limited to 1 and the cell size set to 0.0017 m, which showed to be sufficient to achieve converged and mesh independent results.

Two-component model
The component indicator for the slide component s is initialized to 1 within the square that forms the initial granular column. We assume that hydrostatic pore pressure is negligible (< 2 Pa) and therefore apply Eq. (2.34) to calculate the effective pressure.
The simulation covering a simulation duration of 0.8 s took 6.9 h on eight cores of LEO4 (High Performance Cluster from the University of Innsbruck, consisting of Intel Xeon (Broadwell) compute cores). The total pressure, which is assumed to match the effective pressure, is shown for three-time steps in Fig. 8, alongside the final pile in the experiment. The continuous black line shows the sharp free surface, assumed to be located at s = 0.5. Furthermore, the velocity field is shown as arrows. The collapse takes about 0.4 s and the pile remains in its final shape for the rest of the simulation. The two-component model matches the experiment well, however, the volume of the final pile is slightly underestimated. Results are very robust in terms of mesh refinement or coarsening (see appendix B.2) and

Two-phase model
The two-phase model uses the same parameters as the two-component model, including numerical parameters, such as viscosity threshold and CFL limit. The phase-fraction g was initialized such that effective pressure is in balance with lithostatic vertical stresses, yielding an initial mean phase-fraction of g = 0.608. This procedure ensures that there will be no dilatancy or compaction in stable regions of the pile.
The simulation took 9.1 h under the same conditions as the two-component simulation. A stronger mesh dependency is observed for this model, however, the runout is converging for fine meshes (see appendix B.2). The pore pressure and the effective pressure following the extended ( )-theory are shown for three time steps in Fig. 9, alongside the final pile shape in the experiment. The continuous black line indicates the position of the free surface, assumed to be located at s = 0.25. The average velocity is shown as arrows in Figs. 9a-c, the relative velocity of air with respect to grains in Figs. 9d-f. The relative velocity in the initial phase is considerably high, indicating an inflow of air into the bulk and thus dilatancy. The two-phase model matches the experiment exceptionally well and the dilatancy in the experiment is matched by the simulation to a high degree. Note that the effective pressure at rest is directly linked to the packing density which can be qualitatively estimated from Fig. 9f. 3.3. Discussion and comparison Both models performed well at simulating the subaerial granular collapse. This is in line with previous results of e.g. Lagrée et al. (2011) or Savage et al. (2014. The effective pressure and the total pressure are fairly similar, because excess pore pressure is dissipating quickly through dilatancy or compaction. The magnitude of pore pressure in the two-phase model is smaller than 8 Pa and thus less than 1 % of the effective pressure, validating the assumption of neglectable pore pressure. The runout is similar in both models, the front is slightly elongated in the two-phase model. Further, the two-phase model shows a better match with the experiment at the upper end of the final slope. Both of these minor differences can be attributed to dilatancy effects. The two-component model is intrinsically not able to capture this process. Two mechanisms for dilatancy can be observed in the two-phase model. Firstly, the average effective pressure in the slide is reduced as it is spreading out and the packing density decreases proportionally to the effective pressure, as prescribed by the critical state line. Secondly, shearing can reduce the packing density well below its critical packing density due to the dynamic contribution of the ( )-theory to effective pressure. The loosely packed slide will not return to the critical packing density after shearing but remain in a loose state, forming a hysteresis. The granular material is able to remain in a loose state because the deviatoric stress tensor counteracts one-dimensional settling deformations (known as oedometric compression in soil mechanics). Furthermore, the granular column may have been overconsolidated in the experiment, however, this was not incorporated in the model due to the initialisation in critical state. Dilatancy is rather unimportant under subaerial conditions, as it does not imply changes in rheology or flow dynamics. Therefore, the two-component model is well suited for subaerial granular collapses, where pore pressure is negligibly small and the Stokes number is well above one.
The reduced friction at the smooth basal surface has a small but noticeable effect on the final pile shape. The runout is longer when incorporating the smooth surface and matches the experiment better. Previous works (e.g. Savage et al. 2014) ignored the smooth bottom of the experiment and still obtained accurate final pile shapes by using a constant friction coefficient. The increased friction of the ( )-rheology (in comparison to a constant quasi-static friction coefficient) compensates for the reduced basal friction quite exactly (see appendix B.4).
The two-component model is less sensitive to grid resolution than the two-phase model (see appendix B.2) but more sensitive to the time step duration (see appendix B.3). At the same resolution, both models require roughly the same computational resources and no model shows a substantial advantage in this regard.
It is important to carefully choose the time step duration, as it can have drastic influences on simulation results. Generally, CFL diff has to be limited to one to guarantee satisfying results, while some cases and solver settings allow higher CFL diff numbers. This limitation is much stronger than the traditional CFL criterion and CFL conv is roughly 0.001. Notably, the time step duration is constant in simulations, Δ ≈ 3 · 10 −6 s, because the constant maximum viscosity max in stable regions and the constant cell size Δ controlled the time stepping. 2 · 10 −5 m 2 s −1 (about 10 times higher than water), leading to a very low permeability of ≈ 10 −10 m 2 following Eq. (2.47). Early tests revealed that the two-phase model reacts very sensitively to the critical state line parameters rlp , rcp , and . Parameters from literature (e.g. the critical state line applied by Si et al. 2018a) lead to unrealistic granular pressures at g = 0.60 and could thus not be applied. We set the limiting packing densities to rlp = 0.53 and rcp = 0.63 to allow initial average packing densities between 0.55 and 0.61. The scaling variable was found by matching the peak pore pressure in the dense simulation with the respective measurement (see Fig. 14). The total set of parameters used for both cases is shown in Tab. 2.
Regular meshes of square cells with size 0.0005 m are applied, covering a simulation domain of 0.15 m × 0.105 m (dense case) and 0.25 m × 0.105 m (loose case). The CFL number CFL diff is limited to 10 in order to keep computation times to a reasonable level. A sensitivity study was conducted to proof convergence at this grid size (see appendix B.2) and CFL diff number (see appendix B.3).

Two-component model -dense case
The hydrostatic pore pressure is high under subaquatic conditions and the two-component model applies Eq. (2.36) to consider its influence on the effective pressure. All parameters are taken from Tab. 2. The evolution of the slide geometry, the effective pressure and the velocity are shown in Fig. 11, alongside the final experimental pile shape. The final pile shape of the model corresponds roughly to the experiment. The velocity, on the other hand, is roughly corresponding to the loose case, and the collapse is completed after 1 s, whereas the dense experiment took more than 30 s. The simulation and its failure mechanism are similar to the subaerial case where the free unsupported side of the pile is collapsing until reaching a stable slope inclination. Notably, neither the dense nor the loose experiment showed such a failure mechanism (see Fig. 15). No excess pore pressure is included in this model and a hypothetical pressure sensor at the bottom of the column would measure constantly 0 Pa as indicated in Fig. 14.

Two-component model -loose case
The two-component model provides only a few and ineffective possibilities to consider variations of the packing density. To simulate the loose granular collapse with this model, the average packing density is changed to = 0.55 and the bulk density correspondingly to s = 1825 kg m 3 . Further, the initial column geometry is changed as reported by Rondon et al. (2011). All other parameters match the dense case. Changing rheology parameters, e.g. 1 or 2 (as e.g. Wang et al. 2017) is technically possible but does not help in understanding the physical process or the influence of packing density.
The difference to the dense simulation is very small and thus not shown in here (see, e.g. Bouchut et al. 2017, for similar results). As before, the final pile shape is close to the dense experiment while the simulated velocity is close to the loose experiment. The runout is slightly longer as in the dense simulation because the loose column is slightly taller.

Two-phase model -dense case
The two-phase model allows us to explicitly consider variations in the initial packing density. The dense case is initialized with a homogeneous packing density of 0.60. The evolution of the dense granular column as simulated with the two-phase model is shown in Fig. 12, alongside three states of the experiment at = 3 s, 6 s and 30 s. The simulation covering a duration of 10 s, took 240 h on 8 cores of LEO4. The dense case is dominated by negative excess pore pressure (Fig. 12a-e), meaning that pore pressure within the slide is lower than outside. The effective pressure (Fig. 12f-j) is respectively higher, which increases the shear strength of the column. Initially, the shear strength is high enough to delay the collapse and to keep the column mostly stable. The pore pressure gradient leads to the suction of fluid into the column (Fig. 12g-h) and the granular material is dilating. Dilation reduces the effective pressure and allows the column to collapse. This happens first near the free surface on the unsupported side of the column, leading to a breaching-like flow of grains (Fig. 12g-h). Grains mix with fluid at the breaching edge, reducing packing density, effective pressure, and thus friction to very low values. The resulting mixture behaves like a small turbidity current and reaches long run-outs with shallow slopes, as visible in Fig. 12i-j. The zone of low particle pressure extends towards the centre of the column with time and further mobilisation occurs. At = 0.5 s, we can see the formation of a shear band. The grains above the shear band slide off, first as a triangular cohesive block (note the uniform velocity field in Fig. 12b), which disintegrates between = 1 s and = 3 s (Fig. 12i). The overall process is finished (i.e. end ) in the simulation after roughly 10 s, while the experiment took about 30 s. The final pile form and the failure mechanism match the experiment very well, which can be seen best in a comparison with the videos provided by Rondon et al. (2011), see Fig. 15. Further, a good match with the measured excess pore pressure is achieved, as shown in Fig. 14. The time scale and velocity of the collapse, on the other hand, differ substantially between simulation and experiment. Notably, pore pressure and effective pressure s do not sum up to the total vertical load, as a considerable fraction of the vertical load is transferred to the ground by viscous stresses.

Two-phase model -loose case
The simulation of the loose granular column uses the same parameters as the dense simulation. The packing density in the column is initialized homogeneously to = 0.55 and its height is increased as reported by Rondon et al. (2011). The simulation covering a duration of 6 s took 213 h on 8 cores of LEO4. As a result of the very loose packing, the effective shear strength is low and the column collapses rapidly and entirely, without any static regions. The pore pressure has to support the majority of the weight and is respectively high (Fig. 13a). The effective pressure increases at the rapidly flowing front, at = 0.25 s (Fig. 13g) due to the dynamic increase of effective pressure following the ( )-theory. The increase in effective pressure leads to a proportional increase in friction and the front is slowed down, Fig. 13h-i. Although the effective pressure is low in comparison to the dense case (four times lower), the friction is sufficient to bring the slide to a stop. The final slope inclination is shallow and the low quasi-static particle pressure is sufficient to support the slope, Fig. 13j. The packing density increases slightly during the collapse but the stability is mostly gained by reducing the slope inclination. The final pile shape matches the experiment very well, only a small amount of granular material forms a turbidity current that exceeds the runout of the experiment. The simulated velocity is higher than in the experiment but the difference is less severe than in the dense case. The simulated excess pore pressure differs remarkably from the measured excess pore pressure as shown in Fig. 14. Two stages can be observed in the simulated excess pore pressure history. First, the simulation shows a high peak of excess pore pressure, exceeding the highest experimental pore pressure by a factor of two. The peak dissipates quickly, as the slide and thus overburden pressure leave the region where the pore pressure sensor is installed. This first peak is not appearing in the experiment, where the highest pore pressure is reached in a flatter peak at a later point in time. In a second phase, starting at = 1 s, the pore pressure dissipates much slower. In this phase, the pore pressure dissipation is driven by compaction of the granular material and slightly underestimated by the model.

Discussion and comparison
The subaqueous granular collapse clearly exceeds the capabilities of the two-component model. The high sensitivity to the initial packing density can not be explained with this model and the loose and dense simulation are virtually similar. Results of the two-component model lie between the two extreme cases of the loose and dense experiment, matching the velocity of the loose and the run-out of the dense experiment. This is reasonable, considering that the missing excess pore pressure is stabilizing the dense column and destabilizing the loose column. This model is not sufficient for a practical application, as the runout is substantially underestimated in the loose case. Extremely long run-outs on slopes with 2 • inclination have been observed in nature (e.g. Bryn et al. 2005) and they can not be explained with a granular two-component model.
The two-phase model can take advantage of its ability to capture excess pore pressure. It  (Figs 12f and  13e) and a consistent sensitivity to pore pressure and initial packing density (Fig. 14). The failure mechanism of both, the dense and the loose experiment, are successfully simulated (see Fig. 15), indicating that the two-phase model captures the most important physical phenomena. The model fails in two aspects, as the pore pressure peak in the loose case and the time scale in the dense case differ by a factor of 2 and 3, respectively. It should be noted that no exhausting parameter fitting was required for these results. Solely the critical state line is optimized to yield the correct pore pressure, all other parameters were selected a priori following Rondon et al. (2011), Savage et al. (2014), and Boyer et al. (2011. Notably, some of the issues, e.g. the overestimated velocity of the loose collapse, might be resolvable with fitting parameters. Furthermore, the model allows us to simulate both cases with the same set of parameters with good accuracy. This distinguishes this work from earlier attempts (e.g. Savage et al. 2014;Wang et al. 2017;Si et al. 2018a), where some parameters were fitted individually to the dense and loose case.
Excess pore pressure plays an important role in subaquatic experiments because it controls shear strength and friction. Dilatancy, compaction, and the dynamic particle pressure further influence friction and thus the kinematics of the slide. The dense column is only able to collapse after decreasing its packing density and thus its effective shear strength. The column is dilating until reaching the limiting packing density m . Before this packing density is reached, the shear rate is limited to the creeping shear rate 0 . A relatively high value was used for this parameter and a lower creeping limit would be desirable, especially considering the error of the time scale in the dense simulation (see appendix B.5). However, strong oscillations were observed when choosing lower values for 0 because the shear rate often exceeded 0 before dilating sufficiently. The bottom of the dense column is further compacting in the simulation, up to a packing density of 0.604. This is reasonable as the initial particle pressure of 303.3 Pa at = 0.60 is below the overburden pressure of 370.8 Pa of the pile. At the same time, negative excess pore pressure can be observed at the bottom of the column. Compaction and negative excess pore pressure seem to contradict each other at first glance. However, the negative excess pore pressure in the upper parts of the column is so strong, that fluid is flowing upwards from the bottom of the column. This can be seen in the relative velocity field (Fig. 12h), but also the gradient of pore pressure (Fig. 12b) indicates that pore liquid will flow upwards.
The front speed of the loose collapse is entirely controlled by the dynamic contribution of the ( )-( )-rheology to effective pressure and friction. Simulations with critical state theory (constant friction coefficient and the quasi-static effective pressure model of Johnson & Jackson (1987)) exceed the experimental runout by far (see Appendix B.4). This is a strong contrast to the subaerial case where acceptable results could be achieved with critical state theory.
The dynamic contribution to particle pressure and friction plays also an important role in the dense case, although this pile collapses very slow. The thin layers of grains that are breaching from the unsupported column flank reach packing densities far below rlp = 0.53 due to mixing with the ambient fluid. At this packing density, the quasi-static contribution to effective pressure vanishes, and the runout of these particles is entirely controlled by dynamic particle pressure and friction. The runout of the breaching flank could not be controlled in simulations with critical state theory (see appendix B.4). The pore pressure in the loose case differs qualitatively and quantitatively from the measurement. Within the applied model, it seems reasonable that a high initial peak decreases quickly, as substantial amounts of grains and thus overburden pressure leave the region of the pressure sensor. Similar results with an early, short, and strong peak and a slow further dissipation, close to the measurement, have been obtained with other frameworks, e.g. by Bouchut et al. (2017) or Baumgarten & Kamrin (2019).
The dilatancy of the dense column is substantially faster in the numerical model than in the experiment, although the permeability is underestimated following the comparison of the pore pressure. Therefore it is unlikely that permeability is the cause for this discrepancy and we assume that inaccuracies in the rheology are responsible. The ( )-( )-rheology describes the steady shearing of a fluid-grain mixture very well (Boyer et al. 2011). However, the transient transition towards the steady packing density at a certain effective pressure is not described. This transition depends on the permeability of the granular material but also on its viscosity (shear and bulk viscosity). As mentioned before, the high value for the creeping shear rate 0 could be responsible for this issue but it might also be related to the missing bulk viscosity or a mismatch of constitutive parameters. Bulk viscosity could delay the dilatancy in the early stage of the dense collapse, bringing the time scale of the collapse in the simulation closer to the experiment. Bulk viscosity could further help to decrease the pore pressure peak in the loose case, as some of the pore pressure could be transformed into viscous pressure. Schaeffer et al. (2019) suggests a form for the bulk viscosity which has the potential to improve these aspects. Savage et al. (2014) and Si et al. (2018a) include a cohesive shear strength into their model to correct some of these problems and to fit results to the experiment. However, there is no evidence for cohesive forces in a fully submerged granular flow. Neither electrostatic forces nor cementing have been reported by Rondon et al. (2011). Apparent cohesion can be traced back to negative excess pore pressure, which is directly simulated by the numerical model. Notably, Si et al. (2018a) are able to control the slide velocity very well. However, this is achieved by fitting the cohesion to the respective case and by a strong overestimation of the negative excess pore pressure, reaching values around 500 Pa at the pressure sensor at = 3 s (see Fig. 5 by Si et al. 2018a)).
Baumgarten & Kamrin (2019) applied a similar model (elasto-plastic multiphase model with ( )-( )-scaling) to the same cases. The results show similar problems, i.e. an overestimation of the pore pressure in the loose case and an overestimation of the collapse velocity in the dense case. Notably, we achieve similar results in these test cases with a substantially simpler model.

Conclusions
The Navier-Stokes Equations can be an adequate tool for accurate simulations of granular flows when they are complemented with the correct rheologies. Substantial progress has been made in recent years with the ( )-rheology and its extensions to compressible flows and low Stokes number flows. The incompressible ( )-rheology fits well into the multi-component framework of OpenFOAM and the compressible ( )-( )-rheology fits well into the multiphase framework, as previously shown by e.g. Chauchat et al. (2017). We apply, for the first time, the compressible ( )-( )-rheology to granular collapses and avalanching flows. The superposition with the critical state theory is imperative to get realistic packing densities at rest and a stable solver. For subaerial, i.e. high Stokes number flows, dilatancy plays a minor role and results of the compressible model are similar to the incompressible model. However, the dilatancy predicted by the compressible model is able to close the gap between the experiments and the incompressible model. Further, the compressible model should be well-posed (Barker et al. 2017;Heyman et al. 2017;Schaeffer et al. 2019), in contrast to many incompressible granular flow models (Barker et al. 2015). Note that bulk viscosity, which is imperative for a well-posed rheology (e.g. Schaeffer et al. 2019), was not considered in this study. However, the coupling of the granular phase to the pore fluid has a similar effect as bulk viscosity and might be able to restore a well-posed system. For a guaranteed wellposed compressible rheology that collapses to the ( )-( )-rheology in steady state, the reader is referred to Schaeffer et al. (2019).
The upsides of the compressible two-phase model come at the cost of more parameters and a stronger mesh dependence. Furthermore, code and case setup are more complicated with the two-phase model and simulations are more prone to failure if initial conditions or parameters are not well suited for the case. Therefore, the incompressible model might be better suited for some flows at high Stokes numbers, especially considering regularized rheologies that are well-posed for a wide range of flow regimes (e.g Barker & Gray 2017). Notably, we did not encounter any problems with the partial ill-posedness of the ( )-rheology, which could be related to relatively coarse grids, high numerical diffusion, the short simulation duration or the truncation of the viscosity.
The extension to low Stokes number flows is made possible by the ( )-( )-rheology. At low Stokes numbers, it is imperative to consider excess pore pressure and a two-phase model is required. Therefore, the incompressible ( )-rheology is rather impractical and becomes only applicable after supplementing it with the ( )-curve to the compressible ( )-( )rheology. The dynamic growth of pressure and friction is substantial for accurate results, highlighting the value of the ( )-( )-rheology. The fitting of parameters was reduced to a minimum and only the critical state line had to be optimized to the experiments. It should be noted that these parameters could be determined by measuring the critical packing density at a few pressure levels, making the simulations free of any fitted parameter. The compressible two-phase model reacts sensitive to the packing density, recreating the final runout, pile shape, and failure mechanism of the experiments very well. The model still lacks in some aspects, e.g. the time scale and the velocity of the dense collapse and the pore pressure peak in the loose collapse.
It was shown that the incompressible two-component model can be derived from the compressible two-phase model by neglecting the relative velocity between phases. This simplification yields reasonable results for subaerial granular flows at high Stokes numbers but fails to describe the subaquatic granular flows at low Stokes numbers. This seems to be contradictory, as the relative velocity (which was neglected in the incompressible model) is very small in the subaquatic case (see Figs. 12 and 13) but considerable high in the subaerial case (see Fig. 9). This apparent paradox can be resolved by the fact that unhindered density changes have no notable influence on the flow dynamics. However, if changes in packing density are constrained, pore pressure will build up and the rheology of the material will change drastically. Thus, pore pressure, rather than compressibility is the key factor that allows the two-phase model to accurately capture the flow mechanics. The two-phase model provides many other upsides aside from the inclusion of pore pressure. The continuous transition from dense granular material to pure ambient fluid should be useful for the simulation of granular free streams (Viroulet et al. 2017), turbidity currents (Heerema et al. 2020) and powder snow avalanches (e.g. Sovilla et al. 2015). Other studies showed that the two-phase model is useful for sediment transport (Chauchat et al. 2017) and other dilute particle-fluid mixtures (e.g. Passalacqua & Fox 2011).
OpenFOAM provides a good platform to evaluate concepts (e.g. the multi-component and multi-phase methodology) and models (e.g. ( )-( ) and ( )-( )-rheologies). The implemented rheologies can be further coupled with segregation (Barker et al. 2021) or tsunami simulations (e.g. Si et al. 2018b). However, the segregated semi-implicit solver strategy of OpenFOAM sets limits to models and execution velocity, as (part of the) viscous terms and the particle pressure are included explicitly. This showed to be problematic and a fully implicit solver, that solves all equations simultaneously, might be superior in this regard.
The model can help to understand the extreme outruns of submarine landslides, such as the Storegga landslide (e.g. Bryn et al. 2005) and the big variation in tsunamigenic potentials (e.g. Løvholt et al. 2017). Theories, such as hydroplaning and remoulding (e.g. De Blasio et al. 2004) can be quantitatively described by the critical state theory and its dynamic extension in form of the ( )-( )-rheology. Hydroplaning, formerly described as the flowing of sediment on a thin layer of liquid can be interpreted as a region of low or even zero packing density and vanishing effective pressure. This can be observed in Fig. 15g-i, where the front of the loose slide is lifted by pressure in the surrounding fluid. Remoulding can similarly be explained with critical state theory as an overconsolidated sample that is dilating during shearing (see Fig. 15a-f). The two-phase model and its capability to describe various and realistic failure mechanisms with different time scales are particularly valuable for understanding the tsunamigenic potential of submarine landslides and the respective slopes. The dense column collapses very slowly, reaching velocities of up to 0.1 m s −1 in small layers near the surface. The loose column collapses entirely with velocities up to 0.4 m s −1 . The tsunamigenic potential of a landslide scales with initial acceleration and the mobilized volume (e.g. Løvholt et al. 2017) and a substantial difference in tsunamigenic potential follows for the dense and the loose slide. This shows that packing density, excess pore pressure and permeability are key parameters in controlling stability, failure mechanism, slide acceleration, and tsunamigenic potential.
Many full-scale subaquatic landslide simulations are based on Bingham fluids, a visco-plastic rheology independent of the pressure (e.g. Kim et al. 2019). This seems to stand in strong contradiction to the model applied here. However, the simulation of the loose case shows that packing density changes are small. For a nearly constant packing density, the effective pressure decouples from overburden pressure because the weight is absorbed entirely by pore pressure. As a consequence, overburden pressure and friction will decouple and the microscopic granular friction will appear as cohesion on a macroscopic scale. The macroscopic description as a Bingham fluid is therefore surprisingly consistent with the findings in this work, especially for fine grained marine sediments with low permeabilities.

Summary
This work highlights a path to extend the incompressible ( )-rheology for subaerial granular flows to the compressible ( )-( )-rheology for subaquatic granular flows. The implementation of the ( )-( )-rheology in a multiphase framework and the ( )rheology in a multi-component framework allows us to conduct subaerial granular collapses with two different models. The application shows consistency between the incompressible ( )-rheology (e.g. Lagrée et al. 2011) and the compressible ( )-( )-rheology. Notably, substantial modifications to the ( )-curve are required for a practical application of the rheology. The simulations show that compressibility and dilatancy have a small influence on high Stokes number flows because excess pore pressure is negligibly small.
The implementation of the ( )-( )-rheology extends possible applications to low Stokes number flows, e.g. subaquatic granular collapses. The incompressible model reaches its limitations under these conditions and the compressible model is required for an accurate simulation. Other than previous attempts, we applied the exact same set of parameters to an initially dense and loose granular collapse with satisfying results. Notably, the application of the ( )-( )-rheology does not require an extensive fitting of constitutive parameters. The comparison between the compressible model and experiments uncovered discrepancies in the time scale and the pore pressure. These could be indicators for issues in the rheology, e.g. a missing bulk viscosity or issues with the creeping regime that had to be introduced for numerical stability. The well-posedness of the proposed model is not guaranteed and should be investigated in the future.
The compressible two-phase model has a wide range of applications and the results have implications on many problems in geoscience. Applications to sediment transport and scouring (Cheng et al. 2017) have been shown with a similar model. We further expect the applicability to turbidity currents and all other gravitational mass flows with low and high Stokes numbers. Furthermore, Si et al. (2018b) showed the applicability of a similar model to landslide tsunami simulations by incorporating the free water surface. two-component, subaerial for a solid-like behaviour can be calculated as as a function of respective scales by choosing the magnitude of viscous stresses over other terms, . The required magnitude for can be estimated by conducting a numerical sensitivity analysis.
Variations of max (and thus ) are presented in Fig. 16 for the subaerial case at = 0.8 s, using the two-component model. The value of max = 1 m 2 s −1 is adequate for this example and the left side of the pile stays nearly static as it is the case in the experiment. The respective value for the dimensionless scaling factor follows as 0.1 ( = 0.1 m, | | ≈ 10 m s 2 ), indicating that viscous forces have to be about 10 times higher than all other contributions to the momentum conservation equation. For lower viscosities, the pile is notably deformed and shows no stable regions and no granular characteristics. Rather, the pile shows the characteristics of a visco-plastic fluid (see rheology comparisons by Lagrée et al. 2011), which indicates that the viscosity threshold was dominating the simulation. Notably, cases with high maximum viscosity are stable after = 0.4 s while cases with low maximum viscosity keep flowing beyond = 0.8 s. For an application of granular rheologies in OpenFOAM, we suggest a maximum viscosity following Eq. 2.28 with = 1/10.

B.2. Grid sensitivity
The grid sensitivity is an important issue for complex flow models and we provide a full grid sensitivity analysis for the multi-component and the multi-phase model. The grid sensitivity study for the multi-component model is solely conducted for the subaerial case because the mechanics of this model is similar in all cases. The final pile shape of the investigated case is shown in Fig. 17 for various grid resolutions. This model reacts very robustly to coarse grids and 30 cells along the pile height are sufficient to get accurate results for the final pile shape.
The two-phase model is more complex in terms of grid sensitivity. Three different failure mechanisms of the granular column can be observed in the simulations with the two-phase model. The three mechanisms react differently to a variation of the grid resolution. In some cases, the model reacts very sensitively to coarse grids and a grid refinement study should always be performed when applying this model. Fig. 18 shows the grid convergence analysis for the subaerial, the subaquatic dense, and subaquatic loose case. The two-phase model is slightly more sensitive to coarse girds than the two-component model in the subaerial case, two-component, subaerial see Fig. 18a. The problematic area is the thin flow front and the issue is probably related to the mixed role of the phase-fraction field.
The two-phase model is very sensitive to coarse grids in the subaquatic dense case, see Fig. 18b. Breaching of a thin layer of grains on the unsupported side of the column leads to a reduced phase-fraction in cells that contain the slide surface. This reduces the effective pressure in all of those cells and further the shear strength, accelerating the collapse. The result is a mesh dependency of the final pile shape and the collapse velocity. The mesh had to be refined down to a cell size of 0.0005 m, to achieve accurate results. Notably, the difference between the smallest two cell sizes is still remarkably at the front of the collapse. An additional refinement step would be desirable, but this would have exhausted the available computational resources.
The loose subaquatic case can be simulated with good accuracy on a relatively coarse mesh as shown in Fig. 18c. This is not surprising as the failure mechanism and flow pattern is much simpler. In particular, the effective pressure discontinuity at the free surface is weaker than in the dense case and thus requires a smaller grid resolution.

B.3. Time step duration sensitivity
The time step duration is investigated similarly as the grid resolution. The time step duration is not fixed but adapted to velocity (CFL conv ) and viscosity (CFL diff ), relative to the grid size. The viscous contribution to the CFL number is always much bigger in the investigated cases and the time step sensitivity study was conducted based on this value. Notably, stability can be guaranteed only for CFL diff < 1. However, the CFL diff number allows only an estimation of the stability and in some cases larger time steps are possible.
The two-component solver can operate in two modes, using a full momentum predictor step or a reduced momentum predictor step. The full momentum predictor step solves the full linearised system of the discretized momentum conservation equation. The reduced momentum predictor step calculates an explicit prediction of the velocity field based on the velocity field of the last time step. This has a substantial influence on the stability when viscous stresses are dominating. Fig. 19a shows the final pile shape for various time step durations and the full momentum predictor step in the two-component model. Notably oscillations in pressure can already be seen at CFL diff = 2 (not shown) and they grow substantially for higher CFL diff numbers. The pressure oscillations start to influence the pile shape at CFL diff = 10 and the pile is completely distorted for CFL diff = 100.
The two-component model is more robust to larger time steps when operating with the reduced momentum predictor step, see Fig. 19b. Pressure oscillations start approximately at two-phase, suberial (a) two-phase, subaquatic, dense (b) two-phase, subaquatic, loose (c) Figure 18: Grid sensitivity of the two-phase model for the subaerial case (a), the dense subaquatic case (b) and the subaquatic loose case (c). The black dashed lines show the experimental final pile shapes for reference. CFL diff = 100 and the first influence on the slide geometry can be observed at CFL diff = 1000. Anyway, it is recommended to run the two-component model with small enough time steps to prevent pressure oscillations, ideally at CFL diff = 1, as done in this work.
The two-phase model reacts less sensitive to large time step durations. In fact, simulations were stable up to CFL diff = 1000. No pressure oscillations could be observed and the final pile shape is nearly unaffected for all cases, see Fig. 20. However, the accuracy got worse for large time step durations and we observed a slower initial acceleration for CFL diff = 1000. No subaquatic simulations with the two-phase model and CFL diff = 1 have been conducted, as this would have exhausted the computational resources.

B.4. Influence of dynamic friction and wall friction
The effect of the dynamic increase of friction following the ( )-rheology and the effect of the reduced wall friction were investigated with both models. Results are shown in Fig. 21a for the two-component model and in Fig. 21b for the two-phase model. The models do not react sensitively to the variation of the friction model and the basal friction coefficient. The runout is slightly underestimated in simulations with the ( )-rheology on a rough surface. However, the introduction of the smooth surface elongates the runout and the simulations fit the experiment well. Simulations with a constant friction coefficient on a rough surface fit the experiment also well, simulations with constant friction coefficient on a smooth surface overestimate the runout. The dynamic contribution to effective pressure and friction is imperative for simulations at low Stokes numbers. Fig. 22 shows the dense and loose subaquatic two-phase simulations with critical state theory and ( )-( )-theory. The outrun cannot be controlled without the dynamic contributions of the ( )-( )-theory and exceeds the final runout very quickly.

B.5. Influence of the creep shear rate
The influence of the creep shear rate 0 is shown in Fig. 23 for the subaquatic dense case. The Figure shows an early time step at = 1.0 s and the final pile shape at = 10 s. A reduction of the creeping shear rate from 5 s −1 to 1 s −1 leads to a slower initial collapse of the column. The delayed collapse is desirable and brings the simulation closer to the experiment. However, the low value for 0 leads to oscillations in the particle pressure because the simulation accelerates and reaches shear rates beyond 0 too quickly. The simulation with 0 = 5 s −1 shows no oscillations and has thus been utilized in the main part of this work. Notably, a smaller time step duration will allow smaller values for 0 , however for an increased computational cost. The final pile shape is barley affected by the change in 0 .  two-phase, subaquatic, dense (a) two-phase, subaquatic, loose (b) Figure 22: The dense granular collapse at = 6.0 s (a) and the loose granular collapse at = 0.65 s (b), simulated with critical state theory and ( )-( )-rheology. The dashed black line shows the final experimental pile shape. The simulations with critical state theory clearly exceed the experiment early in the simulation.