faSavageHutterFOAM 1.0: depth-integrated simulation of dense snow avalanches on natural terrain with OpenFOAM

Numerical models for dense snow avalanches have become central to hazard zone mapping and mitigation. Several commercial and free applications, which are used on a regular basis, implement such models. In this study we present a tool based on the open-source toolkit OpenFOAM^® as an alternative to the established solutions. The proposed tool implements a depth-integrated shallow flow model in accordance with current practice. The solver combines advantages of the extensive OpenFOAM infrastructure with popular models from the avalanche community. OpenFOAM allows assembling custom physical models with built-in primitives and implements the numerical solution at a high level. OpenFOAM supports an extendable solver structure, making the tool well-suited for future developments and rapid prototyping. We introduce the basic solver, implementing an incompressible, single-phase model for natural terrain, including entrainment. The respective workflow, consisting of meshing, pre-processing, numerical solution and post-processing, is presented. We demonstrate data transfer from and to a geographic information system (GIS) to allow a simple application in practice. The tool chain is based entirely on open-source applications and libraries and can be easily customised and extended. Simulation results for a well-documented avalanche event are presented and compared to previous numerical studies and historical data.

1 Introduction

Numerical avalanche modelling has become an important and well-accepted ingredient in hazard zone mapping. All popular tools rely on depth-integrated flow models (Pudasaini and Hutter, 2007) and only a few academic exceptions are known (Domnik et al., 2013; Kröner, 2013; von Boetticher et al., 2016, 2017; Barker and Gray, 2017). Depth-integrated flow models, widely known as shallow water equations, have a long tradition in hydraulic modelling (e.g. Vreugdenhil, 1994), dating back to Barré de Saint-Venant (1871). This approach is commonly applied in academia and in practice because it reduces the computational effort to a level at which physical simulations of realistic flows are feasible. The first application to gravitational mass flows is attributed to Grigorian et al. (1967) and the first formal derivation and analysis of the underlying model to Savage and Hutter (1989, 1991). Since then, the mechanical model has been continuously improved and extended to, for example, simple, two-dimensional surfaces (Greve et al., 1994), complex, shallow surfaces (Gray et al., 1999), or curved and twisted flow paths (Pudasaini et al., 2005a, b). Finally, respective models have been adapted to natural, i.e. arbitrary but mildly curved, terrain making simulations of real case avalanches possible. The limitation to mildly curved terrain requires the flow thickness to be small in relation to the curvature radius of the surface. Denlinger and Iverson (2004) proposed a model embedded in an ordinary Cartesian coordinate system as an alternative to the complex curvilinear coordinate system used by Savage and Hutter (1989, 1991). Bouchut and Westdickenberg (2004), Hergarten and Robl (2015), and recently Rauter and Tuković (2018) follow a similar approach. Christen et al. (2010) apply a non-orthogonal local coordinate system (Fischer et al., 2012) but without incorporating the respective correction terms (Hergarten and Robl, 2015). A Lagrangian solution, which has some advantages for natural terrain, has been presented by Hungr (1995) and later on by Sampl and Zwinger (2004) and Sampl and Granig (2009).

Beside improvement of the underlying mechanical model, various physical processes have been added to governing equations, such as multiple phases (e.g. Pudasaini, 2012; Kowalski and McElwaine, 2013; Iverson and George, 2016), entrainment (e.g. Issler, 2014), improved basal friction relations (e.g. Voellmy, 1955; Norem et al., 1987; Pouliquen and Forterre, 2002; Bartelt et al., 2006; Issler and Gauer, 2008; Baker et al., 2016; Rauter et al., 2016), (for a review, see Ancey, 2007), compressibility (e.g. Iverson and George, 2014; Bartelt et al., 2016) or thermodynamic processes Vera Valero et al. (e.g. 2015).

In this work, we strictly distinguish between a mechanical model and process models. The mechanical model consists of basic conservation equations and their reformulation, e.g. in terms of depth integration. Process models, on the other hand, describe the closure of governing equations, for example with constitutive models. The combination of the mechanical model and all closures is called flow model or physical model throughout this work.

There are several numerical methods to solve the respective mathematical equations. Basically, most methods can be classified as finite-difference methods (e.g. Wang et al., 2004), finite-element methods (e.g. Hanert et al., 2005), finite-volume methods (e.g. Christen et al., 2010) or as Lagrangian particle methods (e.g. Sampl and Granig, 2009). Specialised differencing schemes (e.g. upwind, TVD, NVD) prevent oscillations (e.g. Jasak et al., 1999).

Shallow granular flow models have been carefully validated over the last few decades. This includes back-calculations of small-scale experiments (for a review, see Pudasaini and Hutter, 2007), large-scale experiments (e.g. Christen et al., 2010), historic snow avalanches (e.g. Fischer et al., 2015) and rock avalanches (e.g. Mergili et al., 2017). Shallow flow models have various weaknesses, such as the limitation to mildly curved terrain or the missing resolution in surface-normal direction. However, they have proven to be a good trade-off between accuracy and computing time and thus useful for many applications.

Shallow flow models gained popularity through commercial software packages: DAN (Hungr, 1995), SamosAT (Sampl and Zwinger, 2004), FLATModel (Medina et al., 2008) and RAMMS (Christen et al., 2010) implement such models and are used regularly in practice. Open-source alternatives include TITAN2D (Pitman et al., 2003; Patra et al., 2005), r.avaflow (Mergili et al., 2012, 2017) and an extension to the CFD toolkit (computation fluid dynamics) GERRIS (Hergarten and Robl, 2015). From an academic viewpoint, open-source applications have various advantages over their commercial counterparts; for example, users can view and modify the source code to gain a better understanding of the software and adapt the flow model without re-implementing basic models and numerical methods from scratch.

Geographic information systems (GISs) are commonly applied in hazard zone mapping. Therefore numerical simulation tools are usually incorporated or linked to these systems to streamline the respective workflow. GIS allows user-friendly data input, post-processing and the production of publication-quality maps.

Recently, Rauter and Tuković (2018) proposed a shallow granular flow model, expressed in terms of surface partial differential equations (Deckelnick et al., 2005; Tuković and Jasak, 2012) and presented an open-source implementation based on the CFD toolkit OpenFOAM^® (OpenCFD Ltd., 2004). The underlying mechanical model is widely similar to the classic Savage and Hutter (1989, 1991) model and its derivations.

One particular advantage of an OpenFOAM solver is the well-designed, object-oriented source code. This makes the code cleaner than comparable solutions as it hides implementation details, such as numerical schemes, input/output or inter-process communication, behind well-defined interfaces. The top-level solver mimics the tensorial notation of partial differential equations, and specific implementations of, for example, interpolation schemes, are exchangeable without changing the top-level source code. This enables the separation of physical models and numerical solution, which allows a streamlined interdisciplinary development process. Process models, e.g. of entrainment and basal friction, can be incorporated similarly, keeping the source code clean and easy to extend.

The OpenFOAM solver, presented here, implements an incompressible single-phase model including various basal friction and entrainment closures. The solver is called faSavageHutterFoam, indicating that the underlying mechanical model is similar to the one of Savage and Hutter (1989, 1991) but with exchangeable closure models. This model is, to some extent, suitable for dense snow avalanches and constitutes the baseline for complex flow models, as employed by, for example, Bartelt et al. (2016) or Mergili et al. (2017). Moreover, the underlying method has been developed to simplify coupling with three-dimensional ambient flows (Tuković and Jasak, 2012; Marschall et al., 2014; Dieter-Kissling et al., 2015a, b; Pesci et al., 2015), which enables the development of models for mixed snow avalanches (e.g. Sampl and Zwinger, 2004) and turbidity currents (e.g. Huang et al., 2005).

The purpose of this paper is to present the capability of the new OpenFOAM solver and the Rauter and Tuković (2018) model. The solver is evaluated and validated for snow avalanches on natural terrain. We present the basic flow model, as well as methods and tools to incorporate natural terrain and GIS data in OpenFOAM simulations. Also, the export of OpenFOAM results to a GIS for post-processing and visualisation is demonstrated. Results for a well-documented avalanche event are presented and compared to historical records and results of SamosAT. All underlying source code (except SamosAT) and data are available free of charge to encourage reproduction, improvement and cross-validation.

2 Method

2.1 Flow model

Historically, shallow granular flow models have been set up in surface-aligned, curvilinear coordinates, leading to a two-dimensional system of partial differential equations (e.g. Savage and Hutter, 1989, 1991). Rauter and Tuković (2018) follow a different approach (see also Denlinger and Iverson, 2004; Bouchut and Westdickenberg, 2004; Hergarten and Robl, 2015) and formulate the mechanical model in terms of surface partial differential equations (SPDEs; e.g. Deckelnick et al., 2005). Respective SPDEs are defined on a surface Γ, embedded in three-dimensional space, which represents the mountain topography. This approach, popular in the thin liquid-film community (e.g. Craster and Matar, 2009), avoids transformations into the surface-aligned coordinate system and thus complex metric tensors. Considering the relative shallowness of the avalanche, it can be treated as a thin layer flowing along the mountain surface. The governing equations describe the motion of the avalanche in three-dimensional space along this surface. Consequently, velocity is a three-dimensional vector field and contains all information on flow direction and respective effects, such as centrifugal forces. Resulting SPDEs can be solved with various methods, e.g. the finite-element method (e.g. Olshanskii et al., 2009) or the finite-area method, a modified finite-volume method (Tuković and Jasak, 2012).

Figure 1 Definition of velocity u, flow thickness h and basal pressure p_b on a control volume. A hydrostatic and linear pressure distribution is assumed. The shape of the velocity profile is commonly ignored in governing equations (Baker et al., 2016). Flow thickness h is measured normal to the basal surface Γ. The curvature radius of the surface Γ is assumed to be much bigger than the flow thickness.

2.1.1 Mechanical model

A basic shallow granular flow model can be written in terms of surface partial differential equations as¹

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial h}{\partial t}}+\mathbf{\nabla }\cdot \left(h\stackrel{\mathrm{‾}}{\mathit{u}}\right)={\displaystyle \frac{\dot{q}}{\mathit{\rho }}},\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial \left(h\stackrel{\mathrm{‾}}{\mathit{u}}\right)}{\partial t}}+{\mathbf{\nabla }}_{\mathrm{s}}\cdot \left(h\stackrel{\mathrm{‾}}{\mathit{u}}\stackrel{\mathrm{‾}}{\mathit{u}}\right)=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathit{\tau }}_{\mathrm{b}}+h{\mathit{g}}_{\mathrm{s}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\rho }}}{\mathbf{\nabla }}_{\mathrm{s}}\left(h{p}_{\mathrm{b}}\right),\text{(3)}& {\displaystyle }& {\displaystyle }{\mathbf{\nabla }}_{\mathrm{n}}\cdot \left(h\stackrel{\mathrm{‾}}{\mathit{u}}\stackrel{\mathrm{‾}}{\mathit{u}}\right)=h{\mathit{g}}_{\mathrm{n}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\rho }}}{\mathbf{\nabla }}_{\mathrm{n}}\left(h{p}_{\mathrm{b}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathit{n}}_{\mathrm{b}}{p}_{\mathrm{b}}.\end{array}$$

Equations (1) to (3) are equivalent to a Savage–Hutter-like system, consistently extended to complex but mildly curved terrain and entrainment. The notation as SPDE makes extension to complex terrain straightforward and implementation into SPDE environments, e.g. OpenFOAM, possible. A formal derivation is given by Rauter and Tuković (2018). Here, we aim to deliver a short and descriptive introduction.

Equation (1) represents the depth-integrated continuity equation, Eq. (2) the surface-tangential momentum conservation equation and Eq. (3) its surface-normal counterpart, defined at all points x_b on the surface Γ⊂ℝ³, representing the mountain surface. The time is denoted as t. The unknown fields are the surface-normal flow thickness h(x_b) (see Fig. 1), the depth-averaged flow velocity $\stackrel{\mathrm{‾}}{\mathit{u}}\left({\mathit{x}}_{\mathrm{b}}\right)\in {\mathbb{R}}^{\mathrm{3}}$, defined as

$$\begin{array}{}\text{(4)}& \stackrel{\mathrm{‾}}{\mathit{u}}\left({\mathit{x}}_{\mathrm{b}}\right)={\displaystyle \frac{\mathrm{1}}{h\left({\mathit{x}}_{\mathrm{b}}\right)}}\underset{\mathrm{0}}{\overset{h\left({\mathit{x}}_{\mathrm{b}}\right)}{\int }}\mathit{u}({\mathit{x}}_{\mathrm{b}}-{\mathit{n}}_{\mathrm{b}}{z}^{\prime })\mathrm{d}{z}^{\prime },\end{array}$$

and the basal pressure p_b(x_b). The density ρ is assumed to be constant. Note that the earth pressure theory (e.g. Savage and Hutter, 1989, 1991) has been replaced with the hydrostatic pressure assumption, as in most practical applications (e.g. Christen et al., 2010). Moreover, Eqs. (1) to (3) are written in conservative form. Therefore, there is no entrainment term in Eq. (2), which would show up in a non-conservative formulation. The first terms in Eqs. (1) and (2) represent the temporal derivative, i.e. the local change in mass and momentum, respectively. The second terms in Eqs. (1) and (2) are the respective advection terms. The right-hand side of Eq. (1) represents mass growth due to entrainment. The first, second and third terms on the right-hand side of Eq. (2) represent surface-tangential components of basal friction, gravitational acceleration and lateral pressure gradient, respectively. The surface-normal components of these terms appear in the surface-normal momentum conservation equation (Eq. 3). This equation is used to calculate the basal pressure, represented by the last term.

In the framework of SPDEs, the normal vector field n_b(x_b)∈ℝ³ of the surface Γ is sufficient to describe all major curvature effects. This is realised by calculating all contributions to conservation equations in the global coordinate system and projecting results onto the surface and the surface-normal vector. These projections are explained in detail in Appendix A. Surface-tangential and normal components contribute to local acceleration and basal pressure, respectively. This follows from the assumption that movement is constrained in surface-normal direction, which is enforced by a mechanical force, namely the basal pressure. The gravitational acceleration g, for example, is split into a surface-tangential component,

$$\begin{array}{}\text{(5)}& {\mathit{g}}_{\mathrm{s}}=\left(\mathbf{I}-{\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathit{g},\end{array}$$

and a surface-normal component,

$$\begin{array}{}\text{(6)}& {\mathit{g}}_{\mathrm{n}}=\left({\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathit{g}.\end{array}$$

The gradient operator ∇ denotes the three-dimensional derivative along the surface (Deckelnick et al., 2005). If the responding result is a three-dimensional vector field (e.g. gradient of a scalar field or divergence of a tensor field), it can be split, similar to the gravitational acceleration, into a surface-tangential component,

$$\begin{array}{}\text{(7)}& {\mathbf{\nabla }}_{\mathrm{s}}=\left(\mathbf{I}-{\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathbf{\nabla },\end{array}$$

and a surface-normal component,

$$\begin{array}{}\text{(8)}& {\mathbf{\nabla }}_{\mathrm{n}}=\left({\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathbf{\nabla }.\end{array}$$

For simply curved surfaces, the given relation matches the model of Greve et al. (1994), as shown by Rauter and Tuković (2018).

2.1.2 Process models

There are various user-selectable models, describing basal friction τ_b(x_b) and entrainment rate $\dot{q}\left({\mathit{x}}_{\mathrm{b}}\right)$, to close the system of equations. To reassemble the traditional model Christen et al. (often called Voellmy or Voellmy–Salm model; 2010), as applied, for example, by Fischer et al. (2015), the basal friction is described following Voellmy (1955),

$$\begin{array}{}\text{(9)}& {\mathit{\tau }}_{\mathrm{b}}=\mathit{\mu }{p}_{\mathrm{b}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{u}}}{\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|+u_{\mathrm{0}}}}+{\displaystyle \frac{\mathit{\rho }g}{\mathit{\xi }}}\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|\stackrel{\mathrm{‾}}{\mathit{u}}.\end{array}$$

Therein, μ and ξ are constant parameters, although they may depend on avalanche size and surface roughness (Salm et al., 1990) or flow regime (Köhler et al., 2016). The small value u₀ (${\mathrm{10}}^{-\mathrm{7}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ here) avoids divisions by zero and regularises the relation near standstill, where the original function is discontinuous. This regularisation, combined with the employed time integration scheme (implicit three-level second-order; Ferziger and Peric, 2002), leads to well-defined behaviour in the runout zone, where the velocity is nearly zero (Rauter and Tuković, 2018). This allows the avalanche to reach very low velocities in the runout zone, which are lower than the tolerance of the solver and thus virtually zero. For characteristic avalanche velocities, i.e. $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|>\mathrm{100}{u}_{\mathrm{0}}$, this value has no relevant effect on the dynamic behaviour. Previously, this issue has been addressed with operator splitting and explicit stress reduction (e.g. Mangeney-Castelnau et al., 2003; Zhai et al., 2014; Mergili et al., 2017), which is not required in the proposed scheme.

The entrainment rate is calculated, based on an empirical erosive entrainment model, as

$$\begin{array}{}\text{(10)}& \dot{q}=\left\{\begin{array}{cc}{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}\cdot \stackrel{\mathrm{‾}}{\mathit{u}}}{e_{\mathrm{b}}}}& \text{for}{h}_{\mathrm{msc}}>\mathrm{0},\\ \mathrm{0}& \text{for}{h}_{\mathrm{msc}}=\mathrm{0},\end{array}\right.\end{array}$$

where e_b is the specific erosion energy (Fischer et al., 2015). Entrainment is restricted by the available mountain snow cover thickness h_msc. The initial mountain snow cover thickness is calculated following Fischer et al. (2015), using a linear approach,

$$\begin{array}{}\text{(11)}& h_{\mathrm{msc}}\left(z\right)=\left(H_{\mathrm{msc}}\left(z_{\mathrm{0}}\right)+{\displaystyle \frac{\partial H_{\mathrm{msc}}}{\partial z}}\left(z-z_{\mathrm{0}}\right)\right)\cos \left(\mathit{\zeta }\right),\end{array}$$

where z is the surface elevation (corresponding to the vertical coordinate in the numerical model) and z₀ is the elevation of a reference station, which has to be provided by the user, alongside with the base value H_msc(z₀) and the growth rate $\frac{\partial H_{\mathrm{msc}}}{\partial z}$. ζ is the angle between the gravitational acceleration and the surface-normal vector. Its further evolution is described by the conservation equation

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial h_{\mathrm{msc}}}{\partial t}}=-{\displaystyle \frac{\dot{q}}{\mathit{\rho }}}.\end{array}$$

Undershoots, i.e. h_msc<0, are prevented with a regularisation similar to Eq. (9). This can be realised by multiplying the entrainment rate $\dot{q}$ with $\frac{h_{\mathrm{msc}}}{h_{\mathrm{msc}}+h_{\mathrm{0}}}$, where h₀ is a small value, similar to u₀.

2.1.3 Numerical solution

The governing equations are solved with an implicit, conservative, finite-area method (Rauter and Tuković, 2018), using the respective OpenFOAM library (Tuković and Jasak, 2012). The finite-area method is similar to the well-known finite-volume method Jasak (e.g. 1996) but with appropriate differential operators for SPDEs: Eqs. (7) and (8). We apply first- (upwind scheme) and second-order accurate spatial differencing schemes. First-order schemes converge more slowly in terms of mesh refinement due to their high numerical diffusivity. However, they effectively prevent oscillations and increase the stability of the solver. Oscillations in second-order accurate simulations are prevented with a normalised variable diagram (NVD) scheme for unstructured meshes, known as the Gamma scheme (Jasak et al., 1999). NVD schemes blend upwind and a higher-order scheme to combine advantages of both methods.

As mentioned before, OpenFOAM utilises capabilities of C++ to make top-level source code appear similar to the tensor notation of partial differential equations. The conservation equation (Eq. 1), for example, can be solved with the following lines of code using OpenFOAM:

faScalarMatrix hEqn (
    fam::ddt(h)
  + fam::div(phis, h)
 ==
    dqdt/rho
); hEqn.solve();

phis is the velocity edge field (see Rauter and Tuković, 2018, for details), and dqdt is the source term incorporating entrainment. Momentum conservation equations (Eqs. 2 and 3) look similar (see Rauter and Tuković, 2018), and conservation equations for arbitrary fields (e.g. random kinetic energy, Bartelt et al., 2016) can be added with the same syntax.

2.2 Simulation evaluation

We use an established implementation of the same flow model, SamosAT (version 2017_07_05) (Sampl and Zwinger, 2004; Sampl and Granig, 2009), for comparison. The main difference between SamosAT and the presented OpenFOAM solver is the solution method. SamosAT solves similar governing equations, slightly adapted to fit into the respective framework, with smoothed-particle hydrodynamics (SPH). This approach follows a Lagrangian description, making the handling of complex terrain simpler (Sampl and Zwinger, 2004). Therefore, SamosAT provides an excellent reference to validate avalanche models for complex terrain. The second term on the right-hand side of Eq. (3) was deactivated in OpenFOAM computations to reassemble the mechanical model as implemented in SamosAT. This term is usually small and can be safely neglected (Rauter and Tuković, 2018). However, it is shown in equations to preserve the similarity between Eqs. (2) and (3).

We compare simulations using the 1 kPa isoline of the dynamic peak pressure, defined as

$$\begin{array}{}\text{(13)}& p_{\mathrm{dyn}}\left({\mathit{x}}_{\mathrm{b}}\right)=\underset{t}{max}\left(\mathit{\rho }\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right({\mathit{x}}_{\mathrm{b}},t){|}^{\mathrm{2}}\right).\end{array}$$

Definitions of hazard zones are based on this threshold in many European countries (Jóhannesson et al., 2009) and are therefore often used for the evaluation of relevant models (e.g. Fischer et al., 2015; Rauter et al., 2016).

In addition to the comparison with a reference implementation, we present a comparison with historical records from a catastrophic event. A common method to document avalanches is the delineation of deposition. This information is also available for the presented case study. Deposition processes are not explicitly included in the flow model due to depth integration. However, the general form and size of the deposition should be reproduced by the model to be useful for hazard zone mapping. This is problematic in some implementations, e.g. SamosAT, due to missing regularisation of the friction term, but it is possible with the proposed method.

We apply model parameters (μ, ξ, e_b) optimised for SamosAT (Fischer et al., 2015), and the comparison is conducted on a qualitative level.

Finally, we evaluate OpenFOAM simulations with regard to convergence during mesh refinement to give a quantitative estimation of numerical uncertainties as recommended by Roache (1997). The numerical solution should converge to the unknown analytical solution with increasing grid resolution, and the numerical uncertainty should decay with the order of the applied method. Richardson extrapolation allows us to estimate the numerical uncertainty, using results of three different meshes. This way, the expected convergence can be verified and the numerical uncertainty quantified.

2.3 Simulation set-up

The precondition to conduct simulations in OpenFOAM is a mesh, describing the geometry of the problem. For SPDEs, e.g. shallow flow models, a surface mesh, matching the slope topography, is sufficient and no volume mesh is required. In practice, however, three-dimensional meshing tools can be used to create a volume mesh, the boundary of which can be used as surface mesh.

Figure 2 Simulation set-up and tool chain. The tool chain consists exclusively of open-source applications. Individual applications and process models can be replaced with custom ones. Parameters for python scripts are provided via a command line interface. Parameters for OpenFOAM applications are provided through OpenFOAM dictionaries. Domain decomposition and reconstruction, which is handled by separate applications, is not shown. OpenFOAM reads initial conditions from the folder “0” and writes results to folders named after the corresponding time step (“1”, “2”, etc.). Details on OpenFOAM formats can be found in OpenCFD Ltd. (2004).

Topography is usually available as a digital elevation model (DEM) in GIS formats, yielding elevation on a regular two-dimensional grid. The relevant part of the topography is re-sampled with cubic splines, triangulated and stored as an STL file (e.g. Kai et al., 1997) to prepare it for meshing. We chose the meshing application cfMesh (Juretić, 2015) because of its good integration in OpenFOAM and its clean boundary meshes. cfMesh requires a closed triangulated surface to create a volume mesh. This is the case for all general purpose meshing tools, and cfMesh can be replaced easily in our tool chain, for example with Netgen (Schöberl, 1997) (see Rauter and Tuković, 2018, for an application). Various other meshing tools can be applied and OpenFOAM provides a large range of mesh conversion tools. The closed surface can be assembled from a triangulation of the mountain surface, sidewalls and the respective top boundary. The resulting surface and volume mesh are presented in Fig. 3b and c. Refinement near the mountain surface reduces the amount of required volume cells, while keeping the number of surface cells high. The resulting mesh is also valid for three-dimensional simulations with, for example, Navier–Stokes equations, as conducted by, for example, Sampl and Zwinger (2004), Dutykh et al. (2011), Kröner (2013), von Boetticher et al. (2016), von Boetticher et al. (2017) and Huang et al. (2005). The boundary mesh, describing the mountain surface, is shown in Fig. 3d. The shallow flow model is solved on this surface mesh. We used polygonal-dominated (or volumetric polyhedral-dominated) meshes for simulations for stability and accuracy reasons (Juretić, 2005). Triangular (or volumetric tetrahedral) meshes have been evaluated as well. However, second-order accurate simulations on triangular meshes failed, while first-order accurate simulations are virtually identical to the respective simulations on polygonal-dominated meshes.

The release area, acting as an initial condition, is provided as a polygon in an ESRI shape file format (ESRI, 1998). To find all surface cells within the given polygon, the Hormann and Agathos (2001) algorithm as implemented in OpenFOAM is applied. The mountain snow cover h_msc of the corresponding cells is then transferred to the flow thickness h to create a suitable initial condition. The release area for our case study, taken from Fischer et al. (2015), is shown in Fig. 7a as a polygon and as a set of surface cells in Fig. 3d.

Figure 3 Meshing tool chain: the terrain data are usually available as raster data (a). Triangulation of the relevant area and adding walls and a top boundary yields a closed triangulated surface (b; sharp edges are highlighted black). This surface can be processed by most meshing tools; here, we apply cfMesh to get a polyhedral-dominated finite-volume mesh (c). The bottom boundary surface of the finite-volume mesh builds the foundation for the finite-area mesh used for simulations (d). Note that we show a very coarse mesh for the sake of visibility of the edges. Terrain data: Amt der Tiroler Landesregierung (AdTLR). EPSG coordinate reference system: 31254.

Figure 4 Time series of an OpenFOAM simulation with mean cell size Δ=7.45 m and first-order interpolations in ParaView. The colour scale represents flow thickness, which is clipped at 0.5 m. Terrain data: AdTLR.

The solver reads the surface mesh and initial conditions, as well as physical models, numerical schemes and constants to initialise the simulation (see Fig. 2). The respective entries can be found in the designated locations, according to the usual practice in OpenFOAM (OpenCFD Ltd., 2004). The solver can run on multiple processors using domain decomposition (Weller et al., 1998) and message passing interface (MPI).

User-defined friction and entrainment models can be loaded at run-time, meaning that the user does not have to recompile the solver to add a custom friction or entrainment model. The same is the case for general purpose functions which are triggered at the end of every time step. Here, we used this interface to calculate and record the dynamic peak pressure at run-time, without the necessity to save multiple time steps or to change solver source code. Similar functions can be used to check mass, momentum or energy conservation, record specific data (e.g. time line at a certain point), or to manipulate fields during run-time, e.g. to trigger secondary slabs.

Simulation results are written to hard disk in the usual OpenFOAM file format (OpenCFD Ltd., 2004) for post-processing, evaluation and simulation restart. The simulation set-up, all involved applications, and all intermediate and final files are presented in Fig. 2. The tool chain is modularly assembled from various open-source applications. Single modules, such as mesher, solver or friction model, can be replaced easily.

2.4 Post-processing and visualisation

Post-processing and visualisation of OpenFOAM simulations is commonly performed using ParaView^® (Ahrens et al., 2005; Ayachit, 2015) (see Figs. 3, 4 and 5). ParaView is an open-source data analysis and visualisation application. It can read and visualise OpenFOAM files, and they can be used for further operations, such as the calculation of contour lines. To integrate GIS applications in post-processing, results can be exported to common GIS file formats. Contour lines can be exported to ESRI shape file format with a custom python extension based on the library pyshp (Figs. 6c, 7 and 9). Alternatively, individual cells and respective field values can be exported as polygons (Fig. 6a) or points (Fig. 6b) to ESRI shape files.

Figure 5 Perspective view on the OpenFOAM simulation with mean cell size Δ=7.45 m and first-order interpolations in ParaView. The colour scale represents flow thickness, which is clipped at 0.5 m. Terrain data: AdTLR.

Figure 6 The flow thickness field h at time t=40 s for a simulation with mean cell size Δ=7.45 m; first-order interpolations. The figure shows four methods to export and analyse results in GIS: export of cells as polygons (a); export of cell centres as points (b); export of contour lines as polygons (c); remapping of the unstructured finite-area mesh to a regular raster (d). The raster has been created by converting point data to a raster file in QGIS. The resolution of the DEM is 10 m, results have been mapped to a 5 m grid. Terrain data: AdTLR.

To generate regular raster files, the unstructured OpenFOAM mesh and associated fields have to be mapped to a structured Cartesian grid (Figs. 6d and 9). These and other approaches allow an almost seamless integration into general purpose GIS applications, as shown in the following case study. Here, we utilise foam-extend-4.0 with a custom solver, python 2.7.12 with numpy 1.11.0, scipy 0.17.0 and pyshp 1.2.3 for shape file export, ParaView 5.0.1, and QGis 2.8.6.

3 Case study

In this work we focus on the Wolfsgruben avalanche. The event from 13 March 1988, when the avalanche struck inhabited areas, has been repeatedly used as benchmark for avalanche simulations, most recently by Fischer et al. (2015). We chose this example because the relevant data are freely available, making reproduction and cross-validation possible.

The mountain snow cover thickness for the specific event can be described with the parameters H_msc(z₀)=1.61 m, z₀=1289 m and $\frac{\partial H_{\mathrm{msc}}}{\partial z}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$. Physical parameters to reassemble the runout properly are μ=0.26, ξ=8650 m s⁻² and $e_{\mathrm{b}}=\mathrm{11}\mathrm{500}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$. These parameters were optimised in a previous study using SamosAT (Fischer et al., 2015).

Numerical parameters for OpenFOAM (see Rauter and Tuković, 2018) have been chosen such that they do not influence the results, while keeping the solver as stable as possible. The appropriate mesh resolution for OpenFOAM has been identified using a mesh refinement study, which is presented alongside the results. The simulation duration has been set to 150 s. This duration is sufficient to reach standstill (i.e. a velocity lower than the solver tolerance, $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|<{\mathrm{10}}^{-\mathrm{5}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) in the runout zone and thus virtually unchanging deposition. We decomposed the simulation domain into four parts for OpenFOAM and all simulations have been conducted on a Quadcore Intel Core i7-7700K @ 4.20 GHz and 32 GB DDR4 Ram @ 2.667 GHz.

SamosAT utilises a grid with 5 m resolution, and we follow recommendations in terms of appropriate particle numbers and other numerical parameters. The interpolation method has been varied between interpolation on a grid (SPH-mode 0) and interpolation on particles (SPH-mode 1) to get an insight into the numerical uncertainty.

ParaView renderings are presented in Fig. 4 for multiple time steps, showing the dynamic behaviour of the avalanche. A perspective ParaView rendering is shown in Fig. 5. The avalanche follows the narrow channel directly beneath the release area. Small portions of the avalanche overflow the left and right humps in some simulations, which can be seen in the peak dynamic pressure (Fig. 7).

The results at time step t=40 s have been exported to QGIS using various methods; see Fig. 6. Affected areas (i.e. 1 kPa isolines), as predicted by OpenFOAM and SamosAT, are shown in Fig. 7. Variations due to different interpolation schemes are shown for both implementations to give an insight into the numerical uncertainty.

The influence of the mesh resolution on the affected area is shown in Fig. 8 for the OpenFOAM solver. Respective mean cell sizes, an estimation of the numerical uncertainty following Roache (1997) and execution times (excluding time for mesh generation, which may take several minutes) are presented in Table 1. Here, the runout is defined as the length of the central avalanche path (see Fig. 8) within the affected area. The central avalanche path has been taken from Fischer et al. (2015). The mean cell size is defined as the square root of the mean cell area. For comparison, execution times for SamosAT are 98 s (SPH-mode 0) and 368 s (SPH-mode 1). One should keep in mind that SamosAT only utilises a single processor core while OpenFOAM utilises all available cores. Moreover, execution times should be seen as rough estimates because they depend on various factors, such as the number of saved time steps, debug messages and compile options.

Table 1 Mesh size, runout, error estimation and execution time for different OpenFOAM simulations. Base cell size and refinements refer to parameters of cfMesh.

Deposition (i.e. flow thickness field h in the last time step) of the OpenFOAM solution is shown in Fig. 9 alongside with the documentation.

Figure 7 Comparison of OpenFOAM first order (blue, dashed), OpenFOAM second order (blue), SamosAT SPH 0 (red, dashed) and SamosAT SPH 1 (red) in terms of 1 kPa isolines (affected area). OpenFOAM results are based on the mesh with cell size Δ=7.45 m. The documented release area (orange area) and documented deposition area (blue area) are shown for orientation. The shape and reach of the main avalanche branch are similar in all simulations; secondary branches differ to some extent. Overview (a) and focus on the runout zone (b). Terrain data: AdTLR.

Figure 8 Mesh refinement and convergence study for the OpenFOAM solver. Four mesh sizes and both interpolation schemes, first-order upwind (dashed line) and second-order Gamma (solid line) have been evaluated. The central avalanche path from Fischer et al. (2015) is shown in black. Terrain data: AdTLR.

Figure 9 Flow thickness field h at t=150 s of the second-order OpenFOAM simulation (Δ=3.75 m) and the documented deposition area. The flow thickness field in the last time step should roughly replicate the deposition. The bulge on the orographic right side of the deposition area is not matched by any simulation. However, some interesting details, such as the tail of the avalanche are represented well in OpenFOAM simulations. Map data: http://basemap.at (last access: 1 March 2018).

4 Discussion and conclusion

Results of the new OpenFOAM solver are very similar to SamosAT. Differences between SamosAT and OpenFOAM are in the range of numerical uncertainty, and differences between interpolation methods are of a comparable size. This uncertainty has to be expected; in fact, it is well known in the CFD community that numerical schemes and implementation details influence results if they are not converged to the analytical solution (e.g. Ferziger and Peric, 2002). In the case of gravitational mass flows, numerical uncertainty plays a minor role, since underlying models, parameters, terrain and snow cover data are affected by substantially higher uncertainty. This is shown by a comparison of the documented deposition with the result of an OpenFOAM simulation in Fig. 9. Although parameters have been optimised to the specific event, all simulations differ significantly from documentation. In particular, the large bulge on the orographic right side of the deposition area is not matched by any simulation. However, some details, such as the form of the tail and the position where the deposition expands, are accurately simulated by the OpenFOAM solver. Significant differences between simulation and documentation are not limited to the presented case and have been observed before, for example, by Rauter et al. (2016). We deduce that numerical errors are much smaller than the expected model error. Under these circumstances, a quantitative comparison between implementations (as, for example, by Rauter et al., 2016, for basal friction models) is not appropriate.

The refinement study shows that in the presented case, the simulated runout reduces with increasing mesh refinement (Fig. 8). Simulations on fine meshes are stopped by the first embankment, simulations on coarser grids overflow it and reach the next embankment. This is reasonable, considering the higher diffusivity and lower curvature of coarser meshes. However, this trend should not be taken for granted for other cases and a refinement study should always be conducted to get an insight into the numerical uncertainty. Results indicate that a cell size of approximately 3.75 m is required in OpenFOAM to achieve convergence with respect to practical applications. The numerical uncertainty cannot be calculated for the coarsest and finest mesh, since three simulations are required to conduct a Richardson extrapolation. It has to be noted that all simulations are based on the same DEM with a grid size of 10 m. The influence of terrain model quality (see, e.g., Bühler et al., 2011) on simulation results is not investigated.

The execution time of the OpenFOAM solver is acceptable for coarse meshes but increases with the square of the number of cells because the time step duration has to be reduced similarly to cell size. The OpenFOAM solver is noticeably slower than SamosAT, especially when considering OpenFOAM's multiprocessing capabilities. For applications where fast execution is imperative, such as parameter studies, SamosAT may be the appropriate choice. There is potential for future optimisation in OpenFOAM; in particular, the implicit time integration scheme is expensive and should be replaced with a simpler explicit one. However, the implicit solution strategy, in combination with the regularised friction relation, leads to satisfying behaviour in the runout zone. In contrast, the simple explicit solution strategy, for example from SamosAT, leads to a continuous creeping of the deposition, meaning that the final flow thickness cannot be compared with the deposition, as noted by Fischer et al. (2015) and Rauter et al. (2016).

The stability of the OpenFOAM solver is strongly influenced by mesh quality. Simulations with polygonal-dominated surface meshes showed an acceptable stability for first- and second-order interpolations. The high influence of the three-dimensional mesh on stability and its computationally expensive creation is the main drawback of the proposed method. This is, however, also a big advantage, allowing simple coupling with three-dimensional ambient flows, as conducted by Sampl and Zwinger (2004).

5 Summary and outlook

This paper shows the application of a finite-area scheme for shallow granular flows (Rauter and Tuković, 2018) to snow avalanches on natural terrain. Specific processes, such as entrainment, have been added to the basic model to replicate the traditional model as implemented in SamosAT (Fischer et al., 2015).

Various simulations with the new OpenFOAM solver have been conducted. Methods and tools to incorporate the OpenFOAM solver in GIS have been presented. These tools allow the integration of OpenFOAM in hazard mapping workflows and thus the validation of the OpenFOAM solver with a reference implementation, herein SamosAT.

The application of three-dimensional Cartesian coordinates allows simple coupling with GIS applications because no coordinate transformations are required. Unstructured meshes, on the other hand, require re-sampling to structured meshes or data transfer in the form of polygons. This incorporates an additional effort compared to simulations on structured meshes, as conducted for example by Christen et al. (2010).

The OpenFOAM solver roughly reproduces the results of SamosAT. Differences are within the expected numerical uncertainty. A comparison of numerical results to a documented event suggests that model uncertainty is substantially higher than numerical uncertainties.

The major advantage of OpenFOAM is the object-oriented open-source code, which can be extended easily. The flexible code structure allows fast application of new models to real case examples. This especially qualifies the proposed method for model development and academic purposes. Moreover, the vast majority of source code is shared within the OpenFOAM community, leading to faster development of core features and higher code quality.

The finite-area scheme allows a description in terms of surface partial differential equations (Deckelnick et al., 2005), which leads to simple and expressive governing equations. However, this comes at the cost of a complex three-dimensional surface mesh. The projection of the governing equations on a plane surface following, for example,  Bouchut and Westdickenberg (2004) may be beneficial for some applications. The three-dimensional surface mesh can also be an advantage, allowing a simple coupling with three-dimensional ambient two-phase models for powder clouds (Sampl and Zwinger, 2004). The presented meshing method, creating a finite-volume and the corresponding finite-area mesh, is viable for such simulations as well.

Future steps will incorporate the optimisation of the solver in terms of stability and execution time. Mesh generation and the integration of geographic information systems will be further streamlined. The limitation to mildly curved terrain should be eliminated, as this assumption is violated in many practical cases. We aim to implement more complex models, suitable for mixed snow avalanches (e.g. Bartelt et al., 2016; Issler et al., 2018) and debris flow (e.g. Iverson and George, 2014; Mergili et al., 2017) in the near future. Coupling of the dense flow model proposed here with three-dimensional two-phase models for the powder cloud regime (e.g. Cheng et al., 2017; Chauchat et al., 2017) is planned as well.

Appendix A: Understanding projections in surface partial differential equations

Here we briefly explain the concept of projections within the framework of surface partial differential equations. These projections are widely used in computational fluid dynamics, usually when surfaces in three-dimensional space are considered. We do not focus on mathematical formalities and this section cannot replace the formal derivation of Rauter and Tuković (2018). We want to emphasise that no surface-aligned coordinate system is required throughout the whole process, and the reader is encouraged to adhere to global Cartesian coordinates. For simplicity we present a discretised finite-area cell, which has been extruded by flow thickness h to present the flowing mass; see Fig. A1.

Figure A1 Splitting gravitational acceleration into a surface-tangential and surface-normal part with simple projections to the surface-normal vector n_b.

We begin by splitting a simple vectorial entity, the gravitational acceleration g∈ℝ³, into a surface-normal component, g_n∈ℝ³, and a surface-tangential component, g_s∈ℝ³, as shown in Fig. A1. The magnitude of the surface-normal component can be calculated using the scalar product and the surface-normal vector:

$$\begin{array}{}\text{(A1)}& \Vert {\mathit{g}}_{\mathrm{n}}\Vert ={\mathit{n}}_{\mathrm{b}}\cdot \mathit{g},\end{array}$$

which corresponds to a projection of g on n_b. The surface-normal component points in the same direction as the surface-normal vector, which allows the calculation of the vectorial surface-normal component. Rearranging of vector multiplications yields the known form

$$\begin{array}{}\text{(A2)}& {\mathit{g}}_{\mathrm{n}}={\mathit{n}}_{\mathrm{b}}\Vert {\mathit{g}}_{\mathrm{n}}\Vert ={\mathit{n}}_{\mathrm{b}}\left({\mathit{n}}_{\mathrm{b}}\cdot \mathit{g}\right)=\left({\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathit{g}.\end{array}$$

The surface-tangential component follows by subtracting the surface-normal component from total gravitational acceleration:

$$\begin{array}{}\text{(A3)}& {\mathit{g}}_{\mathrm{s}}=\mathit{g}-{\mathit{g}}_{\mathrm{n}}=\mathit{g}-\left({\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathit{g}=\left(\mathbf{I}-{\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathit{g}.\end{array}$$

Movement in surface-normal direction is constrained by the basal topography, which yields the basal pressure. Therefore, the surface-normal component g_n has to contribute to basal pressure p_b (Eq. 3), and only the surface-tangential component contributes to local acceleration $\frac{\partial h\stackrel{\mathrm{‾}}{\mathit{u}}}{\partial t}$ (Eq. 2). The total gravitational acceleration can be reconstructed by summing up both components:

$$\begin{array}{}\text{(A4)}& \mathit{g}={\mathit{g}}_{\mathrm{n}}+{\mathit{g}}_{\mathrm{s}}=\left({\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathit{g}+\left(\mathbf{I}-{\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathit{g}=\mathbf{I}\cdot \mathit{g}=\mathit{g},\end{array}$$

ensuring perfect conservation of three-dimensional momentum.

Figure A2 Splitting the divergence of a flux tensor ∇⋅M into a surface-tangential and surface-normal part with simple projections to the surface-normal vector n_b.

The same concept can be applied to fluxes through the boundary of the control volume, leading to the concept of surface partial differential operators (∇_s and ∇_n). Figure A2 shows the divergence of a tensor, ∇⋅M, which could represent convective momentum transport $\mathbf{\nabla }\cdot \left(h\stackrel{\mathrm{‾}}{\mathit{u}}\stackrel{\mathrm{‾}}{\mathit{u}}\right)$ or the lateral pressure gradient $\frac{\mathrm{1}}{\mathrm{2}\mathit{\rho }}\mathbf{\nabla }\left(p_{\mathrm{b}}h\right)=\frac{\mathrm{1}}{\mathrm{2}\mathit{\rho }}\mathbf{\nabla }\cdot \left(\mathbf{I}{p}_{\mathrm{b}}h\right)$. Using Gauss' theorem, the divergence can be reformulated in terms of the surface integral of face fluxes, which are defined as the scalar product of the flux tensor M with the normal vector on the face (Ferziger and Peric, 2002). In the discretised form, integrals are replaced with sums over faces and in the case of SPDEs, volumes collapse to surfaces, faces to edges and face fluxes to edge fluxes. For the simple case, as shown in Fig. A2, we can write

$$\begin{array}{}\text{(A5)}& \mathbf{\nabla }\cdot \mathbf{M}={\displaystyle \frac{\mathrm{1}}{S_{\mathrm{b}}}}\left({\mathit{m}}_{\mathrm{out}}-{\mathit{m}}_{\mathrm{in}}\right),\end{array}$$

with area of the cell S_b and edge fluxes m_in and m_out. For the exact formulation in terms of finite areas, the reader is refereed to Rauter and Tuković (2018). Note that ∇⋅M is a three-dimensional vector without any particular direction in relation to the basal surface. Hence, it has a surface-tangential and a surface-normal component which can be treated similarly to gravitational acceleration, yielding the surface-normal component

$$\begin{array}{ll}\text{(A6)}& {\displaystyle }{\mathbf{\nabla }}_{\mathrm{n}}\cdot \mathbf{M}=& {\displaystyle }{\mathit{n}}_{\mathrm{b}}\Vert {\mathbf{\nabla }}_{\mathrm{n}}\cdot \mathbf{M}\Vert ={\mathit{n}}_{\mathrm{b}}\left({\mathit{n}}_{\mathrm{b}}\cdot \mathbf{\nabla }\cdot \mathbf{M}\right){\displaystyle }=& {\displaystyle }\left({\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathbf{\nabla }\cdot \mathbf{M},\end{array}$$

and the surface-tangential component

$$\begin{array}{ll}\text{(A7)}& {\displaystyle }{\mathbf{\nabla }}_{\mathrm{s}}\cdot \mathbf{M}=& {\displaystyle }\mathbf{\nabla }\cdot \mathbf{M}-{\mathbf{\nabla }}_{\mathrm{n}}\cdot \mathbf{M}=\mathbf{\nabla }\cdot \mathbf{M}-\left({\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathbf{\nabla }\cdot \mathbf{M}{\displaystyle }=& {\displaystyle }\left(\mathbf{I}-{\mathit{n}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{b}}\right)\cdot \mathbf{\nabla }\cdot \mathbf{M}.\end{array}$$

Surface-normal and tangential components contribute to local acceleration and basal pressure for reasons discussed in terms of gravitational acceleration. Three-dimensional conservation is ensured for fluxes as well if ∇⋅M is calculated conservatively. Finally, we want to note that velocity is a three-dimensional vector field and its direction is not fixed a priori. However, velocity will always be aligned with the surface because only surface-tangential components are present in the respective conservation equation.