the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
AvaFrame com1DFA (v1.3): a thicknessintegrated computational avalanche module – theory, numerics, and testing
Matthias Tonnel
Anna Wirbel
JanThomas Fischer
Simulation tools are important to investigate and predict mobility and the destructive potential of gravitational mass flows (e.g., snow avalanches). AvaFrame – the open avalanche framework – offers wellestablished computational modeling approaches, tools for data handling and analysis, and readytouse modules for evaluation and testing. This paper presents the theoretical background, derivation, and model verification for one of AvaFrame's core modules, the thicknessintegrated computational model for avalanches with flow or mixed form of movement, named com1DFA. Particular emphasis within the description of the utilized numerical particle–grid method is given to the computation of spatial gradients and the accurate implementation of driving and resisting forces. The implemented method allows us to provide a time–space criterion connecting the numerical particles, grid, and time discretization. The convergence and robustness of the numerical implementation is checked with respect to the spatiotemporal evolution of the flow variables using tests with a known analytical solution. In addition, we present a new test for verifying the accuracy of the numerical simulation in terms of runout (angle and distance). This test is derived from the total energy balance along the centerofmass path of the avalanche. This article, particularly in combination with the code availability (opensource code repository) and detailed online documentation provides a description of an extendable framework for modeling and verification of avalanche simulation tools.
 Article
(5260 KB)  Fulltext XML

Supplement
(17270 KB)  BibTeX
 EndNote
Simulation tools for gravitational mass flows – with a focus on snow avalanches in this article – are in great demand for operational engineering applications and scientific model development and are gaining increasing attention in academic education. Each of these application requires different outputs. For operational engineering applications, the runout outcome for different scenarios is usually of highest interest. Scientific applications aim at better understanding the processes and will require outputs such as flow variable evolution. Existing tools for simulating snow avalanches cover a wide range of numerical implementations and vary from proprietary (e.g., Christen et al., 2010; Sampl and Zwinger, 2004; Zugliani and Rosatti, 2021; Li et al., 2021) to opensource software (e.g., Hergarten and Robl, 2015; Mergili et al., 2017; Rauter et al., 2018). The latter ones are generally more focused towards scientific and academic issues, whereas the first are more geared towards operational applications. AvaFrame – the open avalanche framework – strives to fill the gap between operations and scientific development combining over a decade of operational application (using, e.g., SamosAT; see Sampl and Zwinger, 2004; Fischer et al., 2014) with an opensource scientific development environment. Using a modular structure AvaFrame adds indepth testing and analysis modules to the core flow modules. Further modules provide interfaces for visualization and geodata handling for all kinds of existing and emerging simulation tools. It enables us to combine, further develop, and extend the different tools to best suit the users needs.
At their core, avalanche simulation tools are based on a large variety of flow models, differing in their basic assumptions (what physical processes they include, degree of complexity), mathematical derivation, and/or their numerical implementation. These range from Eulerian methods (Christen et al., 2010; Mergili et al., 2017; Rauter and Tuković, 2018; Zugliani and Rosatti, 2021) using spatially fixed meshes to Lagrangian methods (Sampl and Zwinger, 2004) where the mass is discretized. Some methods are combinations of Eulerian and Lagrangian approaches, such as the Material Point method (Stomakhin et al., 2013). AvaFrame's com1DFA denseflow avalanche (DFA) module is based on a flow model that is derived from the thickness integration of the conservation equations of mass and momentum. Classic shallowwater models, e.g., SaintVenant, integrate in the direction of gravity, often called “depth” integration. Other approaches, such as Savage–Hutter models (Savage and Hutter, 1989) integrate in the slopenormal direction but also call it “depth” integrated. To be more consistent with operational terminology, we propose calling the integration in the slopenormal direction “thickness” integration. In this way, we can highlight the special case of gravitational mass flows in steep terrain.
The resulting equations are solved using a mixed particle–grid method, in which mass is tracked using particles. Pressure gradients are computed using a smooth particle hydrodynamics (SPH) method adapted to steep terrain and flow thickness is computed on the grid (Sampl and Zwinger, 2004; Monaghan, 2005; Liu and Liu, 2010; Granig et al., 2016). To avoid nonphysical behavior at starting and stopping, com1DFA applies the method proposed in MangeneyCastelnau et al. (2003), allowing for a friction balanced starting and stopping of the flow.
Verifying and validating the methods applied in our implementation is a challenging but crucial step, as it is for all simulation tools. Verification is done by comparing the numerical model results to an analytical solution (e.g., Zugliani and Rosatti, 2021; Rauter and Tuković, 2018). Validation can be tackled in different ways, either by comparing the model results to observations (e.g., Christen et al., 2010) or by comparing them to the results of already trusted numerical models. In this article, the focus is on model verification, and the numerical model results are compared to tests with known (semi)analytical solutions using two different approaches.
The flow variable tests (Hutter et al., 1993; Faccanoni and Mangeney, 2013) allow us to investigate the local spatiotemporal evolution of flow thickness and velocity. This enables us to verify the proper implementation of the pressure gradient and friction force computation among others. In contrast, the energy line test is based on the investigation of the total energy balance (Körner, 1980). It focuses on the accuracy of the global kinetic energy (velocity) along the path and the corresponding centerofmass runout. Thereby, the proper implementation of stopping behavior, that is to say the proper balancing of driving and friction forces, can be assessed. It also provides a test where the runout, a quantity which is important for operational applications, is tested. We explore and explain the limitations of these two approaches.
The article is structured as follows. Section 2 summarizes the underlying flow model, including fundamental assumptions and derivations of the thicknessintegrated equations, building the foundation for the gravitational mass flow simulations. In Sect. 3 the temporal and spatial discretizations of the model equations, employing a particle–grid approach, are described. The implementation in the AvaFrame computational module com1DFA is outlined in Sect. 4. Model verification tests are presented in Sect. 5, targeting the correct implementation of the mathematical model and the convergence and robustness of the numerical model code. In addition to employing test cases with a known (semi)analytical solution, a new energy line test is introduced in Sect. 5.2.
Besides the indepth description within this article, Oesterle et al. (2022) provide a combination of code and corresponding documentation. Users find more information according to their individual scientific, operational, or educational focus, and the reader is invited to contribute to the future development. It is important to note that this article presents the latest development state of the com1DFA module. It differs slightly from the implementation of the com1DFA module used for operational purposes, which is described in the online AvaFrame documentation. For example, differences include improvements of the SPH gradient computation method and how friction forces are taken into account.
In this section, the mathematical model and associated equations used to simulate DFA processes are presented. The derivation is based on the thickness integration of the threedimensional Navier–Stokes equations, using a Lagrangian approach with a terrainfollowing coordinate system. The equations are simplified using the assumption of shallow flow on a mildly curved topography, meaning flow thickness is considerably smaller than the length and width of the avalanche, and it is considerably smaller than the topography curvature radius.
We consider snow as the material; however, the choice of material does not influence the derivation in the first place. We assume constant density ρ_{0} and a flow on a surface 𝒮, subjected to the gravity force and friction on the surface 𝒮. If needed, additional processes such as entrainment or other external effects can be taken into account. These processes are included in the following derivations but will not be considered for model verification (Sect. 5), as test cases with an analytical solution are only available for simplified conditions where entrainment or any additional forces are discarded. The mass conservation equation applied to a Lagrangian volume of material V(t) reads:
where the source term q^{ent} represents the snow entrainment rate (incoming mass flux). The momentum conservation equation applied to the same volume of material V(t) reads as follows:
where u is the fluid velocity and g the gravity acceleration. $\mathit{\sigma}=p\text{\U0001d7d9}+\mathbf{T}$ represents the stress tensor, where 𝟙 is the identity tensor, p the pressure, and T the deviatoric part of the stress tensor. n is the normal vector on ∂V(t). F^{ext} represents additional forces due to snow entrainment (force needed to break and compact the entrained snow) or due to added resistance (trees, obstacles, etc.). F^{ext} is assumed to be surface parallel. Problemspecific assumptions are needed to solve the mass and momentum conservation equations (closure equation on the stress tensor in Sect. 2.2). These constitutive equations are introduced in the following sections alongside a local coordinate system and boundary conditions.
2.1 Natural coordinate system (NCS) and thicknessintegrated quantities
In order to solve the previously described equations, a local coordinate system is introduced. The avalanche flows on a surface 𝒮, a 2D manifold embedded in the 3D Euclidian space. Different approaches exist to define a coordinate system on this curved surface. Some define a tangent space in every point on 𝒮 based on the coordinate lines, which leads to an orthogonal but not orthonormal coordinate system for a curved surface 𝒮 (e.g., Luca et al., 2016). Instead of this and because of the Lagrangian approach used here, we define a local coordinate system in the tangent plane to 𝒮 at any point using the velocity direction and the normal to the surface 𝒮 at this position. This results in a timedependent orthonormal coordinate system that is advected along with the flow, referred to as the natural coordinate system (NCS).
A control volume V(t) is assumed to be a small truncated pyramid shape extending from the bottom surface 𝒮_{b} (lying on the topography 𝒮) up to the free surface in the surfacenormal direction N^{b} as illustrated in Fig. 1. With the assumption of moderately curved surfaces, this is close to being a prism shape since the normals of the lateral surfaces are almost parallel to the bottom. Note that the bottom surface 𝒮_{b} of area A^{b} has no predefined shape. The octagonal shape used in Fig. 1 is just one possible option.
The normal vector N^{b} to the bottom surface is pointing upwards, whereas ${\mathit{n}}^{\mathrm{b}}={\mathit{N}}^{\mathrm{b}}$ is the bottom normal vector to the Lagrangian control volume (pointing out of the volume).
where J is the transformation matrix from the Cartesian coordinate system to the NCS. The NCS is an orthonormal coordinate system $({\mathit{v}}_{\mathrm{1}},{\mathit{v}}_{\mathrm{2}},{\mathit{v}}_{\mathrm{3}})$ aligned with the bottom surface. ${\mathit{v}}_{\mathrm{3}}={\mathit{N}}^{\mathrm{b}}={\mathit{n}}^{\mathrm{b}}$ is the normal vector to the bottom surface pointing upwards. v_{1} is pointing in the direction of the thicknessintegrated fluid velocity $\stackrel{\mathrm{\u203e}}{\mathit{u}}$ (introduced below).
In the case of shallow flow on weakly curved surfaces, $\mathrm{det}\left(\mathbf{J}\right)=(\mathrm{1}{\mathit{\kappa}}_{\mathrm{1}}{x}_{\mathrm{3}})(\mathrm{1}{\mathit{\kappa}}_{\mathrm{2}}{x}_{\mathrm{3}})\approx \mathrm{1}$. κ_{{1,2}} represents the surface curvature in v_{{1,2}} directions, while x_{3} is the elevation from the bottom surface in the direction N^{b}. This approximation is valid if the curvature radius is much larger than the flow thickness h. In this case, the control volume reads as follows:
The following volume (indicated by the superscript $\stackrel{\mathrm{\u0303}}{\square}$), area (indicated by the superscript $\widehat{\square}$), and thickness (indicated by the superscript $\stackrel{\mathrm{\u203e}}{\square}$) averaged quantities are introduced (where f is a scalar or vector function on Ω⊂ℛ^{3}):
Note that $\stackrel{\mathrm{\u0303}}{f}\left(\mathrm{0}\right)=\stackrel{\mathrm{\u0303}}{f}$ and $\stackrel{\mathrm{\u203e}}{f}\left(\mathrm{0}\right)=\stackrel{\mathrm{\u203e}}{f}$. When the control volume goes to 0, i.e., basal area goes to a point, $\stackrel{\mathrm{\u0303}}{f}\stackrel{{A}^{\mathrm{b}}\to \mathrm{0}}{\u27f6}\stackrel{\mathrm{\u203e}}{f}$ and $\widehat{f}\stackrel{{A}^{\mathrm{b}}\to \mathrm{0}}{\u27f6}f$. Also, note that we assume integration on a tangent plane being equivalent to integration on a small surface in the manifold defined by the terrain (or 3D space in which the terrain is embedded). This is justified by the relative smallness of the basal particle surface compared to the curvature.
The NCS has some interesting properties that will be useful for projecting and solving the equations. Because of the orthogonality of this NCS, we have ${\mathit{v}}_{i}\cdot {\mathit{v}}_{j}={\mathit{\delta}}_{ij},\phantom{\rule{0.125em}{0ex}}\mathit{\{}i,j\mathit{\}}\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3}{\mathit{\}}}^{\mathrm{2}}$, which gives after time derivation:
meaning that
It is possible to express $\frac{\mathrm{d}{\mathit{v}}_{\mathrm{1}}}{\mathrm{d}t}$ in terms of $({\mathit{v}}_{\mathrm{1}},{\mathit{v}}_{\mathrm{2}},{\mathit{v}}_{\mathrm{3}})$ and using orthogonality of $\frac{\mathrm{d}{\mathit{v}}_{i}}{\mathrm{d}t}$ and v_{i}:
The derivative of the thicknessintegrated velocity decomposes to
2.2 Boundary conditions
To complete the conservation Eqs. (1) and (2) the following boundary conditions at the bottom (𝒮_{b}) and free (𝒮_{fs}) surfaces are introduced. σ_{s} and σ_{b} represent the restriction of σ to the free surface 𝒮_{fs} and bottom surface 𝒮_{b}, respectively:

tractionfree top surface –
$$\begin{array}{}\text{(11)}& {\mathit{\sigma}}_{\mathbf{s}}\cdot {\mathit{n}}_{\mathrm{s}}=\mathbf{0}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathcal{S}}_{\mathrm{fs}};\end{array}$$ 
impenetrable bottom surface without detachment –
$$\begin{array}{}\text{(12)}& {\mathit{u}}^{\mathrm{b}}\cdot {\mathit{n}}^{\mathrm{b}}=\mathbf{0}\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathcal{S}}_{b};\end{array}$$ 
bottom friction law –
$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathit{\tau}}^{\mathrm{b}}& ={\mathit{\sigma}}^{\mathrm{b}}\cdot {\mathit{n}}^{\mathrm{b}}\left(\right({\mathit{\sigma}}^{\mathrm{b}}\cdot {\mathit{n}}^{\mathrm{b}})\cdot {\mathit{n}}^{\mathrm{b}}){\mathit{n}}^{\mathrm{b}}\\ & =f({\mathit{\sigma}}^{\mathrm{b}},\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{\mathit{u}},\phantom{\rule{0.125em}{0ex}}h,\phantom{\rule{0.125em}{0ex}}{\mathit{\rho}}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}t,\phantom{\rule{0.125em}{0ex}}\mathit{x})\\ & =f({\mathit{\sigma}}^{\mathrm{b}},\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{\mathit{u}},\phantom{\rule{0.125em}{0ex}}h,\phantom{\rule{0.125em}{0ex}}{\mathit{\rho}}_{\mathrm{0}},t,\phantom{\rule{0.125em}{0ex}}\mathit{x}){\mathit{v}}_{\mathrm{1}}\phantom{\rule{1em}{0ex}}\text{on}{\mathcal{S}}_{b}.\end{array}\end{array}$$
2.3 Constitutive relation: friction force
To close the momentum equation, a constitutive equation describing the basal shear stress tensor τ^{b} as a function of the avalanche flow state is required:
The model verification tests (Sect. 5) are carried out with a Mohr–Coulomb friction model, which describes the friction interaction between two solids. The bottom shear stress reads as follows:
where δ is the friction angle and μ=tan δ is the friction coefficient. The bottom shear stress linearly increases with the normal stress component p^{b}.
With Mohr–Coulomb friction an avalanche starts to flow if the slope inclination exceeds the friction angle δ. In the case of an infinite slope of constant inclination, the avalanche velocity would increase indefinitely. However, because of its relative simplicity, the Mohr–Coulomb friction model is convenient for deriving analytical solutions and testing numerical implementations. For a more detailed discussion of friction laws and their applicability, we refer to Salm and Gubler (1985); Gubler (1987); Gauer (2014).
Different friction models accounting for the influence of flow velocity, flow thickness, etc., have been proposed. Three friction models are available in the com1DFA module. First, a Coulomb one, which is used in this paper, second a Voellmy friction model (Voellmy, 1955), and third the samosAT friction model, which is the one used for hazard mapping by Austrian federal agencies (Sampl, 2007).
2.4 Expression of surface forces in the NCS
Taking advantage of the NCS and using the boundary conditions, it is possible to split the surface forces into bottom, lateral, and free surface forces and perform further simplifications:
Using the notations introduced in Sect. 2.1 and the decomposition of the stress tensor, the bottom force can be expressed as a surfacenormal component and a surfacetangential one:
where τ^{b} is the basal friction term (introduced in Sect. 2.2). Applying Green's theorem, the lateral force reads as follows:
Equations (17) and (18) represent the thicknessintegrated form of the surface forces and can now be used to write the thickness integrated momentum equation.
2.5 Thicknessintegrated momentum equation
Using the definitions of average values given in Sect. 2.1 and the decomposition of the surface forces given by Eqs. (17) and (18), the momentum equation reads as follows:
which leads to
The lateral shear stress term is neglected because of its relative smallness in comparison to the other terms as shown by the dimensional analysis carried out in Gray and Edwards (2014). The mass conservation reads as follows:
Using the approximations from Sect. 2.1, the equation of motion becomes
where all quantities are evaluated at the center of the basal area (point O in Fig. 1). This equation is projected in the normal direction v_{3}=N^{b} to get the expression of the basal pressure p^{b}. The projection of this same equation on the tangential plane leads to the differential equations satisfied by $\stackrel{\mathrm{\u203e}}{\mathit{u}}$.
2.5.1 Pressure distribution, thicknessintegrated pressure, and pressure gradient
We can project the momentum equation (Eq. 22), using the volume between x_{3} and the surface h, in the normal direction (${\mathit{v}}_{\mathrm{3}}={\mathit{N}}^{\mathrm{b}}={\mathit{n}}^{\mathrm{b}}$). Applying the properties of the NCS (Eq. 10) the surfacenormal component of Eq. (22) reads as follows:
Neglecting the normal component of the pressure gradient gives the expression for pressure. Under the condition that $\stackrel{\mathrm{\u203e}}{\mathit{u}}$ is independent of x_{3}, pressure follows a linear profile from the bottom surface to the free surface:
Note that the bottom pressure should always be positive. A negative pressure is nonphysical and means that the material is not in contact with the bottom surface anymore. This can happen in the case of large velocities on convex topography. If so, the material should be in a freefall state until it gets back in contact with the topography. A description on how this is handled within the numerical implementation can be found in Sect. 4.3.
Using Eq. (24), it is possible to express the thicknessintegrated pressure $\stackrel{\mathrm{\u203e}}{p}$:
where g^{eff} is the effective normal acceleration acting on the volume, including the normal component of gravity and a curvature component. Because of the utilized Lagrangian approach, the curvature terms are expressed as the temporal derivative of the normal vector $\frac{\mathrm{d}{\mathit{N}}^{\mathrm{b}}}{\mathrm{d}t}$, effectively computing the curvature along the particle trajectories. The resulting curvatures are equivalent to the ones obtained through Eulerian approaches (Fischer et al., 2012) but does not require the computation of the related κ_{{1,2}}.
The expression of the thicknessintegrated pressure is used to derive the pressure gradient $\mathbf{\nabla}\left(h\stackrel{\mathrm{\u203e}}{p}\right)$. Assuming g^{eff} to be locally constant (i.e., effectively neglecting curvature; otherwise g^{eff} would remain inside the gradient operator) leads to
2.5.2 Tangential momentum equation
Using the derived expression of the thicknessintegrated pressure (Eq. 26), we project the momentum balance (Eq. 22) in the tangent plane, which leads to the following equation:
where ${\mathbf{\nabla}}_{s}=\mathbf{\nabla}(\mathbf{\nabla}\cdot {\mathit{N}}^{\mathrm{b}}){\mathit{N}}^{\mathrm{b}}$ and ${\mathit{g}}_{\mathrm{s}}=\mathit{g}(\mathit{g}\cdot {\mathit{N}}^{\mathrm{b}}){\mathit{N}}^{\mathrm{b}}$ are the tangential component of the gradient operator and of the gravity acceleration, respectively.
After replacing the velocity derivative component in the normal direction with the expression developed in Eq. (23), Eq. (27) reads as follows:
The curvature acceleration is in the normal direction to the tangent plane in order to keep the flow on the surface.
In the previous section, the equation of motion was derived using a Lagrangian approach. In order to solve this set of equations numerically, we employ a mix of particle and grid approaches. We discretize the material into numerical particles and solve the equation of motion, with the total avalanche mass being the sum of the mass associated with each particle. The grid is used to compute several parameters that are required for the computations, e.g., surfacenormal vectors and flow thickness. Combining both approaches allows us to best exploit the advantages of each. The particle approach is used to track the mass, compute the curvature terms and the gradient of the flow thickness, and update the particle positions. The grid is used to handle the topography information and compute the flow thickness and artificial viscosity. We found this to help with numerical stability, and it is more efficient as it decreases the required number of particles. A theoretical convergence criterion is described in the last section.
3.1 Interpolation between particle and grid values
Topography information is usually provided in a raster format which corresponds to a regular rectilinear mesh on the global horizontal x–y plane, henceforth referred to as grid. In order to get information on the surface elevation and normal vectors, the topography information needs to be interpolated at the particle locations, and this needs to be repeated at each time step since the particles are moving. Similarly, the particle properties such as mass or momentum, which translate into flow thickness and velocity, also need to be interpolated onto the grid. Grid velocity is required to compute the artificial viscosity term, ensuring numerical stability; see Sect. 4.2. Grid flow thickness is used to compute the particle flow thickness, which is required for computing the frictional forces. The mesh being regular and rectilinear, we use bilinear interpolation for simplicity and efficiency. This also ensures the conservation of mass or momentum when interpolating from the particles to the grid and back. The properties of the grid and the interpolation method are detailed in the AvaFrame documentation (https://docs.avaframe.org/en/latest/DFAnumerics.html#interpolation, last access: 5 October 2023).
3.2 The particle momentum equation
Discretizing the material into particles (particle quantities are denoted by the subscript k; e.g., m_{k}=ρ_{0}V_{k} is the mass of particle k) leads to the following continuity equation:
where ${A}_{k}^{\text{ent}}$ is the interface area between the particle and the entrainable material while ${q}_{k}^{\text{ent}}$ represents the flux of entrained mass.
By assuming that the Lagrangian control volume V can be represented by a particle, we can derive the particle momentum equation in the normal direction and in the tangent plane:
In this equation, flow thickness gradient, basal friction, and curvature acceleration terms need to be further developed and discretized.
3.3 Flow thickness and its gradient
3.3.1 Flow thickness gradient computation using SPH
To assess the flow thickness gradient, we employ an SPH method (smoothed particles hydrodynamics method (Liu and Liu, 2010)), where the gradient is directly derived from the particles and does not require any mesh. In contrast, a mesh method or an MPM (material point method) would directly use a mesh formulation to approximate the gradient or interpolate the particle properties on an underlying mesh and then approximate the gradient of the flow thickness using a mesh formulation. There are two main advantages of using SPH. First is the computational efficiency, since only the particles of interest are computed without needing to compute the area of the whole terrain (as with Eulerian methods). Second, mass transfer is not required because it is handled by the particles directly, making the implementation easier (Granig et al., 2016).
In theory, an SPH method does not require any mesh to compute the gradient. However, applying this method requires finding neighbor particles. This process can be sped up with the help of an underlying grid; different neighbor search methods are presented in Ihmsen et al. (2014), and a “uniform grid method” is used in this paper.
The SPH method is used to express a quantity (the flow thickness in our case) and its gradient at any particle location as a weighted sum of its neighbors' properties. The principle of the method is described well in Liu and Liu (2010), and the basic formula reads as follows:
where W represents the SPHkernel function (we employ the spiky kernel; see Eq. B2) and the subscript l enumerates the neighboring particles. This kernel function is designed to satisfy the unity condition, be an approximation of the Dirac function and have a compact support domain (Liu and Liu, 2010).
This method is usually either used in a 3D space, in which particles move freely and where the weighting factor for the summation is the volume of the particle, or on a 2D horizontal plane, where the weighting factor for the summation is the area of the particle and the gradient is 2D. Here we want to compute the gradient of the flow thickness on a 2D surface (the topography) embedded in 3D space. The method used is analogous to the SPH gradient computation on the 2D horizontal plane but the gradient is 3D and tangent to the surface (colinear to the local tangent plane). The theoretical derivation in Appendix B2 shows that the SPH computation is equivalent to applying the 2D SPH method in the local tangent plane instead of in the horizontal plane. This leads to the following SPH expression of the flow thickness gradient:
3.3.2 Flow thickness computation
The particle flow thickness is computed with the help of the grid. The mass of the particles is interpolated onto the grid using a bilinear interpolation method (described in Sect. 3.1). Then, dividing the mass at the grid cells by the area of the grid cells, while taking the slope of the cell into account, returns the flow thickness field on the grid. This is interpolated back to the particles, which leads to the particle flow thickness property.
We do not compute the flow thickness directly from the particle properties (mass and position) using an SPH method because it induced instabilities. Indeed, the cases where too few neighbors are found lead to very small flow thickness, which becomes an issue for flowthicknessdependent friction laws. Note that using such an SPH method would lead to a full particle method. But since the flow thickness is only used in some cases for the friction force computation, using the previously described grid method should not affect the computation significantly.
3.4 Friction force discretization
The Coulomb friction force term in Eq. (30) for a particle reads:
This relation stands if the particle is moving. The starting and stopping processes satisfy a different equation and are handled differently in the numerical implementation (using the same equation would lead to nonphysical behavior). This is described in more detail in Sect. 4.5.
3.5 Time discretization
The momentum equation is solved numerically in time using an Euler time scheme. The time derivative of any quantity f is approximated by
where Δt represents the time step and f^{n}=f(t^{n}), t^{n}=nΔt. For the velocity this reads as follows:
The positions of the particles are updated using a centered Euler scheme:
The forces are taken into account in two subsequent steps as forces acting on the particles can be sorted into driving forces and friction forces. Friction forces act against the particle motion, only affecting the magnitude of the velocity. They can in no case become driving forces. This is why in a first step the velocity is updated with the driving forces before updating in a second step the velocity magnitude applying the friction force.
3.6 Convergence
We are looking for a criterion that relates the properties of the spatial and temporal discretization to ensure convergence of the numerical solution. Simply decreasing the time step and increasing the spatial resolution, by decreasing the grid cell size and kernel radius and increasing the number of particles, does not ensure convergence. Ben Moussa and Vila (2000) analyzed hyperbolic one and twodimensional, nonlinear transport equations with a particle and SPH method and showed that the kernel radius size cannot be varied independently from the time step and number of particles. Indeed, they showed that the numerical solution converges towards the solution of the equation under the following condition:
where r_{part} represents the “size” of a particle, r_{kernel} represents the SPH kernel radius, dt is the time step and C a constant. The conditions in Eq. (37) mean that both r_{part} (particle size) and r_{kernel} (SPH kernel radius) need to go to zero but also that the particle size needs to go faster to zero than the SPH kernel radius. Finally, the time step needs to go to zero at the same rate as r_{kernel}. The particle size can be expressed as a function of the SPH kernel radius:
where the particle basal area was assumed to be a circle. Note that this does not affect the results except adding a different shape factor in front of this expression. n_{ppk} is the number of particles per kernel support area and defines the density of the particles when initializing a simulation. Let n_{ppk} be defined by a reference number of particles per kernel support area ${n}_{\text{ppk}}^{\mathrm{0}}>\mathrm{0}$, a reference kernel radius ${r}_{\text{kernel}}^{\mathrm{0}}>\mathrm{0}$, and a real exponent α:
This leads to
Replacing r_{part} with the previous equation in Eq. (37) leads to the following condition:
This brings us to the following choice:
which satisfies the convergence criterion if
Note that this criterion (which is similar to a Courant–Friedrichs–Lewy (CFL) condition) leaves some freedom of choice for the exponent α and that there are no constraints on the reference kernel radius ${r}_{\text{kernel}}^{\mathrm{0}}$ and reference number of particles per kernel radius ${n}_{\text{ppk}}^{\mathrm{0}}$. Nevertheless, it seems appropriate to require a minimum number of particles per kernel radius so that enough neighbors are available to get a reasonable estimate of the gradient. These parameters should be adjusted according to the expected accuracy of the results and/or available compute power. Determining the optimal parameter values for α, ${r}_{\text{kernel}}^{\mathrm{0}}$, and ${n}_{\text{ppk}}^{\mathrm{0}}$, for example, according to a user's needs in terms of accuracy and computational efficiency, requires a specific and detailed investigation of the considered case. In Sect. 5, we will explore model convergence using the condition Eq. (43) with different values of α.
In this section, the numerical implementation and algorithm of the com1DFA module are described. The following sections are organized following the workflow used in the com1DFA code, which is also illustrated in Fig. 2. First, the release mass is discretized into particles and the grid is initialized. As a result of the partial differential equations considered and the time discretization used, stability issues might arise. Hence, artificial viscosity is added in order to ensure the stability of the solution. As a next step, driving forces (including curvature effects) are accounted for. Friction forces are taken into account subsequently in order to ensure proper starting and stopping behavior. Finally, a reprojection step is needed to ensure that particles lie on the topography and that particle velocities are tangential to the topography. For simplicity and because they are not considered in the verification tests in Sect. 5, entrainment and added resistance effects are not included in what follows. Additional information about entrainment or resistance forces is available in the theory section (https://docs.avaframe.org/en/theorypapercode/theoryCom1DFA.html#entrainment, last access: 5 October 2023) of the AvaFrame documentation.
4.1 Initialization
To start a simulation with com1DFA, input information about topography, material properties, and initial conditions is required. Topography is described by a DEM (digital elevation model) using the ESRI ASCII format. It supplies a grid of altitude values, preceded by a header with information about the number of data rows and columns, coordinates of the center or corner of the lowerleft cell, the edge length of the quadratic cells, and the code used for missing values. The material is characterized by its density and some friction properties. The initial condition is given by release areas, polygons describing the initial material location, and the release thickness, in our case measured in the surfacenormal direction. It is possible to provide several polygons with different initial thickness values.
Then the material is discretized into particles. The field of normal vectors to the surface is computed from the input DEM and the different grid fields are initialized. The details of the initialization process are given in the initialization section (https://docs.avaframe.org/en/theorypapercode/com1DFAAlgorithm.html#initialization, last access: 5 October 2023) of the AvaFrame documentation.
4.2 Numerical stability
Because the lateral shear force term was removed when deriving the model equations (because of its relative smallness, Gray and Edwards, 2014), Eq. (22) is hyperbolic. Hyperbolic systems have the characteristic of carrying discontinuities or shocks which will cause numerical instabilities. They would fail to converge if, for example, an Euler forwardintime scheme is used (LeVeque, 1990). Several methods exist to stabilize the numerical integration of a hyperbolic system of differential equations. All add some upwinding in the discretization scheme. Some methods tackle this problem by introducing some upwinding in the discretization of the derivatives (Harten et al., 1983; Harten and Hyman, 1983). Others introduce some artificial viscosity (as in Monaghan, 1992), which is also implemented in com1DFA. The following artificial viscosity force acting on particle k is added to stabilize the momentum equation:
where the velocity difference reads $\mathit{d}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}={\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}{\widehat{\stackrel{\mathrm{\u203e}}{\mathit{u}}}}_{k}$ (${\widehat{\stackrel{\mathrm{\u203e}}{\mathit{u}}}}_{k}$ represents the averaged velocity of the neighbor particles and is practically the grid velocity interpolated at the particle position). C_{Lat} is a coefficient that controls the viscous force. It would be the equivalent of C_{Drag} in the case of the drag force. C_{Lat} is a numerical parameter that depends on the grid size.
In this expression, let ${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}^{n}$ be the velocity at the beginning of the time step and ${{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}^{n+\mathrm{1}}}^{\u25b4}$ be the velocity after adding the numerical viscosity (Fig. 2). In the norm term $\left\left\mathit{d}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}\right\right$ the particle and grid velocity at the beginning of the time step are used. This ensures explicit time discretization with respect to the norm $\left\left\mathit{d}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}\right\right$ and $\widehat{\stackrel{\mathrm{\u203e}}{\mathit{u}}}$ in $\mathit{d}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}$. In contrast, an implicit formulation is used in $\mathit{d}{\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}$ because the new value of the velocity is used there. The artificial viscosity force now reads:
Updating the velocity then gives
with
This approach to stabilize the momentum equation (Eq. 30) is not optimal for different reasons. Firstly, it introduces a new coefficient C_{vis}, which is not a physical quantity and will require calibration. Secondly, it artificially adds a force that should be described physically. So it would be more interesting to take the physical force into account in the first place.
Potential solutions could be taking the physical shear force into account, using for example the μ(I)rheology (Gray and Edwards, 2014; Baker et al., 2016). Another option would be to replace the artificial viscosity with a purely numerical artifact aiming to stabilize the equations, such as an SPH version of the Lax–Friedrichs scheme as presented in Ata and Soulaïmani (2005).
4.3 Curvature acceleration term
The last term of the particle momentum equation (Eq. 30), as well as the effective gravity g^{eff}, are the final terms to be discretized before the time integration. In both of these terms, the remaining unknown is the curvature acceleration term ${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}\cdot \frac{\mathrm{d}{\mathit{v}}_{\mathrm{3},k}}{\mathrm{d}t}$. Using the forward Euler time discretization for the temporal derivative of the normal vector v_{3,k} gives
${\mathit{v}}_{\mathrm{3},k}^{n}$ is a known quantity, the normal vector of the bottom surface at ${\mathit{x}}_{k}^{n}$ which is interpolated from the grid normal vector values at the position of the particle k at time t^{n}. ${\mathit{v}}_{\mathrm{3},k}^{n+\mathrm{1}}$ is unknown since ${\mathit{x}}_{k}^{n+\mathrm{1}}$ is not known yet, hence we estimate ${\mathit{x}}_{k}^{n+\mathrm{1}}$ based on the position ${\mathit{x}}_{k}^{n}$ and the velocity at t^{n}:
This position at t^{n+1} is projected onto the topography, and ${\mathit{v}}_{\mathrm{3},k}^{n+\mathrm{1}}$ can be interpolated from the grid normal vector values.
Note that the curvature acceleration term is needed to compute the bottom pressure (Eq. 24), which is used for the bottom friction computation and for the pressure gradient computation. The curvature acceleration term can lead to a negative value, which means detachment of the particles from the bottom surface. In com1DFA, surface detachment is not allowed, and if pressure becomes negative, it is set back to zero forcing the material to remain in contact with the topography.
4.4 Driving forces
Adding the driving forces is done after adding the artificial viscosity as described in Fig. 2. The velocity is updated as follows ( is the velocity after taking the driving force into account):
where the flow thickness gradient term is computed using the SPH formulation in Eq. (32).
4.5 Friction forces
The friction force related to the bottom shear force needs to be taken into account in the momentum equation and the velocity needs to be updated accordingly. Friction force acts against motion, hence it only affects the magnitude of the velocity and cannot be a driving force (MangeneyCastelnau et al., 2003). Moreover, the friction force magnitude depends on the particle state, i.e., whether it is flowing or at rest. If the velocity of the particle k is after adding the driving forces, adding the friction force leads, depending on the sign of $\frac{{m}_{k}\left\left\textcolor{red}{}\right\right}{\mathrm{\Delta}t}{\left\left{\mathbf{F}}_{k}^{\text{fric}}\right\right}_{\text{max}}$, to

$\left\left{\mathbf{F}}^{\text{fric}}\right\right={\left\left{\mathbf{F}}_{k}^{\text{fric}}\right\right}_{\text{max}}$ and ${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}^{n+\mathrm{1}}=\textcolor{red}{}\left(\mathrm{1}\frac{\mathrm{\Delta}t}{{m}_{k}}\frac{{\left\left{\mathbf{F}}_{k}^{\text{fric}}\right\right}_{\text{max}}}{\left\left\textcolor{red}{}\right\right}\right)$, if $\frac{{m}_{k}\left\left\textcolor{red}{}\right\right}{\mathrm{\Delta}t}>{\left\left{\mathbf{F}}_{k}^{\text{fric}}\right\right}_{\text{max}}$

$\left\left{\mathbf{F}}_{k}^{\text{fric}}\right\right\le {\left\left{\mathbf{F}}_{k}^{\text{fric}}\right\right}_{\text{max}}$, and the particle stops moving (${\stackrel{\mathrm{\u203e}}{\mathit{u}}}_{k}^{n+\mathrm{1}}=\mathrm{0}$) before the end of the time step if $\frac{{m}_{k}\left\left\textcolor{red}{}\right\right}{\mathrm{\Delta}t}\le {\left\left{\mathbf{F}}_{k}^{\text{fric}}\right\right}_{\text{max}}$.
This method prevents the friction force from becoming a driving force and nonphysically changing the direction of the velocity. This would lead to oscillations of the particles instead of stopping. Adding the friction force following this approach (MangeneyCastelnau et al., 2003) allows the particles to start and stop flowing properly.
4.6 Reprojection
The last term in Eq. (30) (accounting for the curvature effects) adds a nontangential component allowing the new velocity to lie in a different plane than the one from the previous time step. This enables the particles to follow the topography. Because the curvature term was only based on an estimation (see Sect. 4.3), the new particle position is not necessarily on the topography, and the new velocity does not necessarily lie in the tangent plane at this new position. Furthermore, in cases with a strong convex curvature and high velocities, the particles can theoretically be in a freefall state (detachment), as mentioned in Sect. 2.5.1. Thus, com1DFA does not allow detachment of the particles, and the particles are forced to stay on the topography. This is a limitation of the model or method, which will lead to nonphysical behavior in special cases (material flowing over a cliff). In both of the previously mentioned cases, the particle positions are projected back onto the topography and the velocity direction is corrected to be tangential to the topography. The position reprojection is done using an iterative method that attempts to conserve the distance traveled by each particle between t^{n} and t^{n+1}. The velocity reprojection changes the direction of the velocity, but its magnitude is conserved.
In this section, the numerical implementation of the mathematical model is tested. We present different tests where, for specific conditions, a (semi)analytical solution exists. The tests described here are implemented in the ana1Tests module from AvaFrame. In the first set of tests, the flow variable tests, we compare the temporal and spatial evolution of the flow thickness (h) and the thicknessintegrated momentum flux ($h\stackrel{\mathrm{\u203e}}{\mathit{u}}$) of the com1DFA simulation results to a (semi)analytical solution. With these tests, we aim to verify the numerical discretization and implementation of the solver and check the validity of the convergence criterion described in Sect. 3.6.
In the second test, the energy line test, we investigate global variables, such as massaveraged position and kinetic energy, that are derived from the DFA simulations. This test is based on energy conservation considerations for simplified topographies. This allows us to verify the accuracy of the DFA simulations in terms of massaveraged runout. All the tests presented and used in what follows are implemented and available in AvaFrame (both data and helper functions). All results and figures can be reproduced using the code available on the AvaFrame GitHub repository (https://github.com/avaframe/AvaFrame/tree/theoryPaperCode, last access: 5 October 2023).
5.1 Flow variable tests
Before performing the abovementioned similarity solution and dam break tests, it is necessary to describe the quantities that are compared and the measures that are used to assess the convergence, accuracy or precision of the numerical model. Both the flow thickness (h) and thicknessintegrated flow momentum flux ($h\stackrel{\mathrm{\u203e}}{\mathit{u}}$) are used to compare the analytical solution to the simulation results. Two different deviation measures are used to quantify the deviation between a reference solution and the simulation result on a domain (one or twodimensional). The first is based on the ℒ_{max} norm (uniform norm), and the second is based on the Euclidean norm (ℒ_{2} norm). Let f_{num} be the numerical solution and f_{ref} the reference solution defined on a domain Ω. The local deviation is defined by $\mathcal{E}\left(x\right)={f}_{\text{num}}\left(x\right){f}_{\text{ref}}\left(x\right)$, and the global deviation is defined by the following terms.

The uniform norm (ℒ_{max}) measures the largest absolute value of the deviation ℰ on Ω:
$${\mathcal{L}}_{\text{max}}\left(\mathcal{E}\right)=\underset{x\in \mathrm{\Omega}}{sup}\left(\right\mathcal{E}\left(x\right)\left\right).$$This norm is applied to one or twodimensional results. It can also be normalized by dividing the uniform norm of the deviation by the uniform norm of the reference. In this case we refer to the relative deviation:
$$\mathcal{R}{\mathcal{L}}_{\text{max}}={\displaystyle \frac{{\mathcal{L}}_{\text{max}}\left(\mathcal{E}\right)}{{\mathcal{L}}_{\text{max}}\left({f}_{\text{ref}}\right)}}.$$ 
The Euclidean norm (ℒ_{2} norm) gives an overall measure of the deviations
$${\mathcal{L}}_{\mathrm{2}}\left(\mathcal{E}\right)=\underset{x\in \mathrm{\Omega}}{\int}{\left\left\mathcal{E}\left(x\right)\right\right}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}x.$$It is useful to normalize the norm of the deviation either by dividing with the norm of the reference solution:
$$\mathcal{R}{\mathcal{L}}_{\mathrm{2}}={\displaystyle \frac{{\mathcal{L}}_{\mathrm{2}}\left(\mathcal{E}\right)}{{\mathcal{L}}_{\mathrm{2}}\left({f}_{\text{ref}}\right)}}$$or by the measure of the interval (${\mathcal{L}}_{\mathrm{2}}\left(\mathrm{1}\right)={\int}_{x\in \mathrm{\Omega}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}x$):
$$\mathcal{R}\mathcal{S}{\mathcal{L}}_{\mathrm{2}}=\sqrt{{\displaystyle \frac{{\mathcal{L}}_{\mathrm{2}}\left(\mathcal{E}\right)}{{\mathcal{L}}_{\mathrm{2}}\left(\mathrm{1}\right)}}}.$$The first normalization approach will give a relative deviation, whereas the second will give an average deviation of f on Ω.
5.1.1 Similarity solution test
The similarity solution is one of the few cases where a semianalytic solution is available for solving the thicknessintegrated equations. This makes it very useful for validating the implementation of denseflow avalanche numerical methods (here com1DFA). In this problem, we consider an avalanche governed by a dryfriction law (Coulomb friction) flowing down an inclined plane. The release mass is initially distributed in an ellipse (diameters of length L_{x} and L_{y}) with a parabolic thickness shape (H in the middle). This mass is released at t=0 and flows down the inclined plane, as illustrated in Fig. 3a for the initial time step and at some later time t. This semianalytic solution can be derived under the following assumptions. First, the ratio $\mathit{\u03f5}=H/L$ of flow thickness over spatial extent in flow and crossflow direction should be smaller than 1. Second, the slope angle must be larger than the bed friction angle. Third, the solution is assumed to retain the symmetry properties of the initial configuration relative to the moving center of mass. The full description of the conditions and assumptions, as well as the derivation of the solution, is presented in detail by Hutter et al. (1993). The term semianalytic is used here because the method enables us to transform the PDE (partial differential equation) of the problem into an ODE (ordinary differential equation) using a similarity analysis method. The ODE can be solved with much less effort, e.g., using an explicit fourthorder Runge–Kutta scheme.
This test is implemented in the ana1Tests module of AvaFrame, which offers functions to compute the semianalytic solution, compare it to the output from the DFA computational module, and visualize the results.
5.1.2 Dam break test
The dam break test is the second test for which an analytical solution of the thicknessintegrated equations is known. In this test, we also consider an avalanche governed by a dryfriction law (Coulomb material), released from rest on an inclined plane (see Fig. 3b). In the case of a thicknessintegrated model as derived by Savage and Hutter (e.g., in Hutter et al., 1993), an analytical solution exists under the assumption of shallowness of the flow. Furthermore, the friction angle has to be smaller than the slope of the inclined plane. This test, in contrast to the similarity solution test, focuses on the very early stages of the flow and not on the evolution over time and lateral spreading. The derivation of the dam break solution is described in Faccanoni and Mangeney (2013) and corresponds to a Riemann problem. It has the following initial conditions:
The dam break assumes an invariance in the y direction, which is achieved using a wide enough domain in the y direction so that lateral effects can be neglected.
5.1.3 Results
DFA simulations are computed using the com1DFA module in AvaFrame varying the different numerical parameters listed in Table 1. The range of α values (called “aPPK” in the code and figures) is determined by the convergence criterion (Eq. 43). The SPH kernel radius r_{kernel} (called “sphKernelRadius” in the code and figures) is varied around 5 m, which is the raster cell size value currently used for operational hazard mapping in Austria. The tested values of n_{ppk0} (called “nPPK0” in the code and figures) and C_{time} (called “cMax” in the code and figures) are listed in Table 1. A large n_{ppk0} and small C_{time} lead to very long computation time, which makes these values unrealistic and impractical to use. Instructions on how to reproduce the results presented below are provided in the Supplement.
For both of the tests, the numerical schemes to apply friction and the method used to compute the SPH gradient are crucial for obtaining a proper starting and stopping behavior of the flow. Some intermediate developments showed that adding the friction with methods differing from as the one presented here and computing the SPH gradient without taking the slope inclination into account leads to unsatisfactory results. This is why the friction force is added as described in Sect. 4.5 and the SPH force is computed as described in Appendix B2. In what follows, artificial viscosity is added with C_{Lat}=10) in the similarity solution test, which stabilizes the solution without degrading the match with the semianalytic solution. For the dam break test, adding artificial viscosity has a negative impact on the solution and the following results were produced with no artificial viscosity (C_{Lat}=0).
Figure 4 shows an example where a DFA simulation is compared to the semianalytical solution of the similarity solution test case. The chosen parameters for this example are r_{kernel}=3 m, $\mathit{\alpha}=\mathrm{3}$, ${n}_{\text{ppk}}^{\mathrm{0}}=\mathrm{15}$ and C_{time}=0.02 s m^{−1} and correspond to the most accurate of the simulations presented here. The reader can find the results of the similarity solution test with the standard parameters in the Supplement. Figure 4a and b show the flow thickness and momentum profiles in and across flow direction after 20 s of flow. Figure 4d shows the evolution of the ℛℒ_{2} and ℛℒ_{max} deviation with time. The deviation at the initial time step (t=0 s) is rather high (this is related to the random process to initialize the particles in the simulation) and then quickly decreases after a few seconds of simulation due to the reorganization of the particles. The deviation then increases again as the numerical inaccuracies accumulate. When varying the numerical parameters in the DFA simulations (according to Eq. 42), the computed ℒ_{2} deviations between DFA results and the semianalytical solution decrease (see Fig. 5). In Fig. 6, the comparison between a DFA simulation and the analytical solution of the dam break test is shown for flow thickness, flow velocity and momentum at t=15 s (upper panel). Figure 6d shows a top view of the flow colored by flow thickness. This panel also shows the domain on which the deviation between analytical and numerical solution is computed. Figure 6e shows the relative deviations ℛℒ_{2} and ℛℒ_{max} on flow thickness and momentum. The same behavior as for the similarity solution test is observed regarding the time evolution of the deviation. Computation was done with r_{kernel}=3 m, $\mathit{\alpha}=\mathrm{3}$, ${n}_{\text{ppk}}^{\mathrm{0}}=\mathrm{15}$ and C_{time}=0.02 s m^{−1}. The reader can find the results of the similarity solution test with the standard parameters in the Supplement.
Results from both the similarity solution tests and the dam break test validate the convergence criterion from Ben Moussa and Vila (2000). Indeed, with an α exponent smaller than −1, decreasing the SPH kernel radius and varying the other parameters according to Eq. (42) leads to a decrease in the deviation, whereas for larger exponents, $\mathit{\alpha}=\mathrm{0.5}$ for both tests or α=0 for the dam break test (Fig. 7), the deviation decreases only slightly or even increases. Moreover, it is observed (not shown in the figure) that decreasing the time step (decreasing the C_{max} parameter) with all other parameters fixed leads to a decreasing deviation. Finally, for these two specific cases, DFA simulation results converge towards the semianalytical or analytical solution.
5.2 Energy and runout testing
The energy line test compares the results of the com1DFA simulation to a geometrical solution derived from the total energy of the system. Solely considering Coulomb friction, this solution is motivated by the first principle of energy conservation along a simplified topography. In this case, the friction force only depends on the slope angle. The analytical runout is the intersection of the path profile with the geometrical line (α line) defined by the friction angle α. From the geometrical line it is also possible to extract information about the massaveraged flow velocity at any time or position along the path profile (see example in Fig. 9). For the detailed theory of this test, please refer to Appendix A.
In this test, we use the α line to evaluate the DFA simulation. Computing the massaveraged path profile for the particles (each particle corresponding to a material point) in the simulation and comparing it to the α line allows us to compute four error indicators. Fig. A1 illustrates the concept.
The first three are related to the analytical runout point defined by the intersection between the α line and the massaveraged path profile. The last one is related to the velocity.

The horizontal distance between the analytical runout point and the end of the path profile defines the ${\mathit{\u03f5}}_{\mathrm{s}}={\stackrel{\mathrm{\u203e}}{s}}_{\mathit{\gamma}}{\stackrel{\mathrm{\u203e}}{s}}_{\mathit{\alpha}}$ error in meters.

The vertical distance between the analytical runout point and the end of the path profile defines the ${\mathit{\u03f5}}_{z}={\stackrel{\mathrm{\u203e}}{z}}_{\mathit{\gamma}}{\stackrel{\mathrm{\u203e}}{z}}_{\mathit{\alpha}}$ error in meters.

The runout angle difference between the α line angle and the DFA simulation runout line defines the ${\mathit{\u03f5}}_{\mathit{\alpha}}=\mathit{\gamma}\mathit{\alpha}$ runout angle error.

The rootmeansquare error (RMSE) between the α line and the DFA simulation energy points defines an error in the velocity altitude $\frac{\stackrel{\mathrm{\u203e}}{{v}^{\mathrm{2}}}}{\mathrm{2}g}$.
5.2.1 Limitations and applicability
It is essential to state where the assumptions of this test hold. One of the important hypotheses for the energy solution is that the inclination of the material point trajectory is equal to the slope angle of the surface, i.e., where $\mathrm{d}l=\frac{\mathrm{d}s}{\mathrm{cos}\mathit{\theta}}$. If this hypothesis fails, e.g., due to a particle trajectory deviating from the steepest slope direction, as illustrated in Fig. 8, it is impossible to derive the analytical energy solution. In the 3D case, the distance traveled by the particles reads $\left\right\mathrm{d}l\left\right=\frac{\mathrm{d}s}{\mathrm{cos}\mathit{\gamma}}$, where γ is the angle between the dl vector and the horizontal plane, which can differ from the slope angle θ (γ≤θ). Here the energy solution is not the solution to the problem and hence cannot be used as reference. In this case, it would not be possible to distinguish what deviation is caused by the numerical error or because of the hypothesis being violated.
The α line can be used to study the effect of terms such as curvature acceleration, artificial viscosity, or pressure gradients. For example, the curvature acceleration modifies the friction term, depending on topography curvature and particle velocity, leading to a mismatch between the energy solution and the DFA simulation. Figure 10 shows this curvature effect. The topography considered here is an inclined plane that smoothly transitions into a horizontal plane, meaning that curvature only occurs in the transition part. The energy line test for this case shows that there is added dissipation only in the transition part, seen by squares not following the α line. Once all particles have reached the horizontal plane, the squares follow the α line again (with a shift in s coordinate or, perhaps more intuitively, in the z coordinate).
Finally, the effects of the pressure force can be studied. For example, with this test it can be shown that adding the pressure forces does not influence the simulation runout (not shown here). This can be explained by the fact that pressure forces do not dissipate any energy and hence should not affect the energy balance. However, pressure forces lead to particle trajectories that do not necessarily follow the steepest direction, which means that the fundamental hypothesis illustrated in Fig. 8 is not satisfied.
5.2.2 Grid orientation effect
The energy line test previously described is also used to test whether the numerical method implemented in com1DFA performs independently of the grid orientation. Indeed we saw in Sect. 3.1 that com1DFA uses a regular grid to update some variables such as flow thickness or flow velocity. In order to quantify the effect of the grid orientation on the simulation results, we perform tests where grid orientation is changed while keeping the grid cell size and topography the same. In the Supplement, we show two examples where we consider a parabolic slope, i.e., the topography varies only in one direction, and a bowl shape, i.e., the topography with a rotational symmetry about its center. The main axis of the flow is not always aligned with the grid, and we provide three cases. First, a 0^{∘} case in which the slope is invariant in the y direction (main flow direction aligned with the grid). Second, a 225^{∘} case, with the slope being invariant in a direction angled 225^{∘} from the y direction (main flow direction aligned with the grid diagonal). Third, a 120^{∘} case meaning that the slope is invariant in a direction angled 120^{∘} from the y direction (main flow direction is neither aligned with the grid nor with the grid diagonal). For each of these test cases, two simulations are performed, i.e., with or without pressure gradients. The results of these tests and the instructions on how to reproduce them are provided in the Supplement. The results are very satisfying. The runout distances differ by only a few centimeters between the 0, 120, and 225^{∘} cases. This is of the same order of magnitude as the numerical errors computed with the energy line tests. This proves the rather low dependence of the simulation results on the grid orientation.
The presented com1DFA module within the opensource simulation framework AvaFrame is aimed at and developed for the simulation of avalanches with flow or a mixed form of movement. These may mostly be classified as A2B7C1DE7F7G1H7J1 according to the morphological classification code (International Association of Hydrological Sciences. International Commission on Snow and Ice, 1981). Effects of snow entrainment or additional sources of resistance can be included. The corresponding theoretical background is not described in this publication but is available in the AvaFrame online documentation (Oesterle et al., 2022).
The default setup of the module targets very large to extremely large avalanches of catastrophic intensity, corresponding to avalanches of size 4–5 of the European Avalanche Warning Service (EAWS) size classification or a 150year return period in the case of Austrian hazard mapping (Johannesson et al., 2009). Additional calibrated parameter sets for small and mediumsized avalanches and wetsnow avalanches (albeit experimental) are also provided within the framework. However, these are exclusively optimized on a subset of Austrian test cases. Other types of avalanches, such as powder snow avalanches or similar, are unsuitable applications for the com1DFA module.
Com1DFA is based on the slopenormal thickness integration of the mass and momentum conservation equations, which are solved with a mixed particle–grid method. Particles are used to track mass and to compute the pressure forces following a smooth particle hydrodynamics (SPH) approach adapted to the special case of steep terrain. An underlying Eulerian grid is used to compute the flow thickness from the distribution of the particles. Velocity is also interpolated from the particles to the grid in order to compute the stabilizing artificial viscosity term. Physical starting and stopping behavior is ensured by considering differences in the friction force, whether the material is flowing or at rest. A time–space criterion ensuring the convergence of the implemented numerical method is provided. Of course limitations and assumptions exist, one of the main ones being the assumption of a moderately curved surface. For currently used resolutions of topographies (about 5–10 m at the time of writing) terrain features remain relatively smooth, but this already requires forcing of particles to stay on the topography. For potential higherresolution applications in the future, this assumption needs to be reevaluated.
To verify the numerical implementation of com1DFA, we apply a series of tests separated into two categories. Flow variable tests, i.e., similarity solution and dam break tests, are used for checking the proper spatiotemporal evolution of flow thickness and velocity. Runout testing checks the accuracy of global variables such as centerofmass runout and kinetic energy. The tests show the validity of the chosen time–space criterion, as well as the accuracy and precision of the com1DFA numerical solution.
Note that the computational efficiency of the com1DFA module is not a topic in this article since simulations with the standard setup compute within seconds to minutes on current computer hardware. To achieve better computational efficiency, we implemented parallel execution of multiple serial simulations. However, topics like computational efficiency, further indepth testing, and application to real topographies will be treated in future publications.
Current and future potential improvements to com1DFA include topics such as the improvement of entrainment and detrainment (e.g., concerning forests), a more numerically sound representation of dams and walls, or a probabilistic approach to uncertainties. These topics will need both conceptual developments and numerical improvements, but some are already being tackled. Possible applications to other types of gravitational mass flows might also be an exciting development, but our current focus lies on snow avalanches. However, since AvaFrame is an opensource framework, we invite everyone to use the presented modules and welcome feedback and contributions, regardless of the topic.
The conservation of energy for a material point (block model) flowing downslope from point O to point B reads as follows (assuming only Coulomb friction):
where δE_{fric} is the energy dissipation due to friction, n the normal vector to the slope surface, and dl is the elementary vector on the path profile traveled by the material point between O and B. The vertical component of the normal vector reads ${\mathit{e}}_{z}\cdot \mathit{n}=\mathrm{cos}\mathit{\theta}$, where θ is the slope angle. Here, m represents the mass of the material point, g the gravity, μ=tan α the friction coefficient and friction angle, z the elevation, and v the velocity of the material point. Note that in the 2D case, $\mathrm{d}l=\frac{\mathrm{d}s}{\mathrm{cos}\mathit{\theta}}$ is only true if the inclination of the material point trajectory is equal to the slope inclination, i.e., the material point is flowing in the steepest slope direction. Now considering O as the origin position (s_{O}=0 and v_{O}=0) leads to the following simplification:
Speaking in terms of altitude, the energy conservation equation can be rewritten (the subscript B is dropped):
Considering a system of material points flowing down a slope with Coulomb friction, we can sum the previous equation (Eq. A3) of each material point after weighting it by its mass. This leads to the massaveraged energy conservation equation:
where the massaveraged value $\stackrel{\mathrm{\u203e}}{a}$ of a quantity a is (k indicates the points indices) as follows:
In this way, we can define the centerofmass path $(\stackrel{\mathrm{\u203e}}{x},\stackrel{\mathrm{\u203e}}{y})$ and the centerofmass path profile $(\stackrel{\mathrm{\u203e}}{s},\stackrel{\mathrm{\u203e}}{z})$. The massaveraged quantities also follow the same energy conservation law when expressed in terms of altitude. This result is illustrated in Fig. A1b and applies to both the material point equation (Eq. A3) and the massaveraged energy conservation equation (Eq. A4). The light blue line in Fig. A1b is obtained by evaluating the massaveraged energy conservation (Eq. A4) at the final time (t_{end}) and position ($({\stackrel{\mathrm{\u203e}}{s}}_{\text{end}}$, ${\stackrel{\mathrm{\u203e}}{z}}_{\text{end}})$), where v^{2}=0. This leads to the α line (also called energy line) equation:
The SPH method used in shallowwater equations is in most applications applied on a horizontal surface. The theoretical development on a horizontal plane is described in Appendix B1. The denseflow avalanche model described in this paper should be expressed on the bottom surface, which is not necessarily horizontal. The SPH gradient computation development is detailed in Appendix B2.
B1 Standard method
Let us start with the computation of the gradient of a scalar function $f:{\mathbb{R}}^{\mathrm{2}}\to \mathbb{R}$ on a horizontal plane. Let ${P}_{k}={\mathit{x}}_{k}=({x}_{k,\mathrm{1}},{x}_{k,\mathrm{2}})$ and ${Q}_{l}={\mathit{x}}_{l}=({x}_{l,\mathrm{1}},{x}_{l,\mathrm{2}})$ be two points in ℝ^{2} defined by their coordinates in the Cartesian coordinate system $({P}_{k},{\mathit{e}}_{\mathrm{1}},{\mathit{e}}_{\mathrm{2}})$. ${\mathit{r}}_{kl}={\mathit{x}}_{k}{\mathit{x}}_{l}$ is the vector from Q_{l} to P_{k}l. Now consider the kernel function W:
r_{0}∈ℝ is the smoothing kernel length (or radius). In the case of the spiky kernel, W reads as follows (2D case):
$\left\left{\mathit{r}}_{kl}\right\right=\left\left{\mathit{x}}_{k}{\mathit{x}}_{l}\right\right$ is the distance between particles k and l, and r_{0} is the smoothing length. Using the chain rule to express the gradient of W in the Cartesian coordinate system (x_{1},x_{2}) leads to
with
and
which leads to the following expression for the gradient:
The gradient of f is then simply
B2 The 2.5D SPH method
We now want to express a function f and its gradient on a curved surface and express this gradient in the threedimensional Cartesian coordinate system $({P}_{k},{\mathit{e}}_{\mathrm{1}},{\mathit{e}}_{\mathrm{2}},{\mathit{e}}_{\mathrm{3}})$. Let us consider a smooth surface 𝒮 and two points ${P}_{k}:{\mathit{x}}_{k}=({x}_{k,\mathrm{1}},{x}_{k,\mathrm{2}},{x}_{k,\mathrm{3}})$ and ${Q}_{l}:{\mathit{x}}_{l}=({x}_{l,\mathrm{1}},{x}_{l,\mathrm{2}},{x}_{l,\mathrm{3}})$ on 𝒮. We can define 𝒯𝒫_{k} the tangent plane to 𝒮 in P_{k}. If u_{k} is the (nonzero) velocity of the particle at P_{k}, it is possible to define the local orthonormal coordinate system $({P}_{k},{\mathit{V}}_{\mathrm{1}},{\mathit{V}}_{\mathrm{2}},{\mathit{V}}_{\mathrm{3}}=\mathit{n})$ with timedependent ${\mathit{V}}_{\mathrm{1}}=\frac{{\mathit{u}}_{k}\left(t\right)}{\left\left{\mathit{u}}_{k}\left(t\right)\right\right}$ and n the normal to 𝒮 at P_{k}. Locally, 𝒮 can be assimilated to 𝒯𝒫_{k} and Q_{l} to its projection ${Q}_{l}^{\prime}$ on 𝒯𝒫_{k} (see Fig. B1). Similarly, we assimilate a function $f:\mathcal{S}\subset {\mathbb{R}}^{\mathrm{3}}\to \mathbb{R}$ to a function ${f}^{\prime}:\mathcal{T}{\mathcal{P}}_{k}\subset {\mathbb{R}}^{\mathrm{3}}\to \mathbb{R}$ satisfying $f\left(Q\right)={f}^{\prime}\left({Q}^{\prime}\right)$.
The vector ${{\mathit{r}}^{\mathbf{\prime}}}_{kl}={\mathit{x}}_{k}{{\mathit{x}}^{\mathbf{\prime}}}_{l}$ from ${Q}_{l}^{\prime}$ to P_{k} lies in 𝒯𝒫_{k} and can be expressed in the local basis of 𝒯𝒫_{k}:
It is important to properly define f and its gradient as follows:
Indeed, since $({x}_{\mathrm{1}},{x}_{\mathrm{2}},{x}_{\mathrm{3}})$ lies in 𝒯𝒫_{k}, x_{3} is not independent of (x_{1},x_{2})
The target is the gradient of g^{′} in terms of the 𝒯𝒫_{k} variables (v_{1},v_{2}). Let us call this gradient ∇_{𝒯𝒫}. It is then possible to apply the Appendix B1 method to compute this gradient:
This leads to
This gradient can now be expressed in the Cartesian coordinate system. It is clear that the change in coordinate system was not needed:
Computing the gradient in the local coordinate system is, however, advantageous if the components (in flow direction or in cross flow direction) need to be treated differently.
The AvaFrame software is publicly available at https://github.com/avaframe/AvaFrame/ (last access: 5 October 2023). The associated online documentation is available at https://docs.avaframe.org/en/latest/ (last access: 23 October 2023). The code is released under the European Union Public license (EUPL version 1.2). The version of the code (including configuration files) and documentation for reproducing the results presented in this paper are available at https://github.com/avaframe/AvaFrame/tree/theoryPaperCode (last access: 5 October 2023) and https://docs.avaframe.org/en/theorypapercode/index.html (last access: 5 October 2023) (DOI: https://doi.org/10.5281/zenodo.7189007, Oesterle et al., 2022).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1670132023supplement.
AW, FO, and MT are the principal programmers and jointly developed the denseflow avalanche numerical model code. MT prepared the article and performed the model simulations to produce the results and figures of this article with contributions from AW. FO and JTF jointly supervised the project and contributed to the discussion and paper.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to acknowledge Peter Sampl, who developed the theory and source code for SamosAT which represented the starting point of the com1DFA module of AvaFrame. Thank you to Marie von Busse and Wolfgang Fellin for their help in developing the energy line test. We also would like to thank our anonymous reviewer and Dieter Issler for the constructive review that immensely helped improve the article's clarity and quality.
AvaFrame is supported by the Austrian Federal Ministry for Agriculture, Forestry, Regions and Water Management through cooperation between the Austrian Avalanche and Torrent Control (WLV) and the Austrian Research Centre for Forests (BFW). Additional support has been provided by the international cooperation project “AvaRange – Particle Tracking in Snow Avalanches” supported by the German Research Foundation (DFG, project no. 421446512) and the Austrian Science Fund (FWF, project no. I 4274N29).
This paper was edited by Ludovic Räss and reviewed by Dieter Issler and Martin Mergili.
Ata, R. and Soulaïmani, A.: A stabilized SPH method for inviscid shallow water flows, Int. J. Numer. Meth. Fl., 47, 139–159, https://doi.org/10.1002/fld.801, 2005. a
Baker, J. L., Barker, T., and Gray, J. M. N. T.: A twodimensional depthaveraged μ(I)rheology for dense granular avalanches, J. Fluid Mech., 787, 367–395, https://doi.org/10.1017/jfm.2015.684, 2016. a
Ben Moussa, B. and Vila, J. P.: Convergence of SPH method for scalar nonlinear conservation laws, SIAM J. Numer. Anal., 37, 863–887, https://doi.org/10.1137/S0036142996307119, 2000. a, b
Christen, M., Kowalski, J., and Bartelt, P.: RAMMS: Numerical simulation of dense snow avalanches in threedimensional terrain, Cold Reg. Sci. Technol., 63, 1–14, https://doi.org/10.1016/j.coldregions.2010.04.005, 2010. a, b, c
Faccanoni, G. and Mangeney, A.: Exact solution for granular flows, Int. J. Numer. Anal. Met., 37, 1408–1433, https://doi.org/10.1002/nag.2124, 2013. a, b
Fischer, J. T., Kowalski, J., and Pudasaini, S.: Topographic curvature effects in applied avalanche modeling, Cold Reg. Sci. Technol., 74–75, 21–30, https://doi.org/10.1016/j.coldregions.2012.01.005, 2012. a
Fischer, J.T., Fromm, R., Gauer, P., and Sovilla, B.: Evaluation of probabilistic snow avalanche simulation ensembles with Doppler radar observations, Cold Reg. Sci. Technol., 97, 151–158, https://doi.org/10.1016/j.coldregions.2013.09.011, 2014. a
Gauer, P.: Comparison of avalanche front velocity measurements and implications for avalanche models, Cold Reg. Sci. Technol., 97, 132–150, https://doi.org/10.1016/j.coldregions.2013.09.010, 2014. a
Granig, M., Sampl, P., Fischer, J., Kofler, A., and Joerg, P.: Adaption and further development of the numerical solution in the avalanche simulation model SamosAT, in: 13th Congress INTERPRAEVENT, Lucerne, Switzerland, 30 May–2 June 2016, 284–289, https://interpraevent2016.ch/wpcontent/ (last access: 23 October 2023), 2016. a, b
Gray, J. and Edwards, A.: A depthaveraged μ(I)rheology for shallow granular freesurface flows, J. Fluid Mech., 755, 297–329, https://doi.org/10.1017/jfm.2014.450, 2014. a, b, c
Gubler, H.: Measurements and modelling of snow avalanche speeds, in: Avalanches Formation, Movement and Effects, edited by: Salm, B. and Gubler, H., IAHS Publication, 162, 405–420, 1987. a
Harten, A. and Hyman, J. M.: Self adjusting grid methods for onedimensional hyperbolic conservation laws, J. Comput. Phys., 50, 235–269, https://doi.org/10.1016/00219991(83)900669, 1983. a
Harten, A., Lax, P. D., and Leer, B. v.: On upstream differencing and Godunovtype schemes for hyperbolic conservation laws, SIAM Rev., 25, 35–61, https://doi.org/10.1137/1025002, 1983. a
Hergarten, S. and Robl, J.: Modelling rapid mass movements using the shallow water equations in Cartesian coordinates, Nat. Hazards Earth Syst. Sci., 15, 671–685, https://doi.org/10.5194/nhess156712015, 2015. a
Hutter, C., Siegel, M., Savage, S., and Nohguchi, Y.: Twodimensional spreading of a granular avalanche down an inclined plane Part I. Theory, Acta Mech., 100, 37–68, https://doi.org/10.1007/BF01176861, 1993. a, b, c
Ihmsen, M., Orthmann, J., Solenthaler, B., Kolb, A., and Teschner, M.: SPH Fluids in Computer Graphics, in: Eurographics 2014 – State of the Art Reports, edited by: Lefebvre, S. and Spagnuolo, M., The Eurographics Association, https://doi.org/10.2312/egst.20141034, 2014. a
International Association of Hydrological Sciences. International Commission on Snow and Ice: Avalanche atlas, in: Illustrated international avalanche classification, Unesco, Paris, 265 pp., https://unesdoc.unesco.org/images/0004/000480/048004MB.pdf (last access: 5 October 2023), 1981. a
Barbolini, M., Domaas, U., Faug, T., Gauer, P., Hakonardottir, K. M., Harbitz, C. B., Issler, D., Johannesson, T., Lied, K., NaaimBouvet, M., and Rammer, L.: The design of avalanche protection dams, in: Recent practical and theoretical developments European Commission, edited by: Jóhannesson, T., Gauer, P., Issler, D., and Lied, K., Project Report: EUR 23339, European Commission, Directorate General for Research, 212 pp., 2009. a
Körner, H. J.: The energyline method in the mechanics of avalanches, J. Glaciol., 26, 501–505, https://doi.org/10.3189/s0022143000011023, 1980. a
LeVeque, R. J.: Numerical methods for conservation laws, Birkhäuser Basel, https://doi.org/10.1007/9783034886291, 1990. a
Li, X., Sovilla, B., Jiang, C., and Gaume, J.: Threedimensional and realscale modeling of flow regimes in dense snow avalanches, Landslides, 18, 3393–3406, https://doi.org/10.1007/s10346021016928, 2021. a
Liu, M. and Liu, G.: Smoothed Particle Hydrodynamics (SPH): an overview and recent developments, Arch. Comput. Method. E., 17, 25–76, https://doi.org/10.1007/s1183101090407, 2010. a, b, c, d
Luca, I., Tai, Y.C., and Kuo, C.Y.: Shallow Geophysical Mass Flows down Arbitrary Topography, Springer International Publishing, https://doi.org/10.1007/9783319026275, 2016. a
MangeneyCastelnau, A., Vilotte, J.P., Bristeau, M.O., Perthame, B., Bouchut, F., Simeoni, C., and Yerneni, S.: Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme, J. Geophys. Res.Sol. Ea., 108, 2527–2544, https://doi.org/10.1029/2002JB002024, 2003. a, b, c
Mergili, M., Fischer, J.T., Krenn, J., and Pudasaini, S. P.: r.avaflow v1, an advanced opensource computational framework for the propagation and interaction of twophase mass flows, Geosci. Model Dev., 10, 553–569, https://doi.org/10.5194/gmd105532017, 2017. a, b
Monaghan, J.: Smoothed particle hydrodynamics, Annu. Rev. Astron. Astr., 30, 543–574, https://doi.org/10.1146/annurev.aa.30.090192.002551, 1992. a
Monaghan, J. J.: Smoothed particle hydrodynamics, Rep. Prog. Phys., 68, 1703–1759, https://doi.org/10.1088/00344885/68/8/r01, 2005. a
Oesterle, F., Tonnel, M., Wirbel, A., and Fischer, J.T.: avaframe/AvaFrame: Version 1.3, Zenodo [code and data set], https://doi.org/10.5281/zenodo.7189007, 2022. a, b, c
Oesterle, F., Tonnel, M., Wirbel, A., and Fischer, J.T.: avaframe/AvaFrame: Version 1.6.1, Zenodo [code], https://doi.org/10.5281/zenodo.8319432, 2023.
Rauter, M. and Tuković, Z.: A finite area scheme for shallow granular flows on threedimensional surfaces, Comput. Fluids, 166, 184–199, https://doi.org/10.1016/j.compfluid.2018.02.017, 2018. a, b
Rauter, M., Kofler, A., Huber, A., and Fellin, W.: faSavageHutterFOAM 1.0: depthintegrated simulation of dense snow avalanches on natural terrain with OpenFOAM, Geosci. Model Dev., 11, 2923–2939, https://doi.org/10.5194/gmd1129232018, 2018. a
Salm, B. and Gubler, H.: Measurement and analysis of the motion of dense flow avalanches, Ann. Glaciol., 6, 26–34, https://doi.org/10.1017/s0260305500009939, 1985. a
Sampl, P.: SamosAT – Modelltheorie und Numerik, AVL List GMBH, Tech. Rep., Zenodo, https://doi.org/10.5281/zenodo.8007681, 2007. a
Sampl, P. and Zwinger, T.: Avalanche simulation with SAMOS, Ann. Glaciol., 38, 393–398, https://doi.org/10.3189/172756404781814780, 2004. a, b, c, d
Savage, S. B. and Hutter, K.: The motion of a finite mass of granular material down a rough incline, J. Fluid Mech., 199, 177–215, https://doi.org/10.1017/s0022112089000340, 1989. a
Stomakhin, A., Schroeder, C., Chai, L., Teran, J., and Selle, A.: A material point method for snow simulation, ACM T. Graphic., 32, 1–10, https://doi.org/10.1145/2461912.2461948, 2013. a
Voellmy, A.: Über die Zerstörungskraft von Lawinen, Schweizerische Bauzeitung, Sonderdruck aus dem 73. Jahrgang, 1–25, 1955. a
Zugliani, D. and Rosatti, G.: TRENT2D: An accurate numerical approach to the simulation of twodimensional dense snow avalanches in global coordinate systems, Cold Reg. Sci. Technol., 190, 103343, https://doi.org/10.1016/j.coldregions.2021.103343, 2021. a, b, c
 Abstract
 Introduction
 Mathematical model: from 3D equations to thicknessintegrated Lagrangian equations
 Numerical method: particle–grid approach
 Numerical implementation and solver
 Model verification
 Conclusions
 Appendix A: Theory: energy line test
 Appendix B: SPH pressure gradients
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Mathematical model: from 3D equations to thicknessintegrated Lagrangian equations
 Numerical method: particle–grid approach
 Numerical implementation and solver
 Model verification
 Conclusions
 Appendix A: Theory: energy line test
 Appendix B: SPH pressure gradients
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement