the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
OpenFOAMavalanche 2312: depthintegrated models beyond denseflow avalanches
Matthias Rauter
Julia Kowalski
Numerical simulations have become an important tool for the estimation and mitigation of gravitational mass flows, such as avalanches, landslides, pyroclastic flows, and turbidity currents. Depth integration stands as a pivotal concept in rendering numerical models applicable to realworld scenarios, as it provides the required efficiency and a streamlined workflow for geographic information systems. In recent years, a large number of flow models were developed following the idea of depth integration, thereby enlarging the applicability and reliability of this family of process models substantially. It has been previously shown that the finite area method of OpenFOAM^{®} can be utilized to express and solve the basic depthintegrated models representing incompressible dense flows. In this article, previous work (Rauter et al., 2018) is extended beyond the denseflow regime to account for suspended particle flows, such as turbidity currents and powder snow avalanches. A novel coupling mechanism is introduced to enhance the simulation capabilities for mixedsnow avalanches. Further, we will give an updated description of the revised computational framework, its integration into OpenFOAM, and interfaces to geographic information systems. This work aims to provide practitioners and scientists with an opensource tool that facilitates transparency and reproducibility and that can be easily applied to realworld scenarios. The tool can be used as a baseline for further developments and in particular allows for modular integration of customized process models.
 Article
(7846 KB)  Fulltext XML

Supplement
(2746 KB)  BibTeX
 EndNote
Runout and impact simulations of gravitational mass flows typically rely on depthintegrated models (e.g. Pitman et al., 2003; Sampl and Zwinger, 2004; Christen et al., 2010; Iverson and George, 2014; Mergili et al., 2017; Eglit et al., 2020). In comparison with fully resolved threedimensional models, this framework provides a range of upsides: the computational expense is substantially reduced, interface and phase tracking are simpler and more reliable, and integration in geographic information systems is straightforward. Depthaveraged models are easier to solve numerically, to set up, to calibrate, and to evaluate. However, depth integration comes at a price: the vertical flow structure (including the profiles of density, velocity, and shear rate) is lost, and all related effects, if needed for closures, have to be reintroduced with additional models. This includes friction, erosion of basal material, and its deposition (e.g. Rauter and Köhler, 2020), as well as layering of varying regimes (e.g. Bartelt et al., 2016). A possibility of overcoming this entails the shallow moment approach (Kowalski and Torrilhon, 2019); however, it has not been applied successfully to realworld granular mass flows yet. Nevertheless, depthintegrated models have proven to be a good compromise between simplicity and complexity, especially for flows of geographic extent from avalanches (Christen et al., 2010) to tsunamis (Løvholt et al., 2015).
Granular flows show a large variety of behaviours. A very strong distinction of properties can be linked to the Stokes number St, expressing the ratio between inertia and drag forces on particles (Boyer et al., 2011; Rauter, 2021). For a flow with shear rate $\dot{\mathit{\gamma}}$ of granules with density ρ_{g} and diameter d, in a medium of viscosity ν_{c}, and density ρ_{c}, the Stokes number can be written as
At high Stokes numbers, drag forces are small and particles move freely through the surrounding fluid or gas. Thus, the bulk motion is dominated by particle–particle interactions, and particles will arrange in a relatively high packing density that only depends on the local shear rate and pressure (e.g. Forterre and Pouliquen, 2008). Furthermore, for many realistic problems, the bulk density can be assumed constant with acceptable accuracy. Denseflow models often take advantage of this fact and are formulated as incompressible nonNewtonian fluids (e.g. Savage and Hutter, 1989; Rauter, 2021).
At low Stokes numbers, drag on particles is substantial and particles are not able to rearrange freely within the carrier medium. Particles and surrounding fluid form a suspension and move like a single fluid, only to be slowly separated by the settling velocity. The packing density or volume fraction depends on various aspects and most importantly on the history of the flow. This is a strong hint that the volume fraction requires an evolution equation to be properly described (as done by, for example, Parker et al., 1986; Kowalski and McElwaine, 2013; Issler et al., 2018; Rauter, 2021).
It can be seen from Eq. (1) that the Stokes number depends on the particle size. In polydisperse granular flows, i.e. flows with particles of various sizes (e.g. Barker et al., 2021), this can lead to vertical segregation of small and large particles and thus a coexistence of both regimes. This can be well observed in snow avalanches (Sovilla et al., 2015), where a dense flow is formed by relatively coarse snow blocks of size 10^{−1} m (Bartelt and McArdell, 2009; Rauter et al., 2018), and a powder cloud is formed by small ice particles of size 10^{−4} m (Rastello et al., 2011; Bartelt et al., 2016); see Fig. 1.
In terms of depthintegrated models, this calls for a twolayer model, capturing the dense flow with an incompressible model and the powder cloud with a suspension model (Issler, 1998; Sampl and Zwinger, 2004; Bartelt et al., 2016).
In this work, we will extend the denseflow model of Rauter et al. (2018) to low Stokes number suspension flows following the model of Parker et al. (1986). We will make and evaluate some adjustments to account for highdensity differences between the carrier medium and the particles. In a further step, we will combine the models for dense flow and suspension into a twolayer model, capable of simulating mixedsnow avalanches, similar to Turnbull and Bartelt (2003) and Bartelt et al. (2016). For this purpose, we have to define a coupling mechanism, i.e. a mass flux term that feeds the powder cloud from the dense core. We develop a novel idealized relation that encapsulates the essential features of this process and deliberately avoids more complex mechanisms (e.g. Sampl and Zwinger, 2004; Bartelt et al., 2016). We focus on clarity, simplicity, and modularity and therefore describe all processes with simple, local relations that can be formulated independently of one another. This is motivated not only by the goal of creating a simple baseline model but also by the observation that complexity does not necessarily lead to better results (Zhao and Kowalski, 2022). The natural terrain is handled as described previously by Rauter et al. (2018). While the main focus of the presented work is snow avalanches, the implementation might very well be useful for the simulation of turbidity currents, as several researchers suspect a dense core in these flows as well (e.g. Heerema et al., 2020).
The naming convention of layers and fluxes follows Bartelt et al. (2016): the dense core is denoted by Φ, the suspension flow by Π, the static bottom layer by Σ, and the stationary ambient fluid by Λ. Flow fields are marked by the respective subscripts and fluxes between layers with two subscripts and an arrow indicating the direction of the flux (see Fig. 1).
The numerical solution and implementation are based on the finite area method (Tuković and Jasak, 2012; Rauter and Tuković, 2018) as implemented in OpenFOAM. Its modular structure and building blocks have proven to be flexible and highly valuable for physical depthintegrated models. Various code parts are reused between all models and various communities, in particular the numerical solver, geometry, and data handling, as well as various code parts related to the physics of the flow, such as friction models. Besides the introduction of the new model and its capabilities, this work highlights the extendability of the basic OpenFOAM solver to complex models.
The toolchain to process the basic terrain data, all the way to the final simulation visualization, has been improved substantially since the work of Rauter et al. (2018), and many external dependencies have been removed in order to facilitate a tight integration into OpenFOAM. Consequently, this paper also provides an updated overview of the toolchain and its practical applications. In this context we will also give a revised introduction into the finite area method and the specific derivations of depth integration. The model caters to practitioners who need a simple mixedsnow avalanche model but mostly to scientists who wish for an open model and framework that can be easily modified and extended to evaluate new concepts and ideas.
The novel model is evaluated with various synthetic test cases and finally applied to two real events: the 1988 Wolfsgruben avalanche and the 2019 Eiskar avalanche.
2.1 Conservation equations and depth integration
The presented method fundamentally relies on balance equations, in particular, the conservation of mass and momentum for fluids. The combination of these two equations is widely known as the Navier–Stokes equations (e.g. Ferziger and Peric, 2002) and can be written as
with the bulk density ρ and the bulk velocity u. (Note that it can also be defined for an individual phase with some modifications, see, for example, Rauter, 2021.) These flow fields are functions of time t and space $\mathit{x}=(x,y,z{)}^{T}$. The model Eqs. (2) and (3) describe their evolution from a known initial state $\mathit{u}\left(\mathrm{0},\mathit{x}\right)={\mathit{u}}_{\mathrm{0}}\left(\mathit{x}\right)$ and similar ρ_{0}(x), under the influence of boundary conditions. The divergence of the stress tensor T has the effect of diffusing momentum, and the volume force f represents additional forces, such as gravity.
Appropriate closure relations that express the stress tensor T as a function of the unknown flow fields yield a wellposed problem that can, in principle, be solved with numerical methods (Barker and Gray, 2017). However, even a wellposed problem is often not practically solvable from a computational perspective. Therefore, multiple simplifications have to be made to make problems of practical relevance accessible. Simplifications often come in the form of averaging over a certain time or over space to get rid of turbulent structures (Reynolds averaging; see, for example, Ferziger and Peric, 2002), to describe the average behaviour of multiple interpenetrating phases (phase averaging; e.g. Rauter, 2021), or to get rid of the vertical dimension (e.g. Savage and Hutter, 1989; Rauter and Tuković, 2018). The latter is referred to as depth averaging or depth integration and avoids the calculation of threedimensional flow details. It yields mean values of, for example, density $\stackrel{\mathrm{\u203e}}{\mathit{\rho}}$ and velocity $\stackrel{\mathrm{\u203e}}{\mathit{u}}$ along the depth.
In the simplest case, where the depth integration is aligned with a spatial axis (e.g. the z axis), the problem can be reduced from three $(x,y,z)$ to two dimensions (x,y). In this case, the depthaveraged value for an arbitrary field ψ is defined as
The newly introduced field $h(x,y,t)$ describes the flow depth, here in terms of the z coordinate of the top boundary of the integration, for a bottom boundary assumed to be aligned with z=0. The bottom and top boundaries are usually defined such that the mass flux through them is zero, meaning that they move with the vertical velocity of the flow at the respective position. The simplest example of such a model is the shallowwater equations (Barré de SaintVenant, 1871). Defining the boundary in any other way will lead to additional source or sink terms, depending on the mass flux through the boundary (e.g. Pudasaini and Hutter, 2007). Examples would be any kind of entrainment and deposition fluxes.
Depthintegrated models are often considered synonymous with twodimensional models. However, real avalanches and landslides travel along paths and surfaces in threedimensional space. The threedimensional nature of the terrain has to be reintroduced by modifying the twodimensional model equations. Most often this is accomplished by abandoning Cartesian coordinate systems and Euclidean geometry, which was described in detail first by Savage and Hutter (1989, 1991) and extended by many others since then (e.g. Bouchut and Westdickenberg, 2004; Denlinger and Iverson, 2004; Pudasaini et al., 2005; Hergarten and Robl, 2015). This introduces various correction terms based on Christoffel formalism that are difficult to handle in complex models. In practice, simpler approximations are frequently employed (e.g. in RAMMS, see Fischer et al., 2012), leading to a disparity between theory and practical implementation. Notably, many of these developments happened in parallel and independently in the Russian avalanche dynamics community (Eglit et al., 2020).
An alternative to twodimensional models with excessive curvature terms is the direct solution of the governing equations in threedimensional space (Craster and Matar, 2009; Hagemeier et al., 2011; Rauter and Tuković, 2018). Depth integration is still compatible with this approach, and it can in principle be conducted in any direction pointing out of the surface. Yet in this work, depth integration is always conducted in the direction of the normal vector n^{Γ} to the flow surface Γ, as shown in Fig. 2. This has formally to be conducted in a surfacealigned coordinate system x^{′}–y^{′}–z^{′}:
The Jacobian matrix J, representing the transformation $\partial {\mathit{x}}^{\prime}/\partial \mathit{x}$, and its determinant det(J) take into account the curvature of the surface and its influence on the volume in a differential volume element of the flow (Bouchut et al., 2003). This effect is of the order of $h/R$ (Bouchut et al., 2003) with the mean curvature radius R and thus is small for mildly curved surfaces (R is large in comparison to the flow height h). As in most other models, the influence of the curvature on depth integration is ignored in this work.
2.2 Surface partial differential equations
Depth integration in terms of Eq. (5) projects all threedimensional flow fields on the surface Γ they are constrained by. The conservation equations can then be expressed as surface partial differential equations (SPDEs) that are defined on the surface Γ and include derivatives of various fields along it (Rauter and Tuković, 2018). These derivatives emerge from depthintegrating the nabla operator ∇ present in the Navier–Stokes equations formulated in a threedimensional Cartesian reference frame. The presented framework differs from the classic approach of handling derivatives in two respects. Both are described in the following.
Depth integration has different effects on derivatives taken along the surface compared to those taken normal to the surface. While derivatives normal to the surface either vanish or appear as local source terms, derivatives along the surface remain in the system. A very common approach is to write the surfacealigned derivative as a twodimensional derivative in local coordinates (e.g. in terms of $\partial {x}^{\prime}$, $\partial {y}^{\prime}$; Savage and Hutter, 1989). In the present framework, however, we express all entities in global Cartesian coordinates. For diffusive processes, this procedure gives rise to the Laplace–Beltrami operator (Dziuk and Elliott, 2013). However, this technique can also be adopted to other differential operators present in the system. The respective surface gradient and divergence operators ∇^{Γ} can be readily calculated in the numerical framework; see Sect. 2.3.
The depth integration of derivatives is then conducted in analogy to ordinary fields (see Eq. 5), in the surface aligned coordinate system x^{′}–y^{′}–z^{′}:
where x_{b} is a point on the bottom of the flow (and thus the flow surface Γ) and x_{t} the corresponding point on the free surface of the flow. The second term on the righthand side of Eq. (6) represents an additional sink or source term, which arises if ψ is not zero at the bottom, x_{b}, or the top of the flow, x_{t}, such as due to entrainment or basal friction.
Further, our approach does not follow, for example, Savage and Hutter (1989) in separating the z component (determining, for example the pressure) from the x and y component (determining, for example, the velocity) in vectorial type balance laws such as the momentum equation. This would not be possible as our coordinate system is not aligned with the surface. We rather project the full threedimensional equation onto the surface and the normal vector. The surface tangential projection of the surface gradient of a scalar is hence given by
and the surface normal projection as
and similar for vectors and divergence operators. The benefit of this approach is that the framework operates entirely in global Cartesian coordinates; thus, similar to an inertial frame, no fictitious centrifugal forces have to be considered. It follows that leadingorder curvature effects are considered in the model by design, without explicitly expressing the curvature (Rauter and Tuković, 2018). For higherorder curvature effects, det(J) needs to be preserved during depth integration (Bouchut and Westdickenberg, 2004).
With these building blocks and some knowledge on how to transform onedimensional shallowflow models (e.g. Savage and Hutter, 1989; Parker et al., 1986), it is possible to extend nearly arbitrary depthintegrated flow models to complex terrain. In particular, ordinary depthintegrated flow models represent the surface–tangential momentum conservation equation and the flow depth equation. The twodimensional ∇ operators have to be replaced with the surface–tangential ${\mathbf{\nabla}}_{\mathrm{s}}^{\mathrm{\Gamma}}$ operators. The surface–normal momentum conservation equation can be applied to replace the usually simplified expression for the basal pressure.
2.3 Finite area method
Partial differential equations, as well as their SPDE counterparts, are rarely solvable in an analytical sense, especially practical problems that represent realworld situations. Therefore, we rely on numerical approximations of SPDEs and the finite area method. This method is a variation of the finite volume method (see Ferziger and Peric, 2002; Jasak, 1996; Moukalled et al., 2016, for details) in N+1 dimensions, where N is the dimension of the control volumes. This means that for twodimensional control volumes (i.e. surfaces), vectorial entities, such as normal vectors, velocities, and fluxes, will be threedimensional. Similar to the conventional finite volume method, the Gaussian surface theorem (Tuković and Jasak, 2012) is applied and discretized by simplifying a control surface S as a flat, convex polygon S_{i}, as shown in Fig. 3. The expressions for the differential operators follow as
and
Index e refers to a discrete number of straight edges that form the polygon with surface S_{i}. ψ_{e} is the average value of the field ψ on the edge e, L_{e} its length, and m_{e} is the Γ–tangential and edgenormal outwardpointing vector. S_{i}, L_{e}, and m_{e} are purely geometrical properties that are defined during mesh generation. Values of fields on edges ψ_{e}, on the other hand, are interpolated from values of edgeadjacent cells, ψ_{P} and ψ_{N}. This introduces flux transport across cells and represents the flow of mass or information from one cell to neighbouring ones. The fluxes can then be associated in a linear system of equations that is solved with a suitable method.
Discretization of nongradient terms (e.g. the temporal derivative or any source term) is done in complete analogy to the finite volume method and obtained from integration over the control surface S_{i}. For details we refer to the large amount of excellent literature on the finite volume method (Ferziger and Peric, 2002; Jasak, 1996; Moukalled et al., 2016; LeVeque, 2002).
The denseflow model describes the flow of incompressible material with density ρ_{Φ} (see Fig. 1). In the case of a granular mass flow, the density follows from the grain density ρ_{g} and the volumetric packing density ϕ_{Φ} as
However, fluids can be simulated with this model as well, in which case ρ_{Φ} is the intrinsic density of the fluid. The depthintegrated mass and momentum conservation equations follow as
The unknown flow fields are the flow depth h_{Φ}, the depthintegrated velocity ${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Phi}}$, and the basal pressure p_{Φ}. The gravitational acceleration is represented by its surface–tangential projection ${\mathit{g}}_{\mathrm{s}}=\left(\mathbf{I}{\mathit{n}}^{\mathrm{\Gamma}}\phantom{\rule{0.125em}{0ex}}{\mathit{n}}^{\mathrm{\Gamma}}\right)\mathit{g}$ and its surfacenormal projection g_{n}=(n^{Γ} n^{Γ})g. Equation (14) represents the surfacenormal component of the momentum conservation equation and yields the basal pressure p_{Φ}.
The shape factor ξ_{Φ} partially compensates for errors introduced by switching integration and multiplication, namely ${\mathit{\xi}}_{\mathrm{\Phi}}\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Phi}}\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Phi}}=\stackrel{\mathrm{\u203e}}{{\mathit{u}}_{\mathrm{\Phi}}\phantom{\rule{0.125em}{0ex}}{\mathit{u}}_{\mathrm{\Phi}}}$. It depends on the velocity profile and as such on the constitutive model and the state of the flow. It is usually neglected or set to a theoretical and constant value, derived, for example, from the Bagnold (1954) velocity profile (${\mathit{\xi}}_{\mathrm{\Phi}}=\mathrm{5}/\mathrm{4}$).
3.1 Friction in the denseflow model
The term τ_{Φ} represents the depthintegrated divergence of the shear stress tensor and thus the constitutive model of the flowing mass. Assuming that the top boundary is stressfree and that surface–tangential derivatives of the deviatoric stress tensor are small, the only remaining entity is the basal friction. In this work, we will use the friction model presented by Rauter et al. (2016), which is closely related to the widely used Voellmy (1955) friction model. It is given as
with dry friction coefficient μ and turbulent friction coefficient χ. A wide range of alternative friction models can be found in the literature, and a number of them are implemented in the presented software.
3.2 Entrainment and deposition in the denseflow model
${S}_{\mathrm{\Phi}}^{\mathit{\varphi}}$ represents the sum of all volumetric source and sink terms of grains (e.g. erosion and entrainment of additional mass or its deposition), and ${S}_{\mathrm{\Phi}}^{u}$ represents its associated momentum. Dividing by the packing density in Eq. (12) simplifies handling of density changes in the different flow regimes. In the simplest case (e.g. laboratory experiments on a nonerodible bed), the source and sink terms are zero.
For snow avalanches and many other realistic gravitational mass flows, entrainment of erodible material along the avalanche path plays an important role. A popular entrainment model can be derived by comparing the dissipated energy in the mass flow with the energy required to mobilize the static material (Fischer et al., 2015):
with the specific erosion energy e_{b} as the single parameter. Here it is assumed that the packing density of the static layer is the same as in the dense flow ϕ_{Φ}.
Rauter and Köhler (2020) presented an extension to account for the deposition of flowing material, ${S}_{\mathrm{\Phi}\to \mathrm{\Sigma}}^{\mathit{\varphi}}$. This aspect is neglected in this work, and the flow height of the last time step is assumed to be the final deposition of the model.
The total flux term between the static layer and the flowing avalanche is determined as the difference between entrainment and deposition:
The related momentum source and sink terms are zero in the case of single layer flows, as both erodible material and deposited material are static.
The height (in surfacenormal direction) of the static material on the topography can be tracked with an additional evolution equation:
with
again under the assumption that the static layer has the same packing density as the flowing avalanche ϕ_{Φ}. Tracking the thickness of the static layer allows the limitation of the available entrainable material; hence we are able to turn off entrainment if the erodible layer is depleted.
The suspension flow model describes the flow of a dynamic mixture of a granular material of density ρ_{g} and the surrounding fluid of density ρ_{c}. It corresponds, to some degree, to a depth integration of the compressible model of Rauter (2021). The mixture density follows as
with the variable packing density or phase fraction ϕ_{Π}. Introducing the buoyant density ratio,
The mixture density can be expressed as
The Boussinesq approximation, an often applied simplification (e.g. Parker et al., 1986), implies that ${\mathit{\varphi}}_{\mathrm{\Pi}}\u2a85{\mathrm{10}}^{\mathrm{2}}$ and r≈1 and thus ρ_{Π}≈ρ_{c}. This is reasonable if ρ_{g} and ρ_{c} are at least similar in order of magnitude (e.g. sand in water). However, this does not hold for snow avalanches, i.e. mixtures of grains or ice (ρ_{g}≈1000 kg m^{−3}) with air (ρ_{c}≈1 kg m^{−3}). Thus, we will omit this assumption and consider the dynamic density as given by Eq. (22) in all terms.
Due to the variable mixture, there will be two phases that have to be described by balance laws. In depthaveraged frameworks, this is usually handled by describing the total volume occupied by the flowing masses (grains and flowing ambient fluid) in terms of the flow depth h_{Π} and the volume of grains, expressed by the depthintegrated volume fraction ${h}_{\mathrm{\Pi}}\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}}_{\mathrm{\Pi}}$ (e.g. Parker et al., 1986). The phases are assumed to move with the same velocity ${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Pi}}$. Differences in velocity (e.g. settling of particles) are considered with empirical corrections.
The depthintegrated mass and momentum conservation equations follow as
All equations and terms are well known from the denseflow model, except for the additional tracking of the grain fraction with Eq. (24). The unknown flow fields are the flow depth h_{Π}, the depthaveraged velocity ${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Pi}}$, and the depthaveraged phase fraction or packing density ${\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}}_{\mathrm{\Pi}}$. Assuming $r\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}}_{\mathrm{\Pi}}\approx \mathrm{0}$ in all terms but the gravitational acceleration (buoyancy assumption) leads to the popular model of Parker et al. (1986). Removing the surface–tangential gravitational acceleration leads to the momentum conservation equation of Bartelt et al. (2016). The effective gravitational acceleration g_{eff} is the surfacenormal gravitational acceleration, corrected for centripetal acceleration due to curved terrain. In terms of surface partial differential equations, it can be easily expressed as (see Appendix A)
This expression replaces the rather complex calculation of the basal pressure in the denseflow model. It is justified here, as the basal pressure has only a weak influence on the flow dynamics of the suspended flow. Further, this notation turns out to be convenient later, as various internal processes in the suspension flow depend on effective gravity. A particle on a streamline of the flow will approximately experience a volume force corresponding to this acceleration, and processes like the terminal settling velocity will depend on this adjusted value.
Considerable attention has to be drawn to the volumetric source and sink terms (${S}_{\mathrm{\Pi}}^{h}$ and ${S}_{\mathrm{\Pi}}^{\mathit{\varphi}}$) and the associated momentum flux (${S}_{\mathrm{\Pi}}^{u}$). These terms are responsible for the varying flow height and the depthaveraged particle volume fraction and influence the flow dynamics substantially.
4.1 Friction in the suspension flow model
Similarly as in the denseflow model, the term τ_{Π} represents the depthintegrated divergence of the shear stress tensor. If the particle fraction in the suspension is low, it can be treated as a simple fluid. The simplest wall friction model can be represented with constant friction coefficient c_{D} (Parker et al., 1986):
However, suspension flows are inherently turbulent, and there is strong evidence that turbulence models with more complex friction relations might be required (e.g. Parker et al., 1986). Nevertheless, we will use the simple model in this work, and we should keep in mind that the wall friction coefficient c_{D} is an empirical parameter that might require adaption to flow conditions. Further, it is assumed that all dissipative processes, such as intergranular friction (e.g. Boyer et al., 2011), are included in this term. Considering the accuracy and uncertainties of the problems at hand, this seems to be a reasonable compromise. Alternative approaches are the turbulence model of Parker et al. (1986), a depth integration of the Einstein viscosity model (e.g. Boyer et al., 2011), and a more complex granular rheology (Boyer et al., 2011).
4.2 Ambient fluid entrainment in the suspension flow model
The volume of the suspension flow will grow due to entrainment of ambient fluid. It is assumed (Parker et al., 1986; Turner, 1986; Ancey, 2004) that ambient fluid entrainment depends solely on the bulk Richardson number, which is given as
In contrast to, for example, Parker et al. (1986), we use the effective surfacenormal acceleration g_{eff} instead of the constant gravitational acceleration $\left\mathit{g}\right$ to account for the influence of centripetal forces on particles in the flow. Adjusting the Richardson number with the centripetal acceleration leads to an increased amount of ambient fluid entrainment if the flow runs over convex terrain and to a decreased amount if the flow runs over concave terrain.
There are various models for the relation between the Richardson number and the entrainment. Parker et al. (1986) use a simple, inverse proportional approach:
with the parameters α=0.00153 and Ri_{0}=0.204.
Turner (1986) provides an alternative formulation:
with the parameters Ri_{0}=0.8, α_{1}=10, and α_{2}=50. Various different parameters were suggested for this empirical relation (see, for example, Ancey, 2004).
Finally, Ancey (2004) suggests yet another relation in the form of an exponential function, here given in the form of Issler et al. (2018):
The parameter α_{1} is supposed to be the only free parameter, with a value of 1.6 following Issler et al. (2018); however, due to different definitions of the entrainment rate, an additional parameter α_{2} is required. In order to be of similar magnitude as the other air entrainment relations, α_{2} has to be roughly 0.05. All relations are shown in Fig. 4.
4.3 Grain entrainment and settlement in the suspension flow model
Suspension flows are, similar to dense flows, able to erode granular material from the bed. It is, in principle, possible to use the same entrainment relations as in the denseflow model, but specialized entrainment relations have been proposed in the literature. An example of subaquatic turbidity currents, is given by Parker et al. (1986) as
with
the settling velocity,
the particle Reynolds number,
the viscosity of the ambient fluid ν_{c}; and two empirical parameters Z_{m}=13.2 and Z_{c}. The parameter Z_{c} was reported to be approximately 5. We found that a value of exactly 0.5 is required to reproduce the examples of Parker et al. (1986) in the examples shown in Sect. 7.2.
The settling of grains is given by Parker et al. (1986) as
with the settling velocity as given in Eq. (34) and the factor r_{0} for the bottom value of the grain concentration,
As before, the total flux term follows as the difference between entrainment and deposition:
The momentum flux into the suspension due to ambient fluid and grain entrainment is zero. The volume occupied by entrained and deposited grains and the respective flux term in the evolution equation of the flow height h_{Π} are neglected at this point.
Granular mass flows can show different regimes, especially in terms of the Stokes number. Sampl and Zwinger (2004) and others (Jóhannesson et al., 2009) describe three regimes: the dense flow, transition or resuspension, and powder snow layer. Sovilla et al. (2015) recognize five regions in mixedsnow avalanches, and Köhler et al. (2018) identified seven regimes. Here we aim to represent the two limit cases of dense flow and suspension in a single model, similar to Bartelt et al. (2016). It is assumed that these regimes are described in appropriate accuracy either by the Savage and Hutter (1989, 1991) model (Eqs. 12 to 14) or the Parker et al. (1986) model (Eqs. 23 to 25). The layers will communicate with mass fluxes S^{ϕ} and S^{h} and momentum fluxes S^{u}. In particular, the fluxes of grains are (see also Fig. 1) for the static layer,
for the denseflow layer,
and for the suspension layer
Entrainment by the suspension layer is assumed to be negligibly small in comparison to the overall mass fluxes and thus is not explicitly accounted for in the simulations. The term ${S}_{\mathrm{\Phi}\to \mathrm{\Pi}}^{\mathit{\varphi}}$ describes the upward mass flux from the dense flow to the suspension flow. It is the remaining term to be specified in the following (see Sect. 5.1). The flux in the opposite direction ${S}_{\mathrm{\Pi}\to \mathrm{\Phi}}^{\mathit{\varphi}}$ is assumed to be equal to the settling flux of the suspension layer ${S}_{\mathrm{\Pi}\to \mathrm{\Sigma}}^{\mathit{\varphi}}$; i.e. the deposition from the suspension is redirected to the dense core and further to the static layer from there, if the deposition model of the denseflow model is active. The corresponding momentum fluxes for the denseflow layer and the suspension layer are
accounting for the momentum that is transferred together with grains between moving layers. The shape factor ξ_{t} takes into account that the velocity at the top boundary of the avalanche, where particles are tossed into the suspension layer, is higher than the depthintegrated velocity. It is related to the previously shown shape factor and can similarly be calculated on the basis of, for example, the Bagnold (1954) velocity profile as $\mathrm{5}/\mathrm{3}$. The particles that fall from the suspension layer onto the denseflow layer, ${S}_{\mathrm{\Pi}\to \mathrm{\Phi}}^{\mathit{\varphi}}$, are assumed to carry the velocity of the suspension layer. While the exact form of this flux is not important, it is vital to remove momentum together with mass, so as to not increase the velocity of the remaining mass to potentially very high values. The momentum fluxes to and from the static layer are zero due to the respective velocity at the interface.
Further we have to account for the volume of fluid that is pushed into the suspension layer with particles. Assuming that particles enter at a packing density of ϕ_{0Π}, we have to add a source term of the form
The value ϕ_{0Π} is set to the phase fraction of the dense core in this work. This avoids unreasonably high grain fractions if a suspension flow is initiated by a denseflow avalanche.
In addition to the particleborne momentum fluxes, we need to consider the shear stress at the interface. This relation is chosen to be identical to the basal shear stress of the suspension layer, τ_{Π}; however, it is no longer proportional to the velocity of the suspension layer but to the relative velocity between the denseflow and the suspension layer:
In areas where the suspension layer detaches from the dense flow, the denseflow velocity is assumed to be zero, and the model collapses to the friction model of the ordinary suspension model. An equal but opposite stress term to τ_{Π} should be applied to the dense core to account for the friction of the top surface of the dense core. However, it is assumed that this stress is already included in the empirical formulation and parameterization of τ_{Φ}, because the top surface friction is also present in pure densesnow avalanches with a stationary or moving air layer above it. The ambient fluid entrainment of the suspension layer stays unchanged.
The mass flux ${S}_{\mathrm{\Phi}\to \mathrm{\Pi}}^{\mathit{\varphi}}$ feeds the suspension layer from the dense core; the associated momentum flux, in combination with the shape factor, propels the suspension flow forwards. This is assumed to be the mayor genesis mechanism for the suspension cloud, as in most other avalanche models.
5.1 Crosslayer coupling
All fluxes of the twolayer model are described relatively well in the literature (see sections above), except for the mass flux from the denseflow layer to the suspension layer, ${S}_{\mathrm{\Phi}\to \mathrm{\Pi}}^{\mathit{\varphi}}$, for which only few suggestions can be found (Eglit, 1998; Issler, 1998; Sampl and Zwinger, 2004; Bartelt et al., 2016). Existing relations do not conceptually fit into the presented framework, either due to missing granular mechanics (Eglit, 1998; Issler, 1998; Sampl and Zwinger, 2004) or due to their dependence on a specific denseflow model (Bartelt et al., 2016). For the purpose of introducing this framework, we choose a simple relation, based on local flow fields of the dense flow.
We assume that the dense flow is composed of small and large particles with diameters d_{Π} and d_{Φ}, respectively. Uptake of particles into the suspension layer requires small particles to be made available by the dense layer that mostly consists of large particles (Bartelt et al., 2016) and the capability of the suspension layer to keep them suspended. The latter is already implemented into the model in the form of the settling model of the suspension flow ${S}_{\mathrm{\Pi}\to \mathrm{\Phi}}^{\mathit{\varphi}}$. This term depends on the particle Reynolds number Re_{g}, which is similar to the Stokes number and a good indicator of the flow regime.
The first step, making small particles available to the suspension, is assumed to be triggered by a fluidized flow that expands in volume, sucks in air, and increases the distance between particles. There are various hints on how this expression should look. At first it is useful to look at dimensionless properties in the dense flow. Besides the nondimensional volumetric mass flux ${S}_{\mathrm{\Phi}\to \mathrm{\Pi}}^{\mathit{\varphi}}/\left{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Phi}}\right$, these are the friction coefficient ${\mathit{\mu}}_{\mathrm{\Phi}}=\left{\mathit{\tau}}_{\mathrm{\Phi}}\right/{p}_{\mathrm{\Phi}}$, the packing density ϕ_{Φ}, and the inertial number (Forterre and Pouliquen, 2008; Rauter, 2021):
with the shear rate at the bottom of the dense flow, which is assumed to follow Bagnold (1954):
It is well established that μ_{Φ} and ϕ_{Φ} can be expressed as a function of only the inertial number I_{Φ}, and it is reasonable to assume that fluidization can be described the same way. This is further emphasized by the linear relationship between the packing density and the inertial number in the denseflow regime (Forterre and Pouliquen, 2008). Finally, Rauter et al. (2016) found a specific relation between the shear rate ${\dot{\mathit{\gamma}}}_{\mathrm{\Phi}}$ and the pressure p_{Φ} in a granular kinetic theory model (Vescovi et al., 2013) at the point where fluidization suddenly occurs:
Comparing this relation to the expression for the inertial number, one can observe a striking resemblance; solely the exponent of the pressure is slightly lower in the relation of Rauter et al. (2016). This strongly indicates that the mass flux from the dense flow to the suspension can be expressed as a function of the inertial number only, starting at a minimum value I_{0} and growing with a specified rate s_{f} from thereon:
The results of Rauter et al. (2016) suggest that the value of I_{0} is close to 0.5, as at this point explosive fluidization starts to occur. The factor s_{f} is expected to be small, as the vertical velocity has to be substantially smaller than the flow velocity. This parameter can be optimized to yield the correct relation between dense flow and powder cloud.
In this model, small particles will be made available to the suspension when the denseflow velocity is high or when the pressure is low (e.g. when an avalanche is running over a bump). If the small particles are sufficiently small, the suspension will be able to keep the particles suspended and a powder cloud will form. Otherwise, the particles will fall back to be reintegrated by the dense core, expressed by the deposition mass flux of the suspension layer, which is stronger for larger particles. The parameters for the suggested model are the small and large particle diameters d_{Π} and d_{Φ}, the minimum value of I at which fluidization occurs I_{0}, the particle density ρ_{g}, and the factor s_{f}. All parameters except for the latter are already used in the model or known otherwise.
Relation (48) finally completes the model and closes the system that will be solved numerically in the following. The model could be improved by tracking and limiting the availability of small particles or by making this property temperaturedependent.
The pre and postprocessing of simulations with the presented models follows the workflow depicted by Rauter et al. (2018). The capabilities of the respective tools have been improved and fully implemented in C++ to allow a seamless integration into OpenFOAM and computational clusters that do not support Python and some of the previously used libraries. Most improvements are based on a native implementation of two common geographic information system (GIS) data types: ESRI^{®} shape files and ESRI^{®} grid files. The native implementation allows all solvers and utilities of the OpenFOAM avalanche module to directly read and write to or from the respective files. This enables many previously difficult tasks that are presented in the following. Generally, all tools are steered with text files that follow the usual OpenFOAM syntax, called dictionary (see Fig. 5). This toolchainbased workflow follows the concept of OpenFOAM, which ensures reproducibility and facilitates the reuse of existing code and the rapid development of new code.
6.1 Mesh generation from terrain data
The mesh generation follows the principles of Rauter et al. (2018). In a first step a triangulation of the terrain and a boundary of the surrounding volume is generated. A new tool for this task, called gridToSTL
, was written entirely in C++ and without any external dependencies. The tool requires input in the form of a polygon that defines the simulation domain and the terrain data in the form of a raster file. In contrast to the previous version, the polygon can be any kind of closed and nonintersecting polygon with an arbitrary number of edges, either convex or concave. This enables flexibility in the choice of the simulation domain, which turned out to be especially useful to cover long and windy submarine canyons.
The finite volume mesh is generated from the triangulated surface with an arbitrary mesh generator. This toolchain can be applied to the depthintegrated models presented here and was also used for the full threedimensional model presented by Rauter et al. (2022). In this study we used the mesh generator pMesh
, while Rauter et al. (2022) used cartesianMesh
, both from the cfMesh toolbox (Juretić, 2015). The finite area mesh is then generated on a dedicated surface of the finite volume mesh using the tool makeFaMesh
, part of the OpenFOAM finite area module.
6.2 Mapping initial conditions
Initial conditions can be set with the tool releaseAreaMapping
. In addition to the functionality of previous versions, this tool is now able to read shape files and grid files and map them directly onto finite area fields to be used by any solver. All input for the tool is read from a dictionary, where further references to shape and grid files can be listed. This tool enables efficient adaption to new scenarios.
6.3 Simulation run
Once the mesh and the initial conditions are defined, the solver of choice can be run. Currently there are three solvers available in the avalanche module, the denseflow solver faSavageHutterFoam
, the suspension flow solver faParkerFukushimaFoam
, and the mixedflow solver faTwoLayerAvalancheFoam
(Fig. 5 shows faTwoLayerAvalancheFoam
only, but it can be replaced with any of the other two models). Physical parameters are read from the file transportProperties
, general simulation settings from controlDict
, and numerical algorithms and parameters from the files faSolution
and faSchemes
. To run the solver in parallel, the tool decomposePar
has to be run before the solver, and the tool reconstructPar
has to be run after the solver. In the common OpenFOAM manner, all steps for a simulation are listed in a script file named Allrun
, which can be launched by the user to automatically execute the pipeline outlined here. Another script, named Allclean
, can be run to clean up the simulation directory.
6.4 Postprocessing and data export
The OpenFOAM architecture allows the execution of customized code, called function objects, in every simulation step. Various function objects are made available in the avalanche module. Most importantly, this includes function objects to export simulation results as either shape or raster files. The export as shape files can be done cellwise (one polygon for each computational cell), or the numerical data can be recombined to generate isolines that are written into the shape file. Function objects can be loaded by placing the respective entry into the control dictionary. As of version v2312, all solvers are able to run in a postprocessing mode, in which old results are read from a hard disc and the function objects are executed. This allows the execution of function objects in a postprocessing workflow without rerunning the whole simulation.
7.1 Denseflow model
The denseflow model has been applied to various cases in multiple studies. The interested reader is referred to Rauter and Tuković (2018) for lab scale simulations, Rauter et al. (2018) and Huber et al. (2018) for largescale snow avalanche simulations, Rauter and Köhler (2020) for simulations with the deposition model, and to Shimizu (2022) for an application to pyroclastic flows.
7.2 Suspension flow model
Parker et al. (1986) simulate steady suspension flows on constantly inclined onedimensional slopes with the model presented in Sect. 4. Four cases with uniform model parameters but different boundary conditions give a good overview of the behaviour of the model and a verification (as defined by Roache, 1997, as solving the equations right) of the presented implementation. The four simulations are conducted on onedimensional slopes with a gradient of 5 %. The gravitational acceleration follows as $\mathit{g}=(\mathrm{0.49},\mathrm{0},\mathrm{9.81}{)}^{T}$ m s^{−2} (chosen to match the setup by Parker et al., 1986). The parameters suggest that the suspensions are composed of sediment in water on a scale of a small turbidity current.
Material parameters for this setup are given in Table 1. The left boundary condition (at x=0) prescribes the inflow in terms of the height h_{Π}, velocity ${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Pi}}$, and grain flux ${\stackrel{\mathrm{\u203e}}{\mathit{\psi}}}_{\mathrm{\Pi}}={h}_{\mathrm{\Pi}}\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}}_{\mathrm{\Pi}}\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{\mathrm{\Pi}}$, in particular as shown in Table 2. All parameters are given normalized to reference values H=2 m, U=0.874 m s^{−1}, and Ψ=0.00828 m^{2} s^{−1}.
The right boundary condition is modelled as a zero gradient for all fields, mimicking an outlet boundary condition. For a basic verification of the novel implementation of the suspension model, the respective simulations are repeated and compared to the original results. We will evaluate the buoyancy assumption of Parker et al. (1986), as well as the formulation with the correct density given in here. The simulations are conducted in an unsteady manner until the flow reaches a steady state, comparable to the results reported by Parker et al. (1986). Figure 6 shows results for the four cases.
The first case, starting with a high velocity but low particle fraction increases its particle fraction quickly, as the high velocity is sufficient to erode and pick up sediment. The second case starts with a very high phase fraction, leading to a sudden ignition of the flow at $x/H=\mathrm{60}$. The height of the suspension stays low and even decreases, showing that a high phase fraction can keep the suspension concentrated at the bottom. The third and fourth case start with a low velocity and low particle phase fraction, respectively, and the suspension fades out quickly. The height of the flow increases in both cases where the flow fades out, indicating that the momentum of the flow is diffused over larger volumes of fluid. This is consistent with the expected scaling of fluid entrainment with the Richardson number.
It can be seen that results of Parker et al. (1986) are reproduced with only small deviations. The OpenFOAM solver yields sharper edges than the implementation of Parker et al. (1986), especially visible in Fig. 6b. This small difference is most likely attributable to the numerical solution method or the numerical resolution. The correction of the time derivative and advection term with $(\mathrm{1}+r\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}}_{\mathrm{\Pi}})$ has only a minimal influence on the model results. This is reasonable, considering the low value for the buoyant density ratio r=1.65 in these cases. These simulations suggest that the model of Parker et al. (1986) is implemented correctly. However, this can not be seen as a validation (Roache, 1997) of the model.
7.3 Twolayer model
7.3.1 Synthetic tests and sensitivity study
In order to better understand the twolayer model, we will conduct tests on synthetic topographies. The topography is based on a parabola with a length of L=4000 m and a height of H=2000 m, with an additional flat runout area of 2000 m. The slope has a width of 2000 m, leading to a simulated region of $x=[\mathrm{4000},\mathrm{2000}]$ m and $y=[\mathrm{1000},\mathrm{1000}]$ m. In addition, the influence of topographic structures will be investigated, as terrain features often initialize the formation of suspension flows (e.g. powder snow avalanches). A bump in the surface is created by superposing the parabola with a Secans hyperbolicus $\mathrm{sech}\left(x\right)=\mathrm{2}/\left(\mathrm{exp}\left(x\right)+\mathrm{exp}(x)\right)$ at ${X}_{\mathrm{p}}=\mathrm{2700}$ m with a height of H_{p}=150 m and a length of L_{p}=200 m:
inspired by the experiments of Viroulet et al. (2017). All boundaries are implemented as von Neumann (zero gradient) boundary conditions.
The release area (initial condition) of the slide was formed by a rectangle between $x=[\mathrm{3900},\mathrm{3500}]$ m and $y=[\mathrm{500},\mathrm{500}]$ m and an initial denseflow height of h_{Φ}=5 m within that square. All other flow fields are set to zero. The parameters, roughly corresponding to snow avalanches, are given in Table 3, if not mentioned otherwise. The value for the coupling factor s_{f} is varied, and the sensitivity of the model to this parameter is investigated. Entrainment and deposition to and from the static layer are not included in this section for simplicity. The simulations were run for 90 s.
Besides the flow thickness, velocity, and phase fraction, we can analyse the dynamic pressure, which is an important indicator of the destructive potential of the flow. It is defined as
for the dense flow and as
for the powder cloud (e.g. Jóhannesson et al., 2009). In particular we evaluate the dynamic peak pressure, which is defined as the maximum of the dynamic pressure at a fixed point over time. Important limits that are used in the definition of hazard zones (e.g. in Austria) are 1 kPa (yellow zone) and 10 kPa (red zone) (Jóhannesson et al., 2009). Notably, the shape factor should be applied to the dynamic pressure for consistency, increasing all simulated pressures by 25 %. However, this is neglected in order to be consistent with previous work and the definition of hazard zones.
Results for a simple parabola (without surface bump) are shown in Fig. 7 for three values of s_{f} (10^{−5}, 10^{−4}, 10^{−3}). This set of simulations allows some valuable conclusions on the model and in particular the coupling model. All simulations start with a dense flow that eventually feeds the powder cloud. The feeding rate of the powder cloud varies strongly due to the variation of the parameter s_{f}.
For a low value of s_{f}, the dense flow is not able to generate a strong powder cloud with a considerable particle phase fraction and thus density. A suspension flow develops eventually; however, it consists almost entirely of air, without any ice particles. Basically, this can be seen as a layer of air that is dragged along by the dense flow. The velocity, dynamic pressure, and runout distance of this layer are correspondingly low. As shown before, the flow height of the suspension layer grows strongly for fading flows, indicating a strong diffusion of momentum.
Increasing the value for s_{f} up to 10^{−4} leads to higher phase fractions up to 0.004, roughly corresponding to a density of 4 kg m^{−3}. Further increasing the value to 10^{−3} leads to phase fractions of up to 0.02 and densities of 20 kg m^{−3}, however only for short periods. Notably, these are depthaveraged phase fractions and densities, and the respective values close to the surface might be considerably higher. The corresponding dynamic pressure of the powder cloud is still low, and only the simulation with the highest coupling factor s_{f} is able to generate a red zone that extends beyond the red zone of the dense flow. These results seem reasonable, considering the average slope gradient of 50 % and the absence of any topographic features that might enhance the feed of the powder cloud. More powerful powder snow avalanches can be expected on steeper slopes and on slopes with high topography variations (e.g. steep cliffs or rough terrain). Further simulations (not shown here) revealed that the powder cloud increases substantially with higher slope gradients.
Results for the slope with a bump are shown in Fig. 8. The model shows a high sensitivity to the terrain, and this case represents natural slopes with varying gradients better. All simulations create a considerable powder cloud with high phase fractions. The highest phase fraction is reached shortly after the top of the bump where the negative centrifugal forces are strongest and the basal pressure the lowest. The phase fraction reaches up to 0.05, roughly corresponding to a density of 50 kg m^{−3}. A shock is formed at the bump in the suspension layer due to the high gradient in the phase fraction, leading to a considerable pressure gradient that decelerates the flow. In all simulations the dynamic powder cloud pressure exceeds 10 kPa, and the respective high pressure zone extends beyond the denseflow runout. The 1 kPa zone of the powder cloud reaches considerable runouts beyond the dense flow.
The results on synthetic terrain show a reasonable behaviour of the model, both in terms of parameterization and response to the terrain. The effect of the terrain is well visible and corresponds to the assumptions from which the model was derived. The sensitivity of the model to the parameter s_{f} is well pronounced, and this factor can be utilized to fit the model to realworld observations. A value between 10^{−4} and 10^{−5} seems reasonable for the parameter of the coupling model.
Finally, we use the synthetic cases to showcase the sensitivity of the model to the air entrainment. Figure 9 shows the simulation on the synthetic terrain with the three presented air entrainment models. The differences are small but noticeable. In particular, the entrainment is stronger with the model of Ancey (2004), however, which is just a question of parameterization. More importantly, the model of Turner (1986) shows a more pronounced flow head. The Richardson number is low in the head, and the relation of Turner (1986) predicts the strongest entrainment at low Richardson numbers (see Fig. 4). Generally, all relations appear reasonable and well in line with each other. We will continue with the entrainment model of Parker et al. (1986) from here on. Considering an optimization of air entrainment parameters to real events, it might be useful to apply the model of Ancey (2004) instead, as it provides the clearest parameterization.
7.3.2 Realcase example: the 1988 Wolfsgruben avalanche
The 1988 Wolfsgruben avalanche represents an important event in Austria, as it was the trigger for many developments and used repeatedly as a benchmark. The event, or at least its dense core, was featured by Fischer et al. (2015) and Rauter et al. (2018). Here we revisit the event with the new twolayer model and include the powder cloud into the analysis. The avalanche is characterized by a channellized, steep slope with an angle of 30° that transitions fairly abruptly to the flat valley floor and the opposite slope.
The preprocessing and simulation setup follows Rauter et al. (2018) but with the novel tool chain and an extended simulation domain to cover the full runout of the powder cloud. The initial release areas of the avalanche and the erodible snow covers are the same, following the linear approach:
where z is the surface elevation and z_{0} the elevation of a reference point with the base value h_{Σ}(z_{0}). The growth rate $\frac{\partial {h}_{\mathrm{\Sigma}}}{\partial z}$ defines the change with elevation from that point. θ is the angle between the gravitational acceleration and the surfacenormal vector. For the 1988 Wolfsgruben avalanche, we use the snow cover parameters h_{Σ}(z_{0})=1.61 m, z_{0}=1289 m, and $\frac{\partial {h}_{\mathrm{\Sigma}}}{\partial z}=\mathrm{8}\times {\mathrm{10}}^{\mathrm{4}}$.
The model parameters are shown in Table 3. The denseflow parameters have been optimized in a previous study (Fischer et al., 2015; Rauter et al., 2018), however for a model ignoring the powder cloud. The addition of the powder cloud could lead to different optimal denseflow parameters. However, that is out of the scope of this work, and we assume that the previously used parameters fit sufficiently well. The suspension parameters are deduced from literature where possible (density, grain diameter). The coupling parameter ${s}_{\mathrm{f}}={\mathrm{10}}^{\mathrm{5}}$ was found after running some simulations, starting from the values derived from the simulations on synthetic cases. A higher value leads to an unrealistically short denseflow runout and a lower value to a severe underestimation of the suspension impact pressure. The friction coefficient c_{D} was chosen to be sufficiently large for the powder cloud to not completely decouple from the dense core. Apart from this effect, the simulation is rather insensitive to the friction coefficient c_{D}.
Four time steps of the simulation are shown in Fig. 10, displaying isolines of the denseflow height h_{Φ} (a–d) and the suspension flow height h_{Π} (e–h). The avalanche starts as a dense flow and rapidly accelerates due to the steep release area (Fig. 10a). Shortly after the release a strong suspension layer is formed that further accelerates beyond the velocity of the denseflow layer (Fig. 10b, f). After roughly 40–50 s the avalanche reaches the bottom of the valley (Fig. 10c, g). The powder cloud outruns the dense flow and hits the valley floor first. The dense flow is stopped quickly due to the high granular friction, while the powder cloud keeps running up on the opposite slope for approximately 50 m of elevation. Both flows experience a shock that increases the flow height in the valley floor drastically. The deposition (i.e. the denseflow height in the last time step) reaches up to 15 m, however, which does not account for the difference between flow (≈200 kg m^{−3}) and deposition density (≈600 kg m^{−3}).
Results for the dense flow can be validated by a comparison with the deposition (see Fig. 11b), and they are similar to previous studies with the same model and the model SamosAT (Fischer et al., 2015; Rauter et al., 2018). Results for the powder cloud are more difficult to validate. Traces of the powder cloud that can be identified in the field are limited and not straightforward to interpret, as no clear deposition pattern emerges from suspended flows. Further, the respective deposition can hardly be related to the impact pressure and thus the destructive potential of the flow. Therefore, we compare the simulated dynamic pressure with observed building damages from the respective avalanche (see Fig. 11a). This includes not only the suspension layer but also the dense flow. An evaluation of the dynamic peak pressure and the deposition height at damaged objects is shown in Table 4.
The dense flow does not reach the two destroyed buildings (Point 1 and 2 in Fig. 11a and Table 4) and stops about 20 m short. Points 12 and 18 were only slightly damaged by the suspension flow in reality but severely hit by the dense flow in the simulation, showing that the simulation tends too strongly to the left side (viewed in the flow direction). The deposition height that was recorded at selected points (Table 4) is matched well, assuming a compaction of the avalanche by a factor of 3 after deposition.
The suspension layer shows a very limited zone of high dynamic pressure (>10 kPa) but an extended zone of intermediate dynamic pressure (1–10 kPa). The model predicts dynamic pressures of 1–4 kPa where balconies and roofs have been damaged and 1–3 kPa where windows have been destroyed. This corresponds well with engineering estimations of resistance capabilities of the respective parts: windows are assumed to break at 2–4 kPa and doors, walls, and roofs at 3–6 kPa (Sovilla et al., 2015). The dynamic pressure of the suspension layer at the destroyed buildings (Point 1 and 2 in Fig. 11) is not sufficient to destroy the respective brick structures (25–45 kPa). These high values strongly indicate that the dense flow or a saltation layer must be responsible for these high impact pressures (see pictures in Fischer et al., 2015). Therefore, we conclude that the simulated suspension layer reaches all observed traces of the powder cloud without covering the region where no traces could be observed.
7.3.3 Realcase example: the 2019 Eiskar avalanche
On 15 January 2019, the Eiskar avalanche occurred after intense snowfall and a rapid temperature drop (Oesterle, 2019). The Eiskar avalanche differs substantially from the Wolfsgruben avalanche regarding topography and thus provides a good supplement to that case. The avalanche was initiated by a slab on the righthand side of the avalanche path (looking in the flow direction) that fell onto a larger snow field. From there, the avalanche slope continues with an inclination of approximately 25° for 1500 m until reaching a flatter slope of 10°. The denseflow avalanche ran 1000 m on the flattening slope, and the powder flow exceeded the dense flow by another 500 m, reaching the village of Ramsau. The powder cloud destroyed a wooden building, damaged a hotel, and knocked over a bus. The dynamic pressure required for the damage was estimated at 1–3 kPa. Areal pictures were taken after the event, which allowed the estimation of the initial snow cover, the release area, and the deposition. The data were used to derive parameters for the snow cover function (Eq. 52), h_{Σ}(z_{0})=1.60 m, z_{0}=1275 m, and $\frac{\partial {h}_{\mathrm{\Sigma}}}{\partial z}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{3}}$, to reach a snow cover thickness of approximately 2.7 m at an elevation of 2200 m (Oesterle, 2019). Other aspects of the simulation, such as the preparation of the terrain data, match the simulation of the Wolfsgruben avalanche.
The first simulation (not shown) was conducted with the same parameters as for the Wolfsgruben avalanche. However, these parameters lead to a severe underestimation of the powder cloud, running short by approximately 400 m. Simulations with the model SamosAT (Sampl and Zwinger, 2004) showed similar results with the standard parameters (Oesterle, 2019). Therefore, the friction coefficients and the coefficient for the suspension feed were adjusted (see Table 3) to reach an appropriate runout and dynamic pressure at the observed impacts.
Five time steps of the simulation are shown in Fig. 12 in terms of the denseflow height h_{Φ} (a–e) and the suspension flow height h_{Π} (f–j). The collapsing slab (Fig. 12a) falls down the steep cliff onto a larger snow field where it can entrain additional snow. After around 30 s the avalanche reaches a second cliff, and a powder cloud starts to emerge (Fig. 12b, g). The suspension layer keeps growing substantially in the slope section with an inclination of 25° (Fig. 12c, h) and starts to detach when reaching the flatter slope of the 10° inclination. The suspension layer reaches the village of Ramsau after approximately 90 s (Fig. 12i, j), while the dense flow comes to a halt at the exit of the valley (Fig. 12d, e). Interestingly the dense flow is pushed towards the left by terrain features at the exit of the valley, while the suspension layer is essentially unaffected by these small obstacles.
The corresponding zones of dynamic pressure are shown in Fig. 13a. The 1 kPa isoline of the suspension layer extends far into the village. This fit was used as a benchmark to determine the optimal model parameters and thus matches observations well. The final deposition of the model is shown in Fig. 13b and compared to the observed deposition. The observations distinguish between suspension and denseflow depositions (Oesterle, 2019), and the same can be done in the numerical model. The denseflow layer leaves behind up to 10 m thick deposits (to be corrected by a factor of $\mathrm{1}/\mathrm{3}$ to match the deposition density) with sharp edges, while the suspension generates deposits with 0.1–0.2 m thickness (to be corrected as well) that fade out gradually. Both the position of the respective deposits and their shape match the observations.
Overall, the model is able to reproduce the observed flow traces, from the dynamic pressure to the varying snow deposits, in a single simulation. However, the model parameters had to be fitted to achieve these results. The friction parameters have to be substantially lower than in the Wolfsgruben case, and the coupling factor has to be an order of magnitude higher. This indicated that conditions were substantially different between the two cases or that the model does not cover some important processes.
This work provides an overview of the implementation of the granular denseflow model of Savage and Hutter (1989, 1991) and the suspension flow model of Parker et al. (1986) into OpenFOAM. Further, the models have been combined by means of a novel coupling mechanism to provide a simple yet effective mixedsnow avalanche model. These three models form the core of the OpenFOAM avalanche module. The module is accompanied by a new toolchain that substantially simplifies the practical application of the framework. The integration of geographic information system (GIS) file types into the OpenFOAM framework enables a simple and deep integration in existing workflows. Moreover, the dependencies on thirdparty libraries for GIS support were removed as they were often missing on computational clusters. In comparison to the work of Rauter et al. (2018), the models and all tools are integrated into OpenFOAM to simplify installation. The physical models are highly modular. Tweaking and replacing specific empirical relation or process models are a core feature of the framework and highly encouraged.
The implementation of the suspension flow model of Parker et al. (1986) was verified by repeating published results, assuring the absence of implementation errors. A novel twolayer model was developed and evaluated with simple synthetic cases. The results are reasonable and follow the expectations set in the model. Further investigations have been conducted with two different realcase avalanches. The reach of the denseflow layer and the suspension layer matched the observed runout in both cases with good accuracy, although a quantitative comparison was not conducted. The dense flow of the Wolfsgruben avalanche came short of approximately 20 m. The impact pressure of the suspension flow is reasonable considering the observed damage. Results for the Eiskar avalanche similarly match observations well if the parameters are fitted accordingly.
The good results are strongly linked to the parameterization, which is highly uncertain due to the limited experience with mixedsnow avalanche models in general and this model in particular. A wide variety of results can be achieved by tweaking the parameters of the model, and substantial investigations will be required to find the appropriate parameters for the large number of semiempirical relations embedded in the flow models. Substantially different parameters were required to yield reasonable results in both cases, a wellknown problem in gravitational mass flow modelling (Scheidegger, 1973; Lucas et al., 2014). Further, snow properties and temperatures might have been substantially different between the two avalanche events. In this regard we see a strong opportunity to substantially improve the twolayer model. Temperature has a strong influence on the particle diameter distribution in snow avalanches and will thus have a high effect on the mobility and the ability to generate suspension flows (Steinkogler et al., 2015a, b).
The denseflow runout and especially its dynamic pressure at a specific point are very sensitive to the parameters. This is related to the strong friction that rises rapidly in flat regions, where the driving gravitational acceleration also vanishes. The suspension cloud is less sensitive to such influences, as the friction is lower and independent of the inclination and the basal pressure. Therefore, the suspension runout is less sensitive to the parameters.
For practical applications we advise to use the existing guidelines for the denseflow parameters (e.g. Salm et al., 1990). For snow avalanches with a high potential to generate powder snow clouds, we suggest to apply the suspension and coupling parameters as used in the Eiskar case. It should be noted that the suspension model absorbs mass from the denseflow model, which reduces the respective runout. Therefore, it might be reasonable to simulate scenarios with less powder flow generation to not underestimate the runout of the dense core. Finally, it should be kept in mind that the results of the model are accompanied by a high number of uncertainties and that they should be used accordingly. Nevertheless, the simulations presented here recreate the processes of the events well and provide a considerable amount of additional information.
Generally, the model and the whole framework aim to be very flexible to provide researchers with a strong platform to develop and evaluate novel friction, entrainment, and coupling models. The introduced coupling model represents a reasonable approach that yields promising results, but there might be big opportunities for improvement. We hope that the framework can provide a starting point for other researchers to develop new coupling mechanisms with better performance. Further, new solvers can be implemented on the basis of the framework, e.g. multiphase models for debris flows (e.g. Pudasaini, 2012; Kowalski and McElwaine, 2013; Iverson and George, 2014) as done by Garcés et al. (2023) with faDebrisFoam
or landslide tsunamis (e.g. George et al., 2017). The toolchain and postprocessing routines presented here can be reused with these models, and additional pre and postprocessing utilities can be added to enlarge the functionality of the whole framework.
The basal pressure is computed following Eq. (14) in the model of Rauter and Tuković (2018). For the Parker et al. (1986) model we tried to achieve a simpler model that can also be combined with the empirical process models in the powder cloud but still follows the general approach. Neglecting the small longitudinal pressure gradient term and removing the indices marking the layer, Eq. (14) can be simplified to
We want to compare this equation to the following equation with an effective gravitational acceleration that contains the effects of centrifugal forces:
Setting n^{Γ} p_{Φ} in Eqs. (A1) and (A2) equal to one another yields
After approximating ${\mathbf{\nabla}}_{\mathrm{n}}^{\mathrm{\Gamma}}\cdot \left(h\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{\mathit{u}}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{\mathit{u}}\right)$ as $h\phantom{\rule{0.125em}{0ex}}{\mathbf{\nabla}}_{\mathrm{n}}^{\mathrm{\Gamma}}\cdot \left(\stackrel{\mathrm{\u203e}}{\mathit{u}}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{\mathit{u}}\right)$ and dividing by h ρ, we can write
A further approximation neglects the shape factor ξ, finally leading to the effective gravitational acceleration as described by Eq. (26).
The code is available in the OpenFOAM avalanche repository at https://develop.openfoam.com/Community/avalanche (Rauter, 2024) under the tag v2312. It is further included in the OpenFOAMv2312 builds and releases (https://www.openfoam.com/news/mainnews/openfoamv2312, ESIOpenCFD team, 2024). The 1988 Wolfsgruben avalanche simulation and previous test and validation cases are included as a tutorial in the repository. The code is licensed under GNU General Public License v3. Test data are licensed under CC BY 3.0 by Amt der Tiroler Landesregierung (AdTLR).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1765452024supplement.
MR: conceptualization, methodology, software, and writing. JK: conceptualization, validation, and writing.
The contact author has declared that neither of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank Matthias Granig and Felix Oesterle (WLV) for providing documentation, support, and data for the realcase examples.
This openaccess publication was funded by RWTH Aachen University.
This paper was edited by Thomas Poulet and reviewed by Dieter Issler and one anonymous referee.
Ancey, C.: Powder snow avalanches: Approximation as nonBoussinesq clouds with a Richardson number–dependent entrainment function, J. Geophys. Res.Earth, 109, F01005, https://doi.org/10.1029/2003JF000052, 2004. a, b, c, d, e
Bagnold, R. A.: Experiments on a gravityfree dispersion of large solid spheres in a Newtonian fluid under shear, Proc. Roy. Soc. Lond. A, 225, 49–63, https://doi.org/10.1098/rspa.1954.0186, 1954. a, b, c
Barker, T. and Gray, J. M. N. T.: Partial regularisation of the incompressible μ(I)rheology for granular flow, J. Fluid Mech., 828, 5–32, https://doi.org/10.1017/jfm.2017.428, 2017. a
Barker, T., Rauter, M., Maguire, E., Johnson, C., and Gray, J.: Coupling rheology and segregation in granular flows, J. Fluid Mech., 909, A22, https://doi.org/10.1017/jfm.2020.973, 2021. a
Barré de SaintVenant, A. J. C.: Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leurs lits, CR Acad. Sci., 73, 237–240, 1871. a
Bartelt, P. and McArdell, B. W.: Granulometric investigations of snow avalanches, J. Glaciol., 55, 829–833, https://doi.org/10.3189/002214309790152384, 2009. a
Bartelt, P., Buser, O., Vera Valero, C., and Bühler, Y.: Configurational energy and the formation of mixed flowing/powder snow and ice avalanches, Ann. Glaciol., 57, 179–188, https://doi.org/10.3189/2016AoG71A464, 2016. a, b, c, d, e, f, g, h, i, j, k
Bouchut, F. and Westdickenberg, M.: Gravity driven shallow water models for arbitrary topography, Commun. Math. Sci., 2, 359–389, 2004. a, b
Bouchut, F., MangeneyCastelnau, A., Perthame, B., and Vilotte, J.P.: A new model of Saint Venant and Savage–Hutter type for gravity driven shallow water flows, C. R. Math., 336, 531–536, https://doi.org/10.1016/S1631073X(03)001171, 2003. a, b
Boyer, F., Guazzelli, É., and Pouliquen, O.: Unifying suspension and granular rheology, Phys. Rev. Lett., 107, 188301, https://doi.org/10.1103/PhysRevLett.107.188301, 2011. a, b, c, d
Christen, M., Kowalski, J., and Bartelt, P.: RAMMS: Numerical simulation of dense snow avalanches in threedimensional terrain, Cold Reg. Sci. Technol., 63, 1–14, https://doi.org/10.1016/j.coldregions.2010.04.005, 2010. a, b
Craster, R. V. and Matar, O. K.: Dynamics and stability of thin liquid films, Rev. Modern Phys., 81, 1131–1198, https://doi.org/10.1103/RevModPhys.81.1131, 2009. a
Denlinger, R. P. and Iverson, R. M.: Granular avalanches across irregular threedimensional terrain: 1. Theory and computation, J. Geophys. Res.Earth, 109, F01014, https://doi.org/10.1029/2003JF000085, 2004. a
Dziuk, G. and Elliott, C. M.: Finite element methods for surface PDEs, Acta Numer., 22, 289–396, https://doi.org/10.1017/S0962492913000056, 2013. a
Eglit, M.: Mathematical and physical modelling of powdersnow avalanches in Russia, Ann. Glaciol., 26, 281–284, https://doi.org/10.3189/1998AoG261281284, 1998. a, b
Eglit, M., Yakubenko, A., and Zayko, J.: A review of Russian snow avalanche models–from analytical solutions to novel 3D models, Geosciences, 10, 77, https://doi.org/10.3390/geosciences10020077, 2020. a, b
ESIOpenCFD team: ESI OpenCFD Release OpenFOAM® v2312, https://www.openfoam.com/news/mainnews/openfoamv2312, last access: 18 August 2024. a
Ferziger, J. H. and Peric, M.: Computational methods for fluid dynamics, Springer, 3rd Edn., ISBN 3540420746, 2002. a, b, c, d
Fischer, J.T., Kowalski, J., and Pudasaini, S. P.: Topographic curvature effects in applied avalanche modeling, Cold Reg. Sci. Technol., 74, 21–30, https://doi.org/10.1016/j.coldregions.2012.01.005, 2012. a
Fischer, J.T., Kofler, A., Wolfgang, F., Granig, M., and Kleemayr, K.: Multivariate parameter optimization for computational snow avalanche simulation in 3d terrain, J. Glaciol., 61, 875–888, https://doi.org/10.3189/2015JoG14J168, 2015. a, b, c, d, e
Forterre, Y. and Pouliquen, O.: Flows of Dense Granular Media, Annual Rev. Fluid Mech., 40, 1–24, https://doi.org/10.1146/annurev.fluid.40.111406.102142, 2008. a, b, c
Garcés, A., González, Á., Tamburrino, A., and Montserrat, S.: faDebrisFOAM validation using field data surveyed in Crucecita (Chile) alluvial fan for the event of 13th May 2017, in: E3S Web of Conferences, 415, p. 02007, https://doi.org/10.1051/e3sconf/202341502007, 2023. a
George, D. L., Iverson, R. M., and Cannon, C. M.: New methodology for computing tsunami generation by subaerial landslides: Application to the 2015 Tyndall Glacier landslide, Alaska, Geophys. Res. Lett., 44, 7276–7284, https://doi.org/10.1002/2017GL074341, 2017. a
Hagemeier, T., Hartmann, M., and Thévenin, D.: Practice of vehicle soiling investigations: A review, Int. J. Multiphas. Flow, 37, 860–875, https://doi.org/10.1016/j.ijmultiphaseflow.2011.05.002, 2011. a
Heerema, C. J., Talling, P. J., Cartigny, M. J., Paull, C. K., Bailey, L., Simmons, S. M., Parsons, D. R., Clare, M. A., Gwiazda, R., Lundsten, E., Anderson, K., Maier, K. L., Xu, J. P., Sumner, E. J., Rosenberger, K., Gales, J., McGann, M., Carter, L., and Pope, E.: What determines the downstream evolution of turbidity currents?, Earth Planet. Sc. Lett., 532, 116023, https://doi.org/10.1016/j.epsl.2019.116023, 2020. a
Hergarten, S. and Robl, J.: Modelling rapid mass movements using the shallow water equations in Cartesian coordinates, Nat. Hazards Earth Syst. Sci., 15, 671–685, https://doi.org/10.5194/nhess156712015, 2015. a
Huber, A., Kofler, A., Rauter, M., Fischer, J.T., and Adams, M. S.: Simulation of dense snow avalanches with open source software, in: Proceedings of the International Snow Science Workshop, Innsbruck, Austria, https://arc.lib.montana.edu/snowscience/item/2643 (last access: 18 August 2024), 2018. a
Issler, D.: Modelling of snow entrainment and deposition in powdersnow avalanches, Ann. Glaciol., 26, 253–258, https://doi.org/10.3189/1998AoG261253258, 1998. a, b, c
Issler, D., Jenkins, J. T., and McElwaine, J. N.: Comments on avalanche flow models based on the concept of random kinetic energy, J. Glaciol., 64, 148–164, https://doi.org/10.1017/jog.2017.62, 2018. a, b, c
Iverson, R. M. and George, D. L.: A depthaveraged debrisflow model that includes the effects of evolving dilatancy. I. Physical basis, Proc. Roy. Soc. Lond. A, 470, 2170, https://doi.org/10.1098/rspa.2013.0819, 2014. a, b
Jasak, H.: Error analysis and estimation for the finite volume method with applications to fluid flows, Ph.D. thesis, Imperial College, University of London, https://www.researchgate.net/publication/230605842_Error_Analysis_and_Estimation_for_the_Finite_Volume_Method_With_Applications_to_Fluid_Flows (last access: 18 August 2024), 1996. a, b
Jóhannesson, T., Gauer, P., Issler, P., and Lied, K.: The design of avalanche protection dams. Recent practical and theoretical developments, No. EUR 23339 in Climate Change and Natural Hazard Research Series 2, ISBN 9789279088858, 2009. a, b, c
Juretić, F.: cfMesh User Guide, Creative Fields, Zagreb, https://cfmesh.com/wpcontent/uploads/2015/09/User_GuidecfMesh_v1.1.pdf (last access: 18 August 2024), 2015. a
Köhler, A., McElwaine, J., and Sovilla, B.: GEODAR data and the flow regimes of snow avalanches, J. Geophys. Res.Earth, 123, 1272–1294, https://doi.org/10.1002/2017JF004375, 2018. a
Kowalski, J. and McElwaine, J. N.: Shallow twocomponent gravitydriven flows with vertical variation, J. Fluid Mech., 714, 434–462, https://doi.org/10.1017/jfm.2012.489, 2013. a, b
Kowalski, J. and Torrilhon, M.: Moment Approximations and Model Cascades for Shallow Flow, Commun. Comput. Phys., 25, 669–702, https://doi.org/10.4208/cicp.OA20170263, 2019. a
LeVeque, R. J.: Finite volume methods for hyperbolic problems, vol. 31, Cambridge University Press, 2002. a
Løvholt, F., Pedersen, G., Harbitz, C. B., Glimsdal, S., and Kim, J.: On the characteristics of landslide tsunamis, Philos. T. Roy. Soc. A, 373, 20140376, https://doi.org/10.1098/rsta.2014.0376, 2015. a
Lucas, A., Mangeney, A., and Ampuero, J. P.: Frictional velocityweakening in landslides on Earth and on other planetary bodies, Nat. Commun., 5, 1–9, https://doi.org/10.1038/ncomms4417, 2014. a
Mergili, M., Fischer, J.T., Krenn, J., and Pudasaini, S. P.: r.avaflow v1, an advanced opensource computational framework for the propagation and interaction of twophase mass flows, Geosci. Model Dev., 10, 553–569, https://doi.org/10.5194/gmd105532017, 2017. a
Moukalled, F., Mangani, L., and Darwish, M.: The finite volume method in computational fluid dynamics, Springer, https://doi.org/10.1007/9783319168746, 2016. a, b
Oesterle, F.: EiskarAvalanche event January 2019 (german), Technischer Bericht/Nachrechung, Wildbach und Lawinenverbauung – Fachzentrum Geologie und Lawinen, 2019. a, b, c, d
Parker, G., Fukushima, Y., and Pantin, H. M.: Selfaccelerating turbidity currents, J. Fluid Mech., 171, 145–181, https://doi.org/10.1017/S0022112086001404, 1986. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac
Pitman, E. B., Nichita, C. C., Patra, A., Bauer, A., Sheridan, M., and Bursik, M.: Computing granular avalanches and landslides, Phys. Fluids, 15, 3638–3646, https://doi.org/10.1063/1.1614253, 2003. a
Pudasaini, S. P.: A general twophase debris flow model, J. Geophys. Res.Earth, 117, F03010, https://doi.org/10.1029/2011JF002186, 2012. a
Pudasaini, S. P. and Hutter, K.: Avalanche Dynamics: Dynamics of Rapid Flows of Dense Granular Avalanches, Springer, https://doi.org/10.1007/9783540326878, 2007. a
Pudasaini, S. P., Wang, Y., and Hutter, K.: Rapid motions of freesurface avalanches down curved and twisted channels and their numerical simulation, Philos. T. Roy. Soc. Lond. A, 363, 1551–1571, https://doi.org/10.1098/rsta.2005.1595, 2005. a
Rastello, M., Rastello, F., Bellot, H., Ousset, F., Dufour, F., and Meier, L.: Size of snow particles in a powdersnow avalanche, J. Glaciol., 57, 151–156, https://doi.org/10.3189/002214311795306637, 2011. a
Rauter, M.: The compressible granular collapse in a fluid as a continuum: validity of a NavierStokes model with μ(J),ϕ(J)rheology, J. Fluid Mech., 915, A87, https://doi.org/10.1017/jfm.2021.107, 2021. a, b, c, d, e, f, g
Rauter, M.: OpenFOAM Avalanche Module Repository, https://develop.openfoam.com/Community/avalanche, last access: 18 August 2024. a
Rauter, M. and Köhler, A.: Constraints on entrainment and deposition models in avalanche simulations from highresolution radar data, Geosciences, 10, 9, https://doi.org/10.3390/geosciences10010009, 2020. a, b, c
Rauter, M. and Tuković, Ž.: A finite area scheme for shallow granular flows on threedimensional surfaces, Comput. Fluids, 166, 184–199, https://doi.org/10.1016/j.compfluid.2018.02.017, 2018. a, b, c, d, e, f, g
Rauter, M., Fischer, J.T., Fellin, W., and Kofler, A.: Snow avalanche friction relation based on extended kinetic theory, Nat. Hazards Earth Syst. Sci., 16, 2325–2345, https://doi.org/10.5194/nhess1623252016, 2016. a, b, c, d
Rauter, M., Kofler, A., Huber, A., and Fellin, W.: faSavageHutterFOAM 1.0: depthintegrated simulation of dense snow avalanches on natural terrain with OpenFOAM, Geosci. Model Dev., 11, 2923–2939, https://doi.org/10.5194/gmd1129232018, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Rauter, M., Viroulet, S., Gylfadóttir, S. S., Fellin, W., and Løvholt, F.: Granular porous landslide tsunami modelling–the 2014 Lake Askja flank collapse, Nat. Commun., 13, 678, https://doi.org/10.1038/s41467022282967, 2022. a, b
Roache, P. J.: Quantification of uncertainty in computational fluid dynamics, Annu. Rev. Fluid Mech., 29, 123–160, https://doi.org/10.1146/annurev.fluid.29.1.123, 1997. a, b
Salm, B., Burkard, A., and Gubler, H. U.: Berechnung von Fliesslawinen: eine Anleitung für Praktiker mit Beispielen, Tech. rep., WSL Institut für Schneeund Lawinenforschung SLF, Davos, https://www.dora.lib4ri.ch/wsl/islandora/object/wsl%3A26106 (last access: 18 August 2024), 1990. a
Sampl, P. and Zwinger, T.: Avalanche simulation with SAMOS, Ann. Glaciol., 38, 393–398, https://doi.org/10.3189/172756404781814780, 2004. a, b, c, d, e, f, g
Savage, S. B. and Hutter, K.: The motion of a finite mass of granular material down a rough incline, J. Fluid Mech., 199, 177–215, https://doi.org/10.1017/S0022112089000340, 1989. a, b, c, d, e, f, g, h
Savage, S. B. and Hutter, K.: The dynamics of avalanches of granular materials from initiation to runout. Part I: Analysis, Acta Mech., 86, 201–223, https://doi.org/10.1007/BF01175958, 1991. a, b, c
Scheidegger, A. E.: On the prediction of the reach and velocity of catastrophic landslides, Rock Mech. Rock Eng., 5, 231–236, https://doi.org/10.1007/BF01301796, 1973. a
Shimizu, H. A.: Numerical Simulations of DomeCollapse Pyroclastic Density Currents Using faSavageHutterFOAM: Application to the 3 June 1991 Eruption of Unzen Volcano, Japan, J. Disaster Res., 17, 768–778, https://doi.org/10.20965/jdr.2022.p0768, 2022. a
Sovilla, B., McElwaine, J. N., and Louge, M. Y.: The structure of powder snow avalanches, C. R. Phys., 16, 97–104, https://doi.org/10.1016/j.crhy.2014.11.005, 2015. a, b, c
Steinkogler, W., Gaume, J., Löwe, H., Sovilla, B., and Lehning, M.: Granulation of snow: From tumbler experiments to discrete element simulations, J. Geophys. Res.Earth, 120, 1107–1126, https://doi.org/10.1002/2014JF003294, 2015a. a
Steinkogler, W., Sovilla, B., and Lehning, M.: Thermal energy in dry snow avalanches, The Cryosphere, 9, 1819–1830, https://doi.org/10.5194/tc918192015, 2015b. a
Tuković, Ž. and Jasak, H.: A moving mesh finite volume interface tracking method for surface tension dominated interfacial fluid flow, Comput. Fluids, 55, 70–84, https://doi.org/10.1016/j.compfluid.2011.11.003, 2012. a, b
Turnbull, B. and Bartelt, P.: Mass and momentum balance model of a mixed flowing/powder snow avalanche, Surv. Geophys., 24, 465–477, https://doi.org/10.1023/B:GEOP.0000006077.82404.84, 2003. a
Turner, J.: Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows, J. Fluid Mech., 173, 431–471, https://doi.org/10.1017/S0022112086001222, 1986. a, b, c, d
Vescovi, D., di Prisco, C., and Berzi, D.: From solid to granular gases: the steady state for granular materials, Int. J. Numer. Anal. Met., 37, 2937–2951, https://doi.org/10.1002/nag.2169, 2013. a
Viroulet, S., Baker, J., Edwards, A., Johnson, C. G., Gjaltema, C., Clavel, P., and Gray, J.: Multiple solutions for granular flow over a smooth twodimensional bump, J. Fluid Mech., 815, 77–116, https://doi.org/10.1017/jfm.2017.41, 2017. a
Voellmy, A.: Über die Zerstörungskraft von Lawinen (On the destructive forces of avalanches), Schweizerische Bauzeitung, 73, 212–217, https://doi.org/10.5169/seals61891, 1955. a
Zhao, H. and Kowalski, J.: Bayesian active learning for parameter calibration of landslide runout models, Landslides, 19, 2033–2045, https://doi.org/10.1007/s1034602201857z, 2022. a
 Abstract
 Introduction
 Foundation and framework
 Denseflow model
 Suspension flow model
 Twolayer granular flow model
 Pre and postprocessing
 Results and discussion
 Conclusions
 Appendix A: Simplified computation of centrifugal forces
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Foundation and framework
 Denseflow model
 Suspension flow model
 Twolayer granular flow model
 Pre and postprocessing
 Results and discussion
 Conclusions
 Appendix A: Simplified computation of centrifugal forces
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement