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 ﬂow model in ac-cordance 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 prim-itives 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 workﬂow, consisting of meshing, pre-processing, numerical solution and post-processing, is presented. We demonstrate data transfer from and to a geographic information system


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., 2016von Boetticher et al., , 2017Barker and Gray, 2017). Depthintegrated 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 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 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).
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.
Shallow granular flow models have been carefully validated over the last few decades. This includes backcalculations 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(Mergili et al., , 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 userfriendly 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 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 faSav-ageHutterFoam, indicating that the underlying mechanical model is similar to the one of 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. 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.
Also, the export of OpenFOAM results to a GIS for postprocessing 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.

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. 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 threedimensional 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 finitevolume method (Tuković and Jasak, 2012).

Mechanical model
A basic shallow granular flow model can be written in terms of surface partial differential equations as 1 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.
(2) the surface-tangential momentum conservation equation and Eq. (3) its surface-normal counterpart, defined at all points x b on the surface ⊂ R 3 , 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 u(x b ) ∈ R 3 , defined as and the basal pressure p b (x b ). The density ρ is assumed to be constant. Note that the earth pressure theory (e.g. 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 ) ∈ R 3 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 surfacetangential component, and a surface-normal component, 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 surfacetangential component, and a surface-normal component, For simply curved surfaces, the given relation matches the model of Greve et al. (1994), as shown by Rauter and Tuković (2018).

Process models
There are various user-selectable models, describing basal friction τ b (x b ) and entrainment rateq(x b ), to close the system of equations. To reassemble the traditional model (often called Voellmy or Voellmy-Salm model; Christen et al., 2010), as applied, for example, by Fischer et al. (2015), the basal friction is described following Voellmy (1955), 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 0 (10 −7 m s −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. |u| > 100 u 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, aṡ 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, where z is the surface elevation (corresponding to the vertical coordinate in the numerical model) and z 0 is the elevation of a reference station, which has to be provided by the user, alongside with the base value H msc (z 0 ) and the growth rate ∂H msc ∂z . ζ is the angle between the gravitational acceleration and the surface-normal vector. Its further evolution is described by the conservation equation Undershoots, i.e. h msc < 0, are prevented with a regularisation similar to Eq. (9). This can be realised by multiplying the entrainment rateq with h msc h msc +h 0 , where h 0 is a small value, similar to u 0 .

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 (e.g. Jasak, 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.

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 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.

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.
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 threedimensional 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 polygonaldominated 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.
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 (a) (b) (c) (d) 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.
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 postprocessing, 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 modu-larly assembled from various open-source applications. Single modules, such as mesher, solver or friction model, can be replaced easily.

Post-processing and visualisation
Post-processing and visualisation of OpenFOAM simulations is commonly performed using ParaView ® (Ahrens et al., 2005;Ayachit,    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.
To generate regular raster files, the unstructured Open-FOAM 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.

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 0 ) = 1.61 m, z 0 = 1289 m and ∂H msc ∂z = 8 × 10 −4 . Physical parameters to reassemble the runout properly are µ = 0.26, ξ = 8650 m s −2 and e b = 11 500 J kg −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, |u| < 10 −5 m s −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 (SPHmode 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.
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.

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 implementa-  tions (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 polygonaldominated 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).

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 Open-FOAM 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 Open-FOAM 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 threedimensional 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 threedimensional 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.
Code and data availability. The OpenFOAM solver, core utilities and the case study presented are available in the OpenFOAM community repository (https://develop.openfoam.com/Community/ avalanche, last access: 1 March 2018) and integrated as a module within OpenFOAM-v1712. The complete code (based on foamextend-4.0), including python scripts for GIS integration and the simulation set-up including the underlying raw data, is included in the Supplement and available at https://bitbucket.org/matti2/ fasavagehutterfoam (last access: 1 March 2018). g g s g n n b x z Figure A1. Splitting gravitational acceleration into a surfacetangential and surface-normal part with simple projections to the surface-normal vector n b .
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 surfacealigned 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. We begin by splitting a simple vectorial entity, the gravitational acceleration g ∈ R 3 , into a surface-normal component, g n ∈ R 3 , and a surface-tangential component, g s ∈ R 3 , as shown in Fig. A1. The magnitude of the surface-normal component can be calculated using the scalar product and the surface-normal vector: which corresponds to a projection of g on n b . The surfacenormal 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 g n = n b g n = n b (n b · g) = (n b n b ) · g.
The surface-tangential component follows by subtracting the surface-normal component from total gravitational acceleration: g s = g − g n = g − (n b n b ) · g = (I − n b n b ) · g. 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 ∂h u ∂t (Eq. 2). The total gravitational acceleration can be reconstructed by summing up both components: g = g n + g s = (n b n b ) · g + (I − n b n b ) · g = I · g = g, (A4) ensuring perfect conservation of three-dimensional momentum.
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 ∇ · (h u u) or the lateral pressure gradient 1 2 ρ ∇ (p b h) = 1 2 ρ ∇·(I p b h). 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 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 surfacetangential and a surface-normal component which can be treated similarly to gravitational acceleration, yielding the surface-normal component ∇ n · M =n b ∇ n · M = n b (n b · ∇ · M) (A6) = (n b n b ) · ∇ · M, and the surface-tangential component 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.
Author contributions. MR developed the model and the respective code. Simulations have been conducted by MR and AK. AH covered GIS aspects and the production of respective figures. WF conceived the presented investigation, verified the underlying analytical models and supervised the project. All authors discussed the results and contributed to the final paper.