2.1 Conservation equations and depth integration

The presented method fundamentally relies on balance equations, in particular, the conservation of mass and momentum for fluids. The combination of these two equations is widely known as the Navier–Stokes equations (e.g. Ferziger and Peric, 2002) and can be written as

with the bulk density ρ and the bulk velocity u. (Note that it can also be defined for an individual phase with some modifications, see, for example, Rauter, 2021.) These flow fields are functions of time t and space $\mathit{x}=(x,y,z{)}^T$. The model Eqs. (2) and (3) describe their evolution from a known initial state $\mathit{u}\left(\mathrm{0},\mathit{x}\right)={\mathit{u}}_{\mathrm{0}}\left(\mathit{x}\right)$ and similar ρ₀(x), under the influence of boundary conditions. The divergence of the stress tensor T has the effect of diffusing momentum, and the volume force f represents additional forces, such as gravity.

Appropriate closure relations that express the stress tensor T as a function of the unknown flow fields yield a well-posed problem that can, in principle, be solved with numerical methods (Barker and Gray, 2017). However, even a well-posed problem is often not practically solvable from a computational perspective. Therefore, multiple simplifications have to be made to make problems of practical relevance accessible. Simplifications often come in the form of averaging over a certain time or over space to get rid of turbulent structures (Reynolds averaging; see, for example, Ferziger and Peric, 2002), to describe the average behaviour of multiple interpenetrating phases (phase averaging; e.g. Rauter, 2021), or to get rid of the vertical dimension (e.g. Savage and Hutter, 1989; Rauter and Tuković, 2018). The latter is referred to as depth averaging or depth integration and avoids the calculation of three-dimensional flow details. It yields mean values of, for example, density $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ and velocity $\stackrel{\mathrm{‾}}{\mathit{u}}$ along the depth.

In the simplest case, where the depth integration is aligned with a spatial axis (e.g. the z axis), the problem can be reduced from three $(x,y,z)$ to two dimensions (x,y). In this case, the depth-averaged value for an arbitrary field ψ is defined as

The newly introduced field $h(x,y,t)$ describes the flow depth, here in terms of the z coordinate of the top boundary of the integration, for a bottom boundary assumed to be aligned with z=0. The bottom and top boundaries are usually defined such that the mass flux through them is zero, meaning that they move with the vertical velocity of the flow at the respective position. The simplest example of such a model is the shallow-water equations (Barré de Saint-Venant, 1871). Defining the boundary in any other way will lead to additional source or sink terms, depending on the mass flux through the boundary (e.g. Pudasaini and Hutter, 2007). Examples would be any kind of entrainment and deposition fluxes.

Depth-integrated models are often considered synonymous with two-dimensional models. However, real avalanches and landslides travel along paths and surfaces in three-dimensional space. The three-dimensional nature of the terrain has to be reintroduced by modifying the two-dimensional model equations. Most often this is accomplished by abandoning Cartesian coordinate systems and Euclidean geometry, which was described in detail first by Savage and Hutter (1989, 1991) and extended by many others since then (e.g. Bouchut and Westdickenberg, 2004; Denlinger and Iverson, 2004; Pudasaini et al., 2005; Hergarten and Robl, 2015). This introduces various correction terms based on Christoffel formalism that are difficult to handle in complex models. In practice, simpler approximations are frequently employed (e.g. in RAMMS, see Fischer et al., 2012), leading to a disparity between theory and practical implementation. Notably, many of these developments happened in parallel and independently in the Russian avalanche dynamics community (Eglit et al., 2020).

An alternative to two-dimensional models with excessive curvature terms is the direct solution of the governing equations in three-dimensional space (Craster and Matar, 2009; Hagemeier et al., 2011; Rauter and Tuković, 2018). Depth integration is still compatible with this approach, and it can in principle be conducted in any direction pointing out of the surface. Yet in this work, depth integration is always conducted in the direction of the normal vector n^Γ to the flow surface Γ, as shown in Fig. 2. This has formally to be conducted in a surface-aligned coordinate system x^′–y^′–z^′:

The Jacobian matrix J, representing the transformation $\partial {\mathit{x}}^{\prime }/\partial \mathit{x}$, and its determinant det(J) take into account the curvature of the surface and its influence on the volume in a differential volume element of the flow (Bouchut et al., 2003). This effect is of the order of $h/R$ (Bouchut et al., 2003) with the mean curvature radius R and thus is small for mildly curved surfaces (R is large in comparison to the flow height h). As in most other models, the influence of the curvature on depth integration is ignored in this work.

Figure 2Depth integration reduces the full three-dimensional flow field ψ (dashed area) to an average flow field $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ (blue-filled area), which is assigned to a point x_b∈Γ.

Download

2.2 Surface partial differential equations

Depth integration in terms of Eq. (5) projects all three-dimensional flow fields on the surface Γ they are constrained by. The conservation equations can then be expressed as surface partial differential equations (SPDEs) that are defined on the surface Γ and include derivatives of various fields along it (Rauter and Tuković, 2018). These derivatives emerge from depth-integrating the nabla operator ∇ present in the Navier–Stokes equations formulated in a three-dimensional Cartesian reference frame. The presented framework differs from the classic approach of handling derivatives in two respects. Both are described in the following.

Depth integration has different effects on derivatives taken along the surface compared to those taken normal to the surface. While derivatives normal to the surface either vanish or appear as local source terms, derivatives along the surface remain in the system. A very common approach is to write the surface-aligned derivative as a two-dimensional derivative in local coordinates (e.g. in terms of $\partial x^{\prime }$, $\partial y^{\prime }$; Savage and Hutter, 1989). In the present framework, however, we express all entities in global Cartesian coordinates. For diffusive processes, this procedure gives rise to the Laplace–Beltrami operator (Dziuk and Elliott, 2013). However, this technique can also be adopted to other differential operators present in the system. The respective surface gradient and divergence operators ∇^Γ can be readily calculated in the numerical framework; see Sect. 2.3.

The depth integration of derivatives is then conducted in analogy to ordinary fields (see Eq. 5), in the surface aligned coordinate system x^′–y^′–z^′:

where x_b is a point on the bottom of the flow (and thus the flow surface Γ) and x_t the corresponding point on the free surface of the flow. The second term on the right-hand side of Eq. (6) represents an additional sink or source term, which arises if ψ is not zero at the bottom, x_b, or the top of the flow, x_t, such as due to entrainment or basal friction.

Further, our approach does not follow, for example, Savage and Hutter (1989) in separating the z component (determining, for example the pressure) from the x and y component (determining, for example, the velocity) in vectorial type balance laws such as the momentum equation. This would not be possible as our coordinate system is not aligned with the surface. We rather project the full three-dimensional equation onto the surface and the normal vector. The surface tangential projection of the surface gradient of a scalar is hence given by

and the surface normal projection as

and similar for vectors and divergence operators. The benefit of this approach is that the framework operates entirely in global Cartesian coordinates; thus, similar to an inertial frame, no fictitious centrifugal forces have to be considered. It follows that leading-order curvature effects are considered in the model by design, without explicitly expressing the curvature (Rauter and Tuković, 2018). For higher-order curvature effects, det(J) needs to be preserved during depth integration (Bouchut and Westdickenberg, 2004).

With these building blocks and some knowledge on how to transform one-dimensional shallow-flow models (e.g. Savage and Hutter, 1989; Parker et al., 1986), it is possible to extend nearly arbitrary depth-integrated flow models to complex terrain. In particular, ordinary depth-integrated flow models represent the surface–tangential momentum conservation equation and the flow depth equation. The two-dimensional ∇ operators have to be replaced with the surface–tangential ${\mathbf{\nabla }}_{\mathrm{s}}^{\mathrm{\Gamma }}$ operators. The surface–normal momentum conservation equation can be applied to replace the usually simplified expression for the basal pressure.

2.3 Finite area method

Partial differential equations, as well as their SPDE counterparts, are rarely solvable in an analytical sense, especially practical problems that represent real-world situations. Therefore, we rely on numerical approximations of SPDEs and the finite area method. This method is a variation of the finite volume method (see Ferziger and Peric, 2002; Jasak, 1996; Moukalled et al., 2016, for details) in N+1 dimensions, where N is the dimension of the control volumes. This means that for two-dimensional control volumes (i.e. surfaces), vectorial entities, such as normal vectors, velocities, and fluxes, will be three-dimensional. Similar to the conventional finite volume method, the Gaussian surface theorem (Tuković and Jasak, 2012) is applied and discretized by simplifying a control surface S as a flat, convex polygon S_i, as shown in Fig. 3. The expressions for the differential operators follow as

and

Index e refers to a discrete number of straight edges that form the polygon with surface S_i. ψ_e is the average value of the field ψ on the edge e, L_e its length, and m_e is the Γ–tangential and edge-normal outward-pointing vector. S_i, L_e, and m_e are purely geometrical properties that are defined during mesh generation. Values of fields on edges ψ_e, on the other hand, are interpolated from values of edge-adjacent cells, ψ_P and ψ_N. This introduces flux transport across cells and represents the flow of mass or information from one cell to neighbouring ones. The fluxes can then be associated in a linear system of equations that is solved with a suitable method.

Figure 3A finite area cell P and its neighbour N, used to calculate the approximation of surface derivatives in terms of the Gaussian surface theorem, integrating fluxes through cell edges e with length L_e and outward-pointing vector m_e.

Download

Discretization of non-gradient terms (e.g. the temporal derivative or any source term) is done in complete analogy to the finite volume method and obtained from integration over the control surface S_i. For details we refer to the large amount of excellent literature on the finite volume method (Ferziger and Peric, 2002; Jasak, 1996; Moukalled et al., 2016; LeVeque, 2002).

3.1 Friction in the dense-flow model

The term τ_Φ represents the depth-integrated divergence of the shear stress tensor and thus the constitutive model of the flowing mass. Assuming that the top boundary is stress-free and that surface–tangential derivatives of the deviatoric stress tensor are small, the only remaining entity is the basal friction. In this work, we will use the friction model presented by Rauter et al. (2016), which is closely related to the widely used Voellmy (1955) friction model. It is given as

with dry friction coefficient μ and turbulent friction coefficient χ. A wide range of alternative friction models can be found in the literature, and a number of them are implemented in the presented software.

3.2 Entrainment and deposition in the dense-flow model

$S_{\mathrm{\Phi }}^{\mathit{\varphi }}$ represents the sum of all volumetric source and sink terms of grains (e.g. erosion and entrainment of additional mass or its deposition), and $S_{\mathrm{\Phi }}^u$ represents its associated momentum. Dividing by the packing density in Eq. (12) simplifies handling of density changes in the different flow regimes. In the simplest case (e.g. laboratory experiments on a non-erodible bed), the source and sink terms are zero.

For snow avalanches and many other realistic gravitational mass flows, entrainment of erodible material along the avalanche path plays an important role. A popular entrainment model can be derived by comparing the dissipated energy in the mass flow with the energy required to mobilize the static material (Fischer et al., 2015):

with the specific erosion energy e_b as the single parameter. Here it is assumed that the packing density of the static layer is the same as in the dense flow ϕ_Φ.

Rauter and Köhler (2020) presented an extension to account for the deposition of flowing material, $S_{\mathrm{\Phi }\to \mathrm{\Sigma }}^{\mathit{\varphi }}$. This aspect is neglected in this work, and the flow height of the last time step is assumed to be the final deposition of the model.

The total flux term between the static layer and the flowing avalanche is determined as the difference between entrainment and deposition:

The related momentum source and sink terms are zero in the case of single layer flows, as both erodible material and deposited material are static.

The height (in surface-normal direction) of the static material on the topography can be tracked with an additional evolution equation:

with

again under the assumption that the static layer has the same packing density as the flowing avalanche ϕ_Φ. Tracking the thickness of the static layer allows the limitation of the available entrainable material; hence we are able to turn off entrainment if the erodible layer is depleted.

4.1 Friction in the suspension flow model

Similarly as in the dense-flow model, the term τ_Π represents the depth-integrated divergence of the shear stress tensor. If the particle fraction in the suspension is low, it can be treated as a simple fluid. The simplest wall friction model can be represented with constant friction coefficient c_D (Parker et al., 1986):

However, suspension flows are inherently turbulent, and there is strong evidence that turbulence models with more complex friction relations might be required (e.g. Parker et al., 1986). Nevertheless, we will use the simple model in this work, and we should keep in mind that the wall friction coefficient c_D is an empirical parameter that might require adaption to flow conditions. Further, it is assumed that all dissipative processes, such as inter-granular friction (e.g. Boyer et al., 2011), are included in this term. Considering the accuracy and uncertainties of the problems at hand, this seems to be a reasonable compromise. Alternative approaches are the turbulence model of Parker et al. (1986), a depth integration of the Einstein viscosity model (e.g. Boyer et al., 2011), and a more complex granular rheology (Boyer et al., 2011).

4.2 Ambient fluid entrainment in the suspension flow model

The volume of the suspension flow will grow due to entrainment of ambient fluid. It is assumed (Parker et al., 1986; Turner, 1986; Ancey, 2004) that ambient fluid entrainment depends solely on the bulk Richardson number, which is given as

In contrast to, for example, Parker et al. (1986), we use the effective surface-normal acceleration g_eff instead of the constant gravitational acceleration $\left|\mathit{g}\right|$ to account for the influence of centripetal forces on particles in the flow. Adjusting the Richardson number with the centripetal acceleration leads to an increased amount of ambient fluid entrainment if the flow runs over convex terrain and to a decreased amount if the flow runs over concave terrain.

There are various models for the relation between the Richardson number and the entrainment. Parker et al. (1986) use a simple, inverse proportional approach:

with the parameters α=0.00153 and Ri₀=0.204.

Turner (1986) provides an alternative formulation:

with the parameters Ri₀=0.8, α₁=10, and α₂=50. Various different parameters were suggested for this empirical relation (see, for example, Ancey, 2004).

Finally, Ancey (2004) suggests yet another relation in the form of an exponential function, here given in the form of Issler et al. (2018):

The parameter α₁ is supposed to be the only free parameter, with a value of 1.6 following Issler et al. (2018); however, due to different definitions of the entrainment rate, an additional parameter α₂ is required. In order to be of similar magnitude as the other air entrainment relations, α₂ has to be roughly 0.05. All relations are shown in Fig. 4.

Figure 4Comparison of the air entrainment functions, all depending solely on the Richardson number Ri.

Download

4.3 Grain entrainment and settlement in the suspension flow model

Suspension flows are, similar to dense flows, able to erode granular material from the bed. It is, in principle, possible to use the same entrainment relations as in the dense-flow model, but specialized entrainment relations have been proposed in the literature. An example of subaquatic turbidity currents, is given by Parker et al. (1986) as

with

the settling velocity,

the particle Reynolds number,

the viscosity of the ambient fluid ν_c; and two empirical parameters Z_m=13.2 and Z_c. The parameter Z_c was reported to be approximately 5. We found that a value of exactly 0.5 is required to reproduce the examples of Parker et al. (1986) in the examples shown in Sect. 7.2.

The settling of grains is given by Parker et al. (1986) as

with the settling velocity as given in Eq. (34) and the factor r₀ for the bottom value of the grain concentration,

As before, the total flux term follows as the difference between entrainment and deposition:

The momentum flux into the suspension due to ambient fluid and grain entrainment is zero. The volume occupied by entrained and deposited grains and the respective flux term in the evolution equation of the flow height h_Π are neglected at this point.

5.1 Cross-layer coupling

All fluxes of the two-layer model are described relatively well in the literature (see sections above), except for the mass flux from the dense-flow layer to the suspension layer, $S_{\mathrm{\Phi }\to \mathrm{\Pi }}^{\mathit{\varphi }}$, for which only few suggestions can be found (Eglit, 1998; Issler, 1998; Sampl and Zwinger, 2004; Bartelt et al., 2016). Existing relations do not conceptually fit into the presented framework, either due to missing granular mechanics (Eglit, 1998; Issler, 1998; Sampl and Zwinger, 2004) or due to their dependence on a specific dense-flow model (Bartelt et al., 2016). For the purpose of introducing this framework, we choose a simple relation, based on local flow fields of the dense flow.

We assume that the dense flow is composed of small and large particles with diameters d_Π and d_Φ, respectively. Uptake of particles into the suspension layer requires small particles to be made available by the dense layer that mostly consists of large particles (Bartelt et al., 2016) and the capability of the suspension layer to keep them suspended. The latter is already implemented into the model in the form of the settling model of the suspension flow $S_{\mathrm{\Pi }\to \mathrm{\Phi }}^{\mathit{\varphi }}$. This term depends on the particle Reynolds number Re_g, which is similar to the Stokes number and a good indicator of the flow regime.

The first step, making small particles available to the suspension, is assumed to be triggered by a fluidized flow that expands in volume, sucks in air, and increases the distance between particles. There are various hints on how this expression should look. At first it is useful to look at dimensionless properties in the dense flow. Besides the non-dimensional volumetric mass flux $S_{\mathrm{\Phi }\to \mathrm{\Pi }}^{\mathit{\varphi }}/\left|{\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{\Phi }}\right|$, these are the friction coefficient ${\mathit{\mu }}_{\mathrm{\Phi }}=\left|{\mathit{\tau }}_{\mathrm{\Phi }}\right|/p_{\mathrm{\Phi }}$, the packing density ϕ_Φ, and the inertial number (Forterre and Pouliquen, 2008; Rauter, 2021):

with the shear rate at the bottom of the dense flow, which is assumed to follow Bagnold (1954):

It is well established that μ_Φ and ϕ_Φ can be expressed as a function of only the inertial number I_Φ, and it is reasonable to assume that fluidization can be described the same way. This is further emphasized by the linear relationship between the packing density and the inertial number in the dense-flow regime (Forterre and Pouliquen, 2008). Finally, Rauter et al. (2016) found a specific relation between the shear rate ${\dot{\mathit{\gamma }}}_{\mathrm{\Phi }}$ and the pressure p_Φ in a granular kinetic theory model (Vescovi et al., 2013) at the point where fluidization suddenly occurs:

Comparing this relation to the expression for the inertial number, one can observe a striking resemblance; solely the exponent of the pressure is slightly lower in the relation of Rauter et al. (2016). This strongly indicates that the mass flux from the dense flow to the suspension can be expressed as a function of the inertial number only, starting at a minimum value I₀ and growing with a specified rate s_f from thereon:

The results of Rauter et al. (2016) suggest that the value of I₀ is close to 0.5, as at this point explosive fluidization starts to occur. The factor s_f is expected to be small, as the vertical velocity has to be substantially smaller than the flow velocity. This parameter can be optimized to yield the correct relation between dense flow and powder cloud.

In this model, small particles will be made available to the suspension when the dense-flow velocity is high or when the pressure is low (e.g. when an avalanche is running over a bump). If the small particles are sufficiently small, the suspension will be able to keep the particles suspended and a powder cloud will form. Otherwise, the particles will fall back to be reintegrated by the dense core, expressed by the deposition mass flux of the suspension layer, which is stronger for larger particles. The parameters for the suggested model are the small and large particle diameters d_Π and d_Φ, the minimum value of I at which fluidization occurs I₀, the particle density ρ_g, and the factor s_f. All parameters except for the latter are already used in the model or known otherwise.

Relation (48) finally completes the model and closes the system that will be solved numerically in the following. The model could be improved by tracking and limiting the availability of small particles or by making this property temperature-dependent.

6.1 Mesh generation from terrain data

The mesh generation follows the principles of Rauter et al. (2018). In a first step a triangulation of the terrain and a boundary of the surrounding volume is generated. A new tool for this task, called gridToSTL, was written entirely in C++ and without any external dependencies. The tool requires input in the form of a polygon that defines the simulation domain and the terrain data in the form of a raster file. In contrast to the previous version, the polygon can be any kind of closed and non-intersecting polygon with an arbitrary number of edges, either convex or concave. This enables flexibility in the choice of the simulation domain, which turned out to be especially useful to cover long and windy submarine canyons.

The finite volume mesh is generated from the triangulated surface with an arbitrary mesh generator. This toolchain can be applied to the depth-integrated models presented here and was also used for the full three-dimensional model presented by Rauter et al. (2022). In this study we used the mesh generator pMesh, while Rauter et al. (2022) used cartesianMesh, both from the cfMesh toolbox (Juretić, 2015). The finite area mesh is then generated on a dedicated surface of the finite volume mesh using the tool makeFaMesh, part of the OpenFOAM finite area module.

6.2 Mapping initial conditions

Initial conditions can be set with the tool releaseAreaMapping. In addition to the functionality of previous versions, this tool is now able to read shape files and grid files and map them directly onto finite area fields to be used by any solver. All input for the tool is read from a dictionary, where further references to shape and grid files can be listed. This tool enables efficient adaption to new scenarios.

6.3 Simulation run

Once the mesh and the initial conditions are defined, the solver of choice can be run. Currently there are three solvers available in the avalanche module, the dense-flow solver faSavageHutterFoam, the suspension flow solver faParkerFukushimaFoam, and the mixed-flow solver faTwoLayerAvalancheFoam (Fig. 5 shows faTwoLayerAvalancheFoam only, but it can be replaced with any of the other two models). Physical parameters are read from the file transportProperties, general simulation settings from controlDict, and numerical algorithms and parameters from the files faSolution and faSchemes. To run the solver in parallel, the tool decomposePar has to be run before the solver, and the tool reconstructPar has to be run after the solver. In the common OpenFOAM manner, all steps for a simulation are listed in a script file named Allrun, which can be launched by the user to automatically execute the pipeline outlined here. Another script, named Allclean, can be run to clean up the simulation directory.

6.4 Postprocessing and data export

The OpenFOAM architecture allows the execution of customized code, called function objects, in every simulation step. Various function objects are made available in the avalanche module. Most importantly, this includes function objects to export simulation results as either shape or raster files. The export as shape files can be done cellwise (one polygon for each computational cell), or the numerical data can be recombined to generate isolines that are written into the shape file. Function objects can be loaded by placing the respective entry into the control dictionary. As of version v2312, all solvers are able to run in a post-processing mode, in which old results are read from a hard disc and the function objects are executed. This allows the execution of function objects in a post-processing workflow without rerunning the whole simulation.

7.1 Dense-flow model

The dense-flow model has been applied to various cases in multiple studies. The interested reader is referred to Rauter and Tuković (2018) for lab scale simulations, Rauter et al. (2018) and Huber et al. (2018) for large-scale snow avalanche simulations, Rauter and Köhler (2020) for simulations with the deposition model, and to Shimizu (2022) for an application to pyroclastic flows.

7.2 Suspension flow model

Parker et al. (1986) simulate steady suspension flows on constantly inclined one-dimensional slopes with the model presented in Sect. 4. Four cases with uniform model parameters but different boundary conditions give a good overview of the behaviour of the model and a verification (as defined by Roache, 1997, as solving the equations right) of the presented implementation. The four simulations are conducted on one-dimensional slopes with a gradient of 5 %. The gravitational acceleration follows as $\mathit{g}=(\mathrm{0.49},\mathrm{0},-\mathrm{9.81}{)}^T$ m s⁻² (chosen to match the setup by Parker et al., 1986). The parameters suggest that the suspensions are composed of sediment in water on a scale of a small turbidity current.

Table 1Parameters for the small-scale simulations of Parker et al. (1986).

Download Print Version | Download XLSX

Material parameters for this setup are given in Table 1. The left boundary condition (at x=0) prescribes the inflow in terms of the height h_Π, velocity ${\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{\Pi }}$, and grain flux ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{\mathrm{\Pi }}=h_{\mathrm{\Pi }}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{\Pi }}{\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{\Pi }}$, in particular as shown in Table 2. All parameters are given normalized to reference values H=2 m, U=0.874 m s⁻¹, and Ψ=0.00828 m² s⁻¹.

Table 2Inlet boundary conditions for the small-scale simulations of Parker et al. (1986), simulating four scenarios of igniting or fading turbidity currents.

Download Print Version | Download XLSX

The right boundary condition is modelled as a zero gradient for all fields, mimicking an outlet boundary condition. For a basic verification of the novel implementation of the suspension model, the respective simulations are repeated and compared to the original results. We will evaluate the buoyancy assumption of Parker et al. (1986), as well as the formulation with the correct density given in here. The simulations are conducted in an unsteady manner until the flow reaches a steady state, comparable to the results reported by Parker et al. (1986). Figure 6 shows results for the four cases.

The first case, starting with a high velocity but low particle fraction increases its particle fraction quickly, as the high velocity is sufficient to erode and pick up sediment. The second case starts with a very high phase fraction, leading to a sudden ignition of the flow at $x/H=\mathrm{60}$. The height of the suspension stays low and even decreases, showing that a high phase fraction can keep the suspension concentrated at the bottom. The third and fourth case start with a low velocity and low particle phase fraction, respectively, and the suspension fades out quickly. The height of the flow increases in both cases where the flow fades out, indicating that the momentum of the flow is diffused over larger volumes of fluid. This is consistent with the expected scaling of fluid entrainment with the Richardson number.

It can be seen that results of Parker et al. (1986) are reproduced with only small deviations. The OpenFOAM solver yields sharper edges than the implementation of Parker et al. (1986), especially visible in Fig. 6b. This small difference is most likely attributable to the numerical solution method or the numerical resolution. The correction of the time derivative and advection term with $(\mathrm{1}+r{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{\Pi }})$ has only a minimal influence on the model results. This is reasonable, considering the low value for the buoyant density ratio r=1.65 in these cases. These simulations suggest that the model of Parker et al. (1986) is implemented correctly. However, this can not be seen as a validation (Roache, 1997) of the model.

Figure 6Numerical simulation of the four test cases presented by Parker et al. (1986) with OpenFOAM, with and without the buoyancy assumption (corrected and uncorrected, respectively). The results of Parker et al. (1986) are reproduced with good accuracy. The buoyancy assumption fits well to the conditions of these numerical experiments.

Download

7.3 Two-layer model

7.3.1 Synthetic tests and sensitivity study

In order to better understand the two-layer model, we will conduct tests on synthetic topographies. The topography is based on a parabola with a length of L=4000 m and a height of H=2000 m, with an additional flat runout area of 2000 m. The slope has a width of 2000 m, leading to a simulated region of $x=[-\mathrm{4000},\mathrm{2000}]$ m and $y=[-\mathrm{1000},\mathrm{1000}]$ m. In addition, the influence of topographic structures will be investigated, as terrain features often initialize the formation of suspension flows (e.g. powder snow avalanches). A bump in the surface is created by superposing the parabola with a Secans hyperbolicus $\mathrm{sech}\left(x\right)=\mathrm{2}/\left(\exp \left(x\right)+\exp (-x)\right)$ at $X_{\mathrm{p}}=-\mathrm{2700}$ m with a height of H_p=150 m and a length of L_p=200 m:

$$\begin{array}{}\text{(49)}& {\displaystyle }z=H_{\mathrm{p}}\mathrm{sech}\left({\displaystyle \frac{x-X_{\mathrm{p}}}{L_{\mathrm{p}}}}\right)+\left\{\begin{array}{ll}H{\left({\displaystyle \frac{x}{L}}\right)}^{\mathrm{2}}& \text{for }x<\mathrm{0},\\ \mathrm{0}& \text{otherwise,}\end{array}\right.\end{array}$$

inspired by the experiments of Viroulet et al. (2017). All boundaries are implemented as von Neumann (zero gradient) boundary conditions.

The release area (initial condition) of the slide was formed by a rectangle between $x=[-\mathrm{3900},-\mathrm{3500}]$ m and $y=[-\mathrm{500},\mathrm{500}]$ m and an initial dense-flow height of h_Φ=5 m within that square. All other flow fields are set to zero. The parameters, roughly corresponding to snow avalanches, are given in Table 3, if not mentioned otherwise. The value for the coupling factor s_f is varied, and the sensitivity of the model to this parameter is investigated. Entrainment and deposition to and from the static layer are not included in this section for simplicity. The simulations were run for 90 s.

Besides the flow thickness, velocity, and phase fraction, we can analyse the dynamic pressure, which is an important indicator of the destructive potential of the flow. It is defined as

$$\begin{array}{}\text{(50)}& {\displaystyle }{p}_{\mathrm{d}\mathrm{\Phi }}={\mathit{\rho }}_{\mathrm{\Phi }}|{\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{\Phi }}{|}^{\mathrm{2}}\end{array}$$

for the dense flow and as

$$\begin{array}{}\text{(51)}& {\displaystyle }{p}_{\mathrm{d}\mathrm{\Pi }}=\left({\mathit{\rho }}_{\mathrm{g}}{\mathit{\varphi }}_{\mathrm{\Pi }}+{\mathit{\rho }}_{\mathrm{c}}\left(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{\Pi }}\right)\right)|{\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{\Pi }}{|}^{\mathrm{2}}\end{array}$$

for the powder cloud (e.g. Jóhannesson et al., 2009). In particular we evaluate the dynamic peak pressure, which is defined as the maximum of the dynamic pressure at a fixed point over time. Important limits that are used in the definition of hazard zones (e.g. in Austria) are 1 kPa (yellow zone) and 10 kPa (red zone) (Jóhannesson et al., 2009). Notably, the shape factor should be applied to the dynamic pressure for consistency, increasing all simulated pressures by 25 %. However, this is neglected in order to be consistent with previous work and the definition of hazard zones.

Results for a simple parabola (without surface bump) are shown in Fig. 7 for three values of s_f (10⁻⁵, 10⁻⁴, 10⁻³). This set of simulations allows some valuable conclusions on the model and in particular the coupling model. All simulations start with a dense flow that eventually feeds the powder cloud. The feeding rate of the powder cloud varies strongly due to the variation of the parameter s_f.

For a low value of s_f, the dense flow is not able to generate a strong powder cloud with a considerable particle phase fraction and thus density. A suspension flow develops eventually; however, it consists almost entirely of air, without any ice particles. Basically, this can be seen as a layer of air that is dragged along by the dense flow. The velocity, dynamic pressure, and runout distance of this layer are correspondingly low. As shown before, the flow height of the suspension layer grows strongly for fading flows, indicating a strong diffusion of momentum.

Increasing the value for s_f up to 10⁻⁴ leads to higher phase fractions up to 0.004, roughly corresponding to a density of 4 kg m⁻³. Further increasing the value to 10⁻³ leads to phase fractions of up to 0.02 and densities of 20 kg m⁻³, however only for short periods. Notably, these are depth-averaged phase fractions and densities, and the respective values close to the surface might be considerably higher. The corresponding dynamic pressure of the powder cloud is still low, and only the simulation with the highest coupling factor s_f is able to generate a red zone that extends beyond the red zone of the dense flow. These results seem reasonable, considering the average slope gradient of 50 % and the absence of any topographic features that might enhance the feed of the powder cloud. More powerful powder snow avalanches can be expected on steeper slopes and on slopes with high topography variations (e.g. steep cliffs or rough terrain). Further simulations (not shown here) revealed that the powder cloud increases substantially with higher slope gradients.

Results for the slope with a bump are shown in Fig. 8. The model shows a high sensitivity to the terrain, and this case represents natural slopes with varying gradients better. All simulations create a considerable powder cloud with high phase fractions. The highest phase fraction is reached shortly after the top of the bump where the negative centrifugal forces are strongest and the basal pressure the lowest. The phase fraction reaches up to 0.05, roughly corresponding to a density of 50 kg m⁻³. A shock is formed at the bump in the suspension layer due to the high gradient in the phase fraction, leading to a considerable pressure gradient that decelerates the flow. In all simulations the dynamic powder cloud pressure exceeds 10 kPa, and the respective high pressure zone extends beyond the dense-flow runout. The 1 kPa zone of the powder cloud reaches considerable runouts beyond the dense flow.

Table 3Parameters for the two-layer model for synthetic cases on parabolas and for the Wolfsgruben and Eiskar avalanches.

Download Print Version | Download XLSX

Figure 7Numerical simulations of mixed-snow avalanches on a parabolic slope with the two-layer model. The parameter s_f was varied between 10⁻⁵ (a, d, g, j), 10⁻⁴ (b, e, h, k), and 10⁻³ (c, f, i, l). Panels (a)–(i) show the cross section in the middle of the slide. The slope is shown as the lower black line. The flow thickness h_Φ is shown as an offset of the surface magnified by a factor of 20; the flow thickness h_Π is shown above the dense flow magnified by a factor of 10. The powder cloud is coloured according to the phase fraction ϕ_Π. The red and yellow lines below the slope mark the regions of high dynamic peak pressure p_d>10 kPa and intermediate dynamic pressure p_d>1 kPa for the dense flow (top) and the powder cloud (bottom) respectively. Panels (j)–(l) show the regions of high and intermediate dynamic peak pressure (dashed: dense flow, continuous: powder cloud) from the top. One tick on the axis equals 1000 m.

Download

Figure 8Numerical simulations of snow avalanches on a parabolic slope with a bump with the two-layer model. Same as Fig. 7 but with a bump with height 150 m and length 200 m at $x=-\mathrm{2700}$ m.

Download

The results on synthetic terrain show a reasonable behaviour of the model, both in terms of parameterization and response to the terrain. The effect of the terrain is well visible and corresponds to the assumptions from which the model was derived. The sensitivity of the model to the parameter s_f is well pronounced, and this factor can be utilized to fit the model to real-world observations. A value between 10⁻⁴ and 10⁻⁵ seems reasonable for the parameter of the coupling model.

Finally, we use the synthetic cases to showcase the sensitivity of the model to the air entrainment. Figure 9 shows the simulation on the synthetic terrain with the three presented air entrainment models. The differences are small but noticeable. In particular, the entrainment is stronger with the model of Ancey (2004), however, which is just a question of parameterization. More importantly, the model of Turner (1986) shows a more pronounced flow head. The Richardson number is low in the head, and the relation of Turner (1986) predicts the strongest entrainment at low Richardson numbers (see Fig. 4). Generally, all relations appear reasonable and well in line with each other. We will continue with the entrainment model of Parker et al. (1986) from here on. Considering an optimization of air entrainment parameters to real events, it might be useful to apply the model of Ancey (2004) instead, as it provides the clearest parameterization.

Figure 9Numerical simulations of snow avalanches on a parabolic slope with a bump. Same as Fig. 8 but with a variation of the air entrainment model and a fixed parameter s_f.

Download

7.3.2 Real-case example: the 1988 Wolfsgruben avalanche

The 1988 Wolfsgruben avalanche represents an important event in Austria, as it was the trigger for many developments and used repeatedly as a benchmark. The event, or at least its dense core, was featured by Fischer et al. (2015) and Rauter et al. (2018). Here we revisit the event with the new two-layer model and include the powder cloud into the analysis. The avalanche is characterized by a channellized, steep slope with an angle of 30° that transitions fairly abruptly to the flat valley floor and the opposite slope.

The preprocessing and simulation setup follows Rauter et al. (2018) but with the novel tool chain and an extended simulation domain to cover the full runout of the powder cloud. The initial release areas of the avalanche and the erodible snow covers are the same, following the linear approach:

$$\begin{array}{}\text{(52)}& {\displaystyle }\begin{array}{rl}{h}_{\mathrm{\Sigma }}(x,y,z,t=\mathrm{0})& =\left(h_{\mathrm{\Sigma }}\left(z_{\mathrm{0}}\right)+{\displaystyle \frac{\partial h_{\mathrm{\Sigma }}}{\partial z}}\left(z-z_{\mathrm{0}}\right)\right)\\ & \cos \left(\mathit{\theta }\left(x,y,z\right)\right),\end{array}\end{array}$$

where z is the surface elevation and z₀ the elevation of a reference point with the base value h_Σ(z₀). The growth rate $\frac{\partial h_{\mathrm{\Sigma }}}{\partial z}$ defines the change with elevation from that point. θ is the angle between the gravitational acceleration and the surface-normal vector. For the 1988 Wolfsgruben avalanche, we use the snow cover parameters h_Σ(z₀)=1.61 m, z₀=1289 m, and $\frac{\partial h_{\mathrm{\Sigma }}}{\partial z}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$.

The model parameters are shown in Table 3. The dense-flow parameters have been optimized in a previous study (Fischer et al., 2015; Rauter et al., 2018), however for a model ignoring the powder cloud. The addition of the powder cloud could lead to different optimal dense-flow parameters. However, that is out of the scope of this work, and we assume that the previously used parameters fit sufficiently well. The suspension parameters are deduced from literature where possible (density, grain diameter). The coupling parameter $s_{\mathrm{f}}={\mathrm{10}}^{-\mathrm{5}}$ was found after running some simulations, starting from the values derived from the simulations on synthetic cases. A higher value leads to an unrealistically short dense-flow runout and a lower value to a severe underestimation of the suspension impact pressure. The friction coefficient c_D was chosen to be sufficiently large for the powder cloud to not completely decouple from the dense core. Apart from this effect, the simulation is rather insensitive to the friction coefficient c_D.

Figure 10Numerical simulation of the Wolfsgruben avalanche with the two-layer model. The first row (a–d) shows the height of the dense-flow layer; the second row (e–h) shows the height of the powder cloud layer. Each tick on the x and y axes corresponds to 500 m. Terrain data: Amt der Tiroler Landesregierung (AdTLR).

Figure 11(a) The dynamic peak pressure of the suspension layer (solid lines) and the dense-flow layer (dashed line). The yellow line marks the 1 kPa isoline and the red line the 10 kPa isoline. (b) The deposition of the avalanche (at t=180 s). Terrain data: AdTLR.

Table 4Simulated dynamic peak pressure at the location where damage was observed.

Download Print Version | Download XLSX

Four time steps of the simulation are shown in Fig. 10, displaying isolines of the dense-flow height h_Φ (a–d) and the suspension flow height h_Π (e–h). The avalanche starts as a dense flow and rapidly accelerates due to the steep release area (Fig. 10a). Shortly after the release a strong suspension layer is formed that further accelerates beyond the velocity of the dense-flow layer (Fig. 10b, f). After roughly 40–50 s the avalanche reaches the bottom of the valley (Fig. 10c, g). The powder cloud outruns the dense flow and hits the valley floor first. The dense flow is stopped quickly due to the high granular friction, while the powder cloud keeps running up on the opposite slope for approximately 50 m of elevation. Both flows experience a shock that increases the flow height in the valley floor drastically. The deposition (i.e. the dense-flow height in the last time step) reaches up to 15 m, however, which does not account for the difference between flow (≈200 kg m⁻³) and deposition density (≈600 kg m⁻³).

Results for the dense flow can be validated by a comparison with the deposition (see Fig. 11b), and they are similar to previous studies with the same model and the model SamosAT (Fischer et al., 2015; Rauter et al., 2018). Results for the powder cloud are more difficult to validate. Traces of the powder cloud that can be identified in the field are limited and not straightforward to interpret, as no clear deposition pattern emerges from suspended flows. Further, the respective deposition can hardly be related to the impact pressure and thus the destructive potential of the flow. Therefore, we compare the simulated dynamic pressure with observed building damages from the respective avalanche (see Fig. 11a). This includes not only the suspension layer but also the dense flow. An evaluation of the dynamic peak pressure and the deposition height at damaged objects is shown in Table 4.

The dense flow does not reach the two destroyed buildings (Point 1 and 2 in Fig. 11a and Table 4) and stops about 20 m short. Points 12 and 18 were only slightly damaged by the suspension flow in reality but severely hit by the dense flow in the simulation, showing that the simulation tends too strongly to the left side (viewed in the flow direction). The deposition height that was recorded at selected points (Table 4) is matched well, assuming a compaction of the avalanche by a factor of 3 after deposition.

The suspension layer shows a very limited zone of high dynamic pressure (>10 kPa) but an extended zone of intermediate dynamic pressure (1–10 kPa). The model predicts dynamic pressures of 1–4 kPa where balconies and roofs have been damaged and 1–3 kPa where windows have been destroyed. This corresponds well with engineering estimations of resistance capabilities of the respective parts: windows are assumed to break at 2–4 kPa and doors, walls, and roofs at 3–6 kPa (Sovilla et al., 2015). The dynamic pressure of the suspension layer at the destroyed buildings (Point 1 and 2 in Fig. 11) is not sufficient to destroy the respective brick structures (25–45 kPa). These high values strongly indicate that the dense flow or a saltation layer must be responsible for these high impact pressures (see pictures in Fischer et al., 2015). Therefore, we conclude that the simulated suspension layer reaches all observed traces of the powder cloud without covering the region where no traces could be observed.

7.3.3 Real-case example: the 2019 Eiskar avalanche

On 15 January 2019, the Eiskar avalanche occurred after intense snowfall and a rapid temperature drop (Oesterle, 2019). The Eiskar avalanche differs substantially from the Wolfsgruben avalanche regarding topography and thus provides a good supplement to that case. The avalanche was initiated by a slab on the right-hand side of the avalanche path (looking in the flow direction) that fell onto a larger snow field. From there, the avalanche slope continues with an inclination of approximately 25° for 1500 m until reaching a flatter slope of 10°. The dense-flow avalanche ran 1000 m on the flattening slope, and the powder flow exceeded the dense flow by another 500 m, reaching the village of Ramsau. The powder cloud destroyed a wooden building, damaged a hotel, and knocked over a bus. The dynamic pressure required for the damage was estimated at 1–3 kPa. Areal pictures were taken after the event, which allowed the estimation of the initial snow cover, the release area, and the deposition. The data were used to derive parameters for the snow cover function (Eq. 52), h_Σ(z₀)=1.60 m, z₀=1275 m, and $\frac{\partial h_{\mathrm{\Sigma }}}{\partial z}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$, to reach a snow cover thickness of approximately 2.7 m at an elevation of 2200 m (Oesterle, 2019). Other aspects of the simulation, such as the preparation of the terrain data, match the simulation of the Wolfsgruben avalanche.

The first simulation (not shown) was conducted with the same parameters as for the Wolfsgruben avalanche. However, these parameters lead to a severe underestimation of the powder cloud, running short by approximately 400 m. Simulations with the model SamosAT (Sampl and Zwinger, 2004) showed similar results with the standard parameters (Oesterle, 2019). Therefore, the friction coefficients and the coefficient for the suspension feed were adjusted (see Table 3) to reach an appropriate runout and dynamic pressure at the observed impacts.

Figure 12Numerical simulation of the Eiskar avalanche with the two-layer model. The first row (a–e) shows the height of the dense-flow layer; the second row (f–j) shows the height of the powder cloud layer. Each tick on the x and y axes corresponds to 500 m. Terrain data: Land Steiermark/GIS-Steiermark.

Figure 13(a) The dynamic peak pressure of the suspension layer (solid lines) and the dense-flow layer (dashed line). The yellow line marks the 1 kPa isoline and the red line the 10 kPa isoline. The black line marks the limit of the estimated 1 kPa isoline following observations. (b) The deposition of the avalanche (at t=200 s). The black polygons mark regions of dense-flow deposition (above) and powder cloud deposition (below). Terrain data: Land Steiermark/GIS-Steiermark.

Five time steps of the simulation are shown in Fig. 12 in terms of the dense-flow height h_Φ (a–e) and the suspension flow height h_Π (f–j). The collapsing slab (Fig. 12a) falls down the steep cliff onto a larger snow field where it can entrain additional snow. After around 30 s the avalanche reaches a second cliff, and a powder cloud starts to emerge (Fig. 12b, g). The suspension layer keeps growing substantially in the slope section with an inclination of 25° (Fig. 12c, h) and starts to detach when reaching the flatter slope of the 10° inclination. The suspension layer reaches the village of Ramsau after approximately 90 s (Fig. 12i, j), while the dense flow comes to a halt at the exit of the valley (Fig. 12d, e). Interestingly the dense flow is pushed towards the left by terrain features at the exit of the valley, while the suspension layer is essentially unaffected by these small obstacles.

The corresponding zones of dynamic pressure are shown in Fig. 13a. The 1 kPa isoline of the suspension layer extends far into the village. This fit was used as a benchmark to determine the optimal model parameters and thus matches observations well. The final deposition of the model is shown in Fig. 13b and compared to the observed deposition. The observations distinguish between suspension and dense-flow depositions (Oesterle, 2019), and the same can be done in the numerical model. The dense-flow layer leaves behind up to 10 m thick deposits (to be corrected by a factor of $\mathrm{1}/\mathrm{3}$ to match the deposition density) with sharp edges, while the suspension generates deposits with 0.1–0.2 m thickness (to be corrected as well) that fade out gradually. Both the position of the respective deposits and their shape match the observations.

Overall, the model is able to reproduce the observed flow traces, from the dynamic pressure to the varying snow deposits, in a single simulation. However, the model parameters had to be fitted to achieve these results. The friction parameters have to be substantially lower than in the Wolfsgruben case, and the coupling factor has to be an order of magnitude higher. This indicated that conditions were substantially different between the two cases or that the model does not cover some important processes.