MoT-PSA: a two-layer depth-averaged model for simulation of powder snow avalanches on 3-D terrain

Abstract

For snow-avalanche hazard mapping, one needs efficient tools that nevertheless capture the essential physical processes. The code MoT-PSA (Method of Transport – Powder Snow Avalanche) described here is based on the two-layer depth-averaged formulation for mixed snow avalanches developed by Eglit and co-workers in the 1980s but is extended to 3-D terrain and uses a fast numerical scheme based on the method of transport. Compared to previous works, we introduce novel formulations for the suspension and deposition of snow from the dense core. Snow cover and air entrainment are quantified with physics-based models. A sensitivity study of the model parameters on an idealized topography shows that both the dense core and the parameters of the powder snow cloud (PSC) governing particle suspension and settling significantly affect the dynamics. As expected, we observe that snow cover entrainment favours the formation of large PSCs with long runout. The powder-snow avalanche that occurred in Lom (Norway) on 27 February 2020 is back-calculated using MoT-PSA. With plausible parameter values, the model reproduces the dense core stopping at the gully's base and the dilute PSC travelling across the frozen lake for almost 1 km.

1. Introduction

Mixed snow avalanches (henceforth also called powder snow avalanches) are characterized by three distinct flow regimes (Sovilla and others, 2015): a dense core is formed when a snow slab is released and entrains the dense snow cover. Where collisions between snow particles dominate over enduring frictional contacts in parts of the flow, a fluidized flow regime is attained, which is more dilute and faster than the dense core. An upper suspension layer (also termed powder snow cloud, abbreviated as PSC in the following) forms by suspension of fine snow grains from the fluidized flow and may grow in size (up to 100‒200 m high) by entrainment of ambient air if its mass is sufficiently high. The PSC is mainly above the front and body of the avalanche but may detach from it and reach longer runout even on flat terrain and counter-slopes. Thus, extensive damage may be produced by the PSC itself, which warrants its implementation in numerical models for hazard mapping.

Because of their simplicity, most of the existing operational numerical avalanche models (e.g. Bartelt and others, 2017) only describe the dense core and neglect the formation and dynamics of the PSC. The use of such models is only adequate when modelling wet snow avalanches or small slab avalanches on relatively flat terrain where powerful PSCs cannot develop. The dynamics of pure PSCs (i.e. detached from the dense/fluidized core) were explored in the 1970s with centre-of-mass models endowed with dynamically changing height and length (Kulikovskiy and Sveshnikova, 1977). Eglit (1983) and Nazarov (1991) derived and coded depth-averaged equations in one dimension for a two-layer model, in which the bottom layer represents the dense core and the upper layer the PSC. They included snow cover and air entrainment, so that the density of the PSC can evolve. Later, Fukushima and Parker (1990) combined the centre-of-mass approach of Kulikovskiy and Sveshnikova with the four-equation model of turbidity currents by Parker and others (1986) to model non-Boussinesq PSCs. In addition to the volume, mass and momentum conservation equations already identified by Eglit (1983), Fukushima and Parker (1990) included an additional equation for the conservation of turbulent kinetic energy. However, they did not model the dense core. All these approaches were limited to 2-D terrain, in terms of a prescribed cross-section of the avalanche path which must be chosen by the modeller. Naaim and Gurer (1998) formulated a fully 3-D model of the suspension layer including erosion, sedimentation and a two-equation turbulence model, for which a 2-D depth-averaged dense-flow model of the Voellmy-type provides the boundary conditions through one-way coupling, that is, the effect of the suspension layer on the dense flow is neglected. Sampl and Zwinger (2004) developed a similar approach but solved the dense-flow and suspension-layer equations simultaneously to achieve two-way coupling. Both models are computationally very demanding. More recently, Bartelt and others (2016) proposed a two-layer depth-averaged model for powder snow avalanches, whose underlying physical assumptions are, however, debated (Issler and others, 2018).

In this paper, we present a new two-layer depth-averaged code – MoT-PSA (Method of Transport – Powder Snow Avalanche) – to model mixed snow avalanches. The aforementioned regimes in mixed snow avalanches have been observed to produce quite sharp and distinguishable transitions in radargrams, impact pressures, velocities, densities (Sovilla and others, 2015) and deposits (Issler and others, 2020). This justifies modelling the structure of a mixed snow avalanche with distinct layers. For the sake of simplicity, we provisionally model the dense core and fluidized layer as a basal layer with constant density. An upper PSC layer may then form, for which formulas describing its interaction with the basal core and the surrounding air must be specified. The longitudinal extent of the non-suspended part of snow avalanches is typically two orders of magnitude larger than their thickness, which justifies using a depth-averaged approach to reduce the computational cost. The length–to–height ratio is only O(10) for the PSC so that depth-averaging is harder to justify, but the gain in computational efficiency of a code outweighs the loss of accuracy in hazard-mapping applications, as experience with a 1-D depth-averaged code (Issler, 1998) has shown. Moreover, variations of avalanche variables like velocity, density and pressure along the depth may be captured partly through adequate parametrization, although complex 3D phenomena related to turbulence have to be neglected in a depth-averaged approach. MoT-PSA was initially designed to extend Eglit's (1983) and Nazarov's (1991) two-layer model to 3-D terrain (Issler, 2023). It was therefore implemented within the MoT framework, already developed for dense snow avalanches (MoT-Voellmy, Issler, in preparation). Along the way, some of Eglit's closure assumptions concerning the mass exchange rates between the layers were modified, however, to account explicitly for the shear strength in the snow cover and the dense layer, to include deposition, and to make use of more recent experiments on the entrainment of ambient fluid in density currents. In the first part of this paper, we describe the mathematical model of MoT-PSA. In the second part, we perform a sensitivity analysis on an idealized parabolic topography to derive the influence of the model parameters on the flow dynamics. We then back-calculate a powder-snow avalanche that occurred in Norway to verify that the salient features of powder snow avalanches are captured by the model.

2. MoT-PSA model equations

Two depth-averaged flow layers are used to model mixed snow avalanches (Fig. 1). The bottom layer (denoted by index 1) represents the dense core, while the upper layer (denoted by index 2) models the PSC. The fluidized layer (Schaerer and Salway, 1980; Issler and others, 1996, 2020; Sovilla and others, 2015) – earlier called ‘light flow’ or ‘saltation layer’ – is not modelled explicitly but embedded within the dense layer. An erodible snow cover (denoted by index 0 in the following) is also considered. The index a refers to the ambient air.

Figure 1. Cross-section of a two-layer powder-snow avalanche. The bottom layer (1) represents the dense core (and the fluidized layer), and the upper layer (2) represents the powder snow cloud, which entrains air (a). The two layers may entrain the snow cover (0). The primitive variables are indicated with blue vectors and within the shaded white boxes. The boundary stresses are indicated with black vectors, and the volumetric exchange rates are denoted with red arrows.

In the mathematical model, the avalanche flows over a general 3-D topography Σ described by a function Z(X, Y) in a global Euclidian coordinate system, where X and Y are in the horizontal plane. On Σ, we define a curvilinear coordinate system, where the x-coordinate lines project vertically on the X-coordinate lines of the Euclidian system, and analogously for the y-coordinate lines. At each point on Σ, a z-coordinate line normal to Σ is defined, with z = 0 on Σ. The metric tensor on Σ is given by

(1)$$G = \left({\matrix{ {1 + {( {\partial_XZ} ) }^2} & {\partial_XZ\;\partial_YZ} \cr {\partial_XZ\;\partial_YZ} & {1 + {( {\partial_YZ} ) }^2} \cr } } \right).$$

It is used to calculate the scalar products of vectors in this non-orthogonal coordinate system and to derive the curvature (Issler, in preparation).

2.1. Conservation equations

Particles in the suspension layer are typically so small that they have the density of ice, ρ _ice = 917 kg m⁻³. Snow clods in the dense/fluidized layer have considerably smaller density, ρ _s, in the range 200–600 kg m⁻³ (McClung and Schaerer, 1985) but may be embedded in a matrix of snow grains (Issler and others, 2020). The model assumes the components of the granular mixture – interstitial as well as ambient air and snow particles or snow grains – to be incompressible. The bulk density of the mixture can nevertheless change because the volumetric concentration of the snow particles need not be constant.

Incompressibility of the constituents implies that not only mass and momentum but also the volume is conserved in the flow. Consequently, the density of at least some layers must be variable if they exchange mass between each other. In the PSC, particle concentration may change by orders of magnitude so that it is imperative to solve the volume, mass and momentum conservation equations of layer 2. With this, the ambient air density, ρ _a, can consistently be held constant. For simplicity and efficiency, MoT-PSA assumes the density of the snow cover, ρ ₀ to be constant, and the dense/fluidized layer, ρ ₁, to also be constant. This simplification comes, however, at the price of volume or mass not being strictly conserved in some situations: If, say, the snow cover with density ρ ₀ is eroded to a depth Δh ₀ by the dense-flow layer with density ρ ₁ and flow depth h ₁, ρ ₀ remains the same but ρ ₁ changes to $\rho _1^{{\prime}} = ( {\rho_1h_1 + \rho_0{\rm \Delta }h_0}) / ({ h_1 + {\rm \Delta }h_0} ) =$ $\rho _1 \left[ 1 + ( {\rho_0/\rho_1} ) \cdot ( {{\rm \Delta }h_0/h_1 )} \right] / ({1 + {\rm \Delta }h_0/h_1} ) .$ In the situation where the snow cover and the dense flow have similar density (ρ ₁ ≈ ρ ₀), the density of layer 1 does not change after erosion ($ \rho_1^{\prime}$ ≈ ρ ₁), and mass and volume are hence exactly conserved. However, with typical values for the dense core, ρ ₀ ≈ 0.5ρ ₁, Δh ₀ ≈ 0.3h ₁, one obtains $ \rho_1^{\prime}$ ≈ 0.85ρ ₁; similarly, for a dilute, fluidized flow with ρ ₀ ≈ 3ρ ₁, $ \rho_1^{\prime}$ ≈ 1.5ρ ₁ results: Erosion has the effect of changing the density of layer 1. Conversely, if one enforces the density of layer 1 to remain constant ($ \rho_1^{\prime}$ = ρ ₁), the flow depth is distorted, and mass not strictly conserved. Physically, the most reasonable approximation would be to enforce conservation of ice mass and accept that layer 1 expels some air into the ambient or ingests air from it as needed.

MoT-PSA solves conservation equations for the following depth-integrated variables: mass hold-ups ρ ₀h ₀(x, y, t), ρ ₁h ₁(x, y, t), ρ ₂h ₂(x, y, t), momentum hold-ups ρ ₁h ₁u₁(x, y, t) = (ρ ₁h ₁u ₁, ρ ₁h ₁v ₁) and ρ ₂h ₂u₂(x, y, t) = (ρ ₂h ₂u ₂, ρ ₂h ₂v ₂), and bulk density of the PSC, ρ ₂(x, y, t). The density of the snow cover, ρ ₀, and of the dense core, ρ ₁, are instead assumed to be constant.

In the presence of an erodible bed, the evolution of the bed depth (h ₀) can be written as an Exner equation:

(2)$$\partial _{\rm t}h_0 = {-}E_{01}-E_{02} + D_{10} + D_{20\;}, \;$$

where E ₀₁ and E ₀₂ represent the entrainment rates (m s^–1) of the bed into layer 1 and layer 2, respectively, D ₁₀ represents the deposition rate of the dense layer, and D ₂₀ represents the settling rate of snow particles from the PSC onto the snow cover. If ρ ₀ is assumed constant, Eqn (2) also describes mass conservation (snow particles and air) within the snow cover.

The volume and exact mass conservation equations for the dense layer are written as

(3)$$\partial _{\rm t}h_1 + {\boldsymbol \nabla }_\parallel{\cdot} ( {h_1{\boldsymbol u}_1} ) = E_{01}-S_{12} + D_{21}-D_{10}, \;$$

(4)$$\partial _{\rm t}( {\rho_1h_1} ) + {\boldsymbol \nabla }_\parallel{\cdot} ( {\rho_1h_1{\boldsymbol u}_1} ) = \rho _0E_{01}-\rho _1S_{12} + \rho _{\rm s}D_{21}-\rho _{\rm d}D_{10}, \;$$

where S ₁₂ is the suspension rate from the dense layer into the PSC and D ₂₁ is the settling rate of snow particles from the PSC onto the dense layer. ρ _d is the density of the deposited dense core. ${\boldsymbol \nabla }_\parallel$ indicates the slope-parallel (${\parallel}$) component of the gradient operator.

Equations (3) and (4) express the volume and mass conservation for a dense core with variable density. However, in our model, the dense core is simplified to a fluid with constant density ρ ₁. Hence, Eqn (3) is not considered. Instead, using the condition ρ ₁ = const., the mass conservation Eqn (4) simplifies to

(5)$$\partial _{\rm t}h_1 + {\boldsymbol \nabla }_\parallel{\cdot} ( {h_1{\boldsymbol u}_1} ) = \displaystyle{{\rho _0} \over {\rho _1}}\;E_{01}-S_{12} + \displaystyle{{\rho _{\rm s}} \over {\rho _1}}D_{21}-\displaystyle{{\rho _{\rm d}} \over {\rho _1}}D_{10}.$$

The volumes exchanged between layers 0 and 1 differ in the two directions (compare the terms on the right-hand side of Eqns (2) and (5) respectively), while the mass exchanged between the two layers is conserved. This is a consequence of the assumption ρ ₁ = const. When the snow cover is entrained, it loses a volume E ₀₁ per unit time and area, while the dense core gains a volume ρ ₀/ρ ₁ E ₀₁. For instance, in the case ρ ₁ < ρ ₀, when entrained, the snow cover must expand and ingest ambient air to attain the (constant) density of the dense core. Similarly, during the deposition process, the bed gains a volume D ₁₀ per unit time and area, while the flow loses a volume ρ _d/ρ ₁ D ₁₀ per unit time and area. The situation where ρ _d > ρ ₁ physically corresponds to the dense core compressing and sintering to form the final deposit, as is observed in nature (Issler and others, 2020).

The momentum conservation for the dense layer is written as

(6)$$\eqalign{&\partial _{\rm t}( \rho _1 h_1 {\boldsymbol u}_1) + {\boldsymbol \nabla }_\parallel{\cdot} ( {\rho_1h_1{\boldsymbol u}_1{\boldsymbol u}_1} ) \cr &\quad = \rho _1{\boldsymbol g}_\parallel h_1-{\boldsymbol \nabla }_\parallel \left({\displaystyle{1 \over 2}\rho_1g_{z, 1}h_1^2 } \right) - h_1{\boldsymbol \nabla }_\parallel \left[ {( {\rho_2-\rho_{\rm a}} ) g_{z, 2}h_2} \right] \cr & \quad \quad - {\boldsymbol \tau }_{01} + {\boldsymbol \;}{\boldsymbol \tau }_{21\parallel } -\rho _1S_{12}{\boldsymbol u}_1 + \rho _{\rm s}D_{21}{\boldsymbol u}_{2{\rm b}}-\rho _{\rm d}D_{10}{\boldsymbol u}_1, \;} $$

where ${\boldsymbol g}_\parallel$ is the slope-parallel component of the gravitational acceleration. The second term on the right-hand side of Eqn (6) represents the longitudinal pressure gradient that originates from the dense core. The third term on the right-hand side represents the longitudinal pressure gradient that originates from the superposition pressure from the PSC, cf. Eqn (18). To evaluate these longitudinal pressures, we assume for simplicity a hydrostatic stress state and therefore neglect passive and active stress states, which may arise in the avalanche during compressional and dilative motion. Introducing a variable passive/active earth pressure coefficient has been shown to produce more realistic shapes of the avalanche deposit (Gray and others, 1999), but was for simplicity not implemented in the model at this stage.

Introducing the hypothesis ρ ₁ = const., Eqn (6) reduces to

(7)$$\eqalign{&\partial _{\rm t}( h_1{\boldsymbol u}_1) + {\boldsymbol \nabla }_\parallel{\cdot} ( {h_1{\boldsymbol u}_1{\boldsymbol u}_1} ) \cr & \quad\quad= {\boldsymbol g}_\parallel h_1-{\boldsymbol \nabla }_\parallel \left({\displaystyle{1 \over 2}g_{z, 1}h_1^2 } \right) - \displaystyle{{h_1} \over {\rho _1}}{\boldsymbol \nabla }_\parallel \left[ {( {\rho_2-\rho_{\rm a}} ) \,g_{z, 2}h_2} \right] \cr & \quad\quad\quad- \displaystyle{{{\boldsymbol \tau }_{01}} \over {\rho _1}} + \displaystyle{{{}{\boldsymbol \tau }_{21\parallel }} \over {\rho_1}}-S_{12}{\boldsymbol u}_1 + \displaystyle{{\rho_{\rm s}} \over {\rho _1}}D_{21}{\boldsymbol u}_{2{\rm b}} - \displaystyle{{\rho _{\rm d}} \over {\rho_1}}D_{10}{\boldsymbol u}_1.} $$

In the following, we will adopt this hypothesis and hence compute the dynamics of the dense core using the conservation Eqns (5) and (7). However, the code already implements the conservation Eqns (3), (4) and (6), which allow modelling a dense core with variable density if one adds an explicit equation for the density or an evolution equation specifying its rate of change in terms of the other field variables.

The shear stress ${\boldsymbol \tau }_{01}$ acts at the base of the dense layer, while ${\boldsymbol \tau }_{21\parallel }$ acts between layers 1 and 2, which are projected along the flow direction (${\parallel}$). The term g _z,i is the component of the gravitational acceleration normal to the topography, inclusive of curvature effects:

(8)$$g_{z, i} = \{ {g\cos \theta + \kappa_{\boldsymbol u}{\Vert {{\boldsymbol u}_i} \Vert }^2} \} , \;$$

where i = 1, 2 and {} are the Macaulay brackets to denote the ramp function defined as $\{ \eta \} = 0\;( {\forall \eta < 0} ) , \;\;\eta \;( {\forall \eta \ge 0} )$. g is the acceleration due to gravity, θ is the slope angle of the topography, and κ _u is the terrain curvature in the flow direction, which is calculated with the help of the second fundamental form of the surface (Issler, in preparation) and assumed to be the same for the two layers. g _z,i may differ for the two layers if these have different velocities.

The volume, mass and momentum conservation equations for the PSC are given as

(9)$$\partial _{\rm t}h_2 + {\boldsymbol \nabla }_\parallel{\cdot} ( {h_2{\boldsymbol u}_2} ) = E_{02} + S_{12}-D_{21}-D_{20} + E_{{\rm a}2}, \;$$

(10)$$\begin{align}\partial _{\rm t}( \rho _2h_2) + {\boldsymbol \nabla }_\parallel{\cdot} ( {\,f_{\rho {\rm u}}\rho_2h_2{\boldsymbol u}_2} ) = & \rho _0E_{02} + \rho _1S_{12}-\rho _{\rm s}D_{21}-\rho _{\rm s}D_{20} \\ & + \rho _{\rm a}E_{{\rm a}2}, \;\end{align}$$

(11)$$\eqalign{\partial _{\rm t}( \,f_{\rho {\rm u}}\rho _2 & h_2{\boldsymbol u}_2) + {\boldsymbol \nabla }_\parallel{\cdot} ( {\,f_{\rho {\rm uu}}\rho_2h_2{\boldsymbol u}_2{\boldsymbol u}_2} ) \cr &= ( \rho _2-\rho _{\rm a}) {\boldsymbol g}_\parallel h_2 - {\boldsymbol \nabla }_\parallel \left[ \! {\,f_{{\rm c}\zeta }( \rho_2-\rho_{\rm a}) g_{z, 2}h_2^2 } \right] \cr &\quad -{\boldsymbol \tau }_{02} -{\boldsymbol \tau }_{12\parallel {\boldsymbol \;}}-{\boldsymbol \tau }_{{\rm a}2\parallel }-{\boldsymbol \;}{\boldsymbol p}_{{\rm n}\parallel } \cr & \quad+ \rho_1S_{12}{\boldsymbol u}_1 - \rho_{\rm s}D_{21}{\boldsymbol u}_{2{\rm b}} - \rho_{\rm s}D_{20}{\boldsymbol u}_{2{\rm b}}, \;} $$

where $ {\boldsymbol \tau}_{02}$ is the shear stress acting where the base of the PSC directly touches the snow cover; $ {\boldsymbol \tau}_{a2}$ is the shear stress exerted by the surrounding air on the top surface of the PSC; and p_n is the normal stress on the top surface of the PSC due to the stagnation pressure from the surrounding air. Both ${\boldsymbol \tau }_{{\rm a}2\parallel }$ and p_n|| are projected along the flow direction. E _a2 is the air entrainment rate.

The density and velocity profiles of the dense core are assumed to be uniform along the direction normal to the terrain, parametrized by the normalized depth coordinate ζ ₁ = (z − h ₀)/h ₁, that is, ρ _1*(ζ ₁) = ρ ₁ and u_1*(ζ ₁) = u₁ (the subscript * is used to indicate that the variable is evaluated at a certain normalized depth ζ, while the subscript b will later be used to indicate that the variable is evaluated at the base of the layer). In contrast, the snow concentration and velocities of the PSC vary strongly along the normalized flow depth ζ ₂ = (z − h ₀ − h ₁)/h ₂. The density profile of the PSC is assumed as

(12)$$\rho _{2\ast }( {\zeta_2} ) = \rho _{\rm a} + f_{\rm c}( {\zeta_2} ) c\, ( {\rho_{\rm s}-\rho_{\rm a}} ) = \rho _{\rm a} + f_{\rm c}( {\zeta_2} ) ( {\rho_2-\rho_{\rm a}} ) , \;$$

where c = (ρ ₂ − ρ _a)/(ρ _s − ρ _a) is the depth-averaged volumetric concentration of snow in the PSC, and ρ _s is the density of the bulk snow. The shape function of the concentration is modelled as a generic parabolic function, $f_{\rm c}( {\zeta_2} ) = c_0 + c_1\zeta _2 + c_2\zeta _2^2$, where the coefficients c ₀, c ₁, c ₂ may be chosen based on laboratory experiments (e.g. Hermann and Hutter, 1991) and must respect the condition $\int_0^1 {f_{\rm c}( {\zeta_2} ) {\rm d}\zeta _2 = 1}$. Similarly, the velocity profile is assumed as

(13)$${\boldsymbol u}_{2\ast }( {\zeta_2} ) = f_{\rm u}( {\zeta_2} ) {\boldsymbol u}_2, \;$$

with a parabolic velocity shape function $f_{\rm u}( {\zeta_2} ) = s_0 + s_1\zeta _2 + s_2\zeta _2^2$ normalized as $\int_0^1 {f_{\rm u}( {\zeta_2} ) {\rm d}\zeta _2 = 1}$. The coefficients s ₀, s ₁, s ₂ may be chosen based on experimental observations (e.g. Hermann and Hutter, 1991).

The three shape factors f _ρu, f _ρuu, f _cζ in Eqns (10) and (11) account for non-uniform density and velocity profiles, and are respectively defined as

(14)$$f_{\rho {\rm u}} = {{\mathop \smallint \nolimits_0^1 \rho _{2\ast }( \zeta_2 ) {\boldsymbol u}_{2\ast }( \zeta_2 ) {\rm d}\zeta_2 } \over {\rho _2{\boldsymbol u}_2}} = \left({1-\displaystyle{{\rho_{\rm a}} \over {\rho_2}}} \right)\mathop \int \nolimits_0^1 f_{\rm c}( \zeta_2 ) f_{\rm u}( \zeta_2 ) {\rm d}\zeta_2 + \displaystyle{{\rho _{\rm a}} \over {\rho _2}}, \;$$

(15)$$\eqalign{\,f_{\rho {\rm uu}} = & \displaystyle{{\mathop \smallint \nolimits_0^1 \rho_{2\ast }( \zeta_2 ) {\boldsymbol u}_{2\ast }( \zeta_2 ) {\boldsymbol u}_{2\ast }( \zeta_2 ) {\rm d}\zeta_2 } \over {\rho _2{\boldsymbol u}_2{\boldsymbol u}_2}} \cr = & \left({1-\displaystyle{{\rho_{\rm a}} \over {\rho_2}}} \right)\mathop \int \limits_0^1 f_{\rm c}( \zeta_2 ) f_{\rm u}^2 ( \zeta_2 ) {\rm d}\zeta_2 + \displaystyle{{\rho _{\rm a}} \over {\rho _2}}\mathop \int \limits_0^1 f_{\rm u}^2 ( \zeta_2 ) {\rm d}\zeta_2 , \;} $$

(16)$$f_{{\rm c}\zeta } = \displaystyle{{\mathop \smallint \nolimits_0^1 ( \rho _{2\ast }( \zeta_2 ) -\rho _{\rm a}) \zeta_2 {\rm d}\zeta_2 } \over {\rho _2-\rho _{\rm a}}} = \mathop \int \nolimits_0^1 f_{\rm c}( \zeta_2 ) \zeta_2 {\rm d}\zeta_2 .$$

2.2. Closure assumptions

To solve the depth-averaged conservation equations, a total of 12 closure assumptions are required, expressing the stresses acting on the basal boundaries (${\boldsymbol \tau }_{01}, \;\;{\boldsymbol \tau }_{12\parallel }, \;\;{\boldsymbol \tau }_{02}$) and on top of the PSC (${\boldsymbol \tau }_{{\rm a}2\parallel }$, ${\boldsymbol \;}{\boldsymbol p}_{{\rm n}\parallel }$), and the volumetric exchange rates (E ₀₁, E ₀₂, D ₁₀, D ₂₀, S ₁₂, D ₂₁, E _a2). In our notation, the first subscript refers to the source layer exerting stress on, or feeding mass to, the target layer (second subscript). For the boundary stresses, the relationship $ {\boldsymbol \tau}_{ij}$ = $ {\boldsymbol \tau}_{ji}$ is generally not valid because of the jump conditions associated with entrainment.

2.2.1. Boundary stresses

The Voellmy model is used to model the shear stress acting at the base of the dense core:

(17)$${\boldsymbol \tau }_{01} =\left\{ \matrix{ \min ( {\tau}_{\rm c}, \;\;\mu \sigma_{zz} + k_{01}\rho_1{\Vert {\boldsymbol u_1} \Vert }^2 ) {\displaystyle{{\boldsymbol u}_1 \over \Vert {\boldsymbol u}_1 \Vert }} & {\rm if}\;h_0 > 0 \hfill \cr ( \mu \sigma_{zz} + k_{01}\;\rho_1 {\Vert {\boldsymbol u}_1 \Vert }^2 ) { \displaystyle{{\boldsymbol u}_1 \over \Vert {\boldsymbol u}_1 \Vert }} \hfill & {\rm if}\;h_0 = 0, } \right. $$

where μ is the Coulomb friction coefficient and k ₀₁ is the drag coefficient. The basal normal stress is defined as

(18)$$\sigma _{zz} = g_{z, 1}\rho _1h_1 + g_{z, 2}( {\rho_2-\rho_{\rm a}} ) h_2, \;$$

where the first term on the right-hand side represents the overburden weight of layer 1, while the second term on the right-hand side is the buoyant overburden weight of layer 2. As ρ ₁ ≫ ρ _a, the buoyancy effect is for simplicity neglected in the first term. Both terms are projected normal to the topography. The snow cover shear strength, $ { \tau}_{\rm c}$, is introduced in Eqn (17) for the case where the snow cover becomes entrained by the dense layer, that is, the upper expression applies only when h ₀ > 0 (more details will be given in Section 2.2.2).

The PSC is modelled as a turbulent fluid. A drag term is introduced to model the shear stress at the base of the PSC:

(19)$${\boldsymbol \tau }_{02} = \left\{ \matrix{\min ( {\tau}_{\rm c}, \;k_{02}\rho_2 \Vert {\boldsymbol u}_2 \Vert ^2 ) {\displaystyle{{\boldsymbol u}_2 \over \Vert {\boldsymbol u}_2 \Vert }} \hfill & {\rm if}\;h_0 > 0 \hfill \cr k_{02}\rho_2\Vert {\boldsymbol u}_2 \Vert {\boldsymbol u}_2 \hfill & {\rm if}\;h_0 = 0, \;} \right.$$

(20)$${\boldsymbol \tau }_{12\parallel } = \min ( { {\tau}_{{\rm su}}, \;\;k_{12}\rho_2{\Vert {{\boldsymbol u}_2-{\boldsymbol u}_1} \Vert }^2} ) \displaystyle{{{\boldsymbol u}_2-{\boldsymbol u}_1} \over {\Vert {{\boldsymbol u}_2-{\boldsymbol u}_1} \Vert }}, \;$$

where $ {\tau}_{\rm su}$ is the snow shear strength at the top of the dense layer, which is introduced in Eqn (20) to consistently model the suspension of snow from the dense core into the PSC (more details will be provided in Section 2.2.3, where the suspension model is described). Note that the shear stresses are evaluated at the layer boundary. However, by definition, the drag model makes use of the far-field velocities and densities (here assumed to correspond to the depth-averaged variables), and not of the variables evaluated at the boundaries.

The shear stress exerted on top of the PSC by the surrounding air is modelled as a drag term. Its projection in the flow direction is given by:

(21)$${\boldsymbol \tau }_{{\rm a}2\parallel } = k_{{\rm a}2}\rho _{\rm a}\left[ {1 + {( {\partial_xh} ) }^2 + {( {\partial_yh} ) }^2} \right]\Vert {{\boldsymbol u}_2} \Vert {\boldsymbol u}_2, \;$$

where h = h ₁ + h ₂. The term (1 + (∂_xh)² + (∂_yh)²) is introduced to provisionally model the speed-up of the air flow relative to the suspension-layer speed u₂ as the air gets deflected at the PSC–air interface (i.e. the deflected speed is enhanced by a factor (1 + (∂_xh)² + (∂_yh)²)^1/2). The speed-up of the air flow derives from conservation of the ambient air mass at the upper boundary of the PSC. The shear stress $ {\boldsymbol \tau}_{a2}$ is therefore projected along the flow direction of layer 2 (i.e. $ {\boldsymbol \tau}_{a2}$ is multiplied by a factor [1 + (∂_xh)² + (∂_yh)²]^−1/2) which acts on an enhanced upper surface compared to the basal area by a factor [1 + (∂_xh)² + (∂_yh)²]^1/2.

The stagnation pressure is defined as (Fig. 1)

(22)$${\boldsymbol p}_{\rm n} = {-} {1 \over 2} \rho _{\rm a} \left[ {{\boldsymbol u}_2\,{\cdot}\, \bigg( {{\boldsymbol \nabla }h \over \Vert {{\boldsymbol \nabla }h} \Vert } \bigg) } \right]^2 \bigg( {{\boldsymbol \nabla }h \over \Vert {{\boldsymbol \nabla }h} \Vert } \bigg) \left[ {1 + {( {\partial_xh} ) }^2 + {( {\partial_yh} ) }^2} \right]^{1/2}$$

where the term $-{\boldsymbol \nabla }h/\Vert {{\boldsymbol \nabla }h} \Vert$ defines the outward direction normal to the upper surface and $-{\boldsymbol u}_2\cdot ( {{\boldsymbol \nabla }h/\Vert {{\boldsymbol \nabla }h} \Vert } )$ is the speed of the PSC in the direction normal to the upper surface. The term [1 + (∂_xh)² + (∂_yh)²]^1/2 is again introduced to account for the enhanced upper surface compared to the basal area. p_n is therefore projected along the flow direction (u₂/‖u₂‖) to obtain

(23)$${\boldsymbol p}_{{\rm n}\parallel } = \displaystyle{1 \over 2}\rho _{\rm a}\left\{{-\displaystyle{{{\boldsymbol u}_2} \over {\Vert {{\boldsymbol u}_2} \Vert }}\cdot \displaystyle{{{\boldsymbol \nabla }h} \over {\Vert {{\boldsymbol \nabla }h} \Vert }}} \right\}^3 \left[ {1 + {( {\partial_xh} ) }^2 + {( {\partial_yh} ) }^2} \right]^{1/2} \Vert {{\boldsymbol u}_2} \Vert {\boldsymbol u}_2.$$

The Macaulay brackets are introduced to ensure that the stagnation pressure is only active in the flow regions impacting the air frontally and not on the avalanche tails.

2.2.2. Entrainment model

The tangential jump entrainment model – a slightly modified version of the formula derived by Fraccarollo and Capart (2002) – is used to compute the volumetric entrainment rate (per unit area) from the snow cover onto the dense core (i = 1) or onto the PSC (i = 2):

(24)$$E_{0i} = \displaystyle{{\{ {\Vert {{\boldsymbol \tau }_{i0}} \Vert -{\tau}_{\rm c}} \} } \over {\rho _0\Vert {{\boldsymbol u}_i} \Vert }}, \;$$

where ‖$ {\boldsymbol \tau}_{10}$‖ = μσ _zz + k ₀₁ρ ₁‖u₁‖² and ‖$ {\boldsymbol \tau}_{20}$‖ = k ₀₂ρ ₂‖u₂‖² are the shear stresses inside the flow just above the interface to the snow cover. If there is entrainment (i.e. when h ₀ > 0 and E _0i > 0), the shear stress in the snow cover just below the interface ($ {\boldsymbol \tau}_{01}$ or $ {\boldsymbol \tau}_{02}$) equals the snow cover shear strength $ {\tau}_{\rm c}$ (cf. Eqns (17) and (19)), as the momentum conservation equations are formulated for the mechanical system comprising the flow and the just-eroded bed layer. Entrainment hence requires ‖$ {\boldsymbol \tau}_{i0}$‖ > ‖$ {\boldsymbol \tau}_{0i}$‖ = $ { \tau}_{\rm c}$, the difference ‖$ {\boldsymbol \tau}_{i0}$‖ − $ {\tau}_{\rm c}$ serving to accelerate the eroded snow to the depth-averaged speed. In contrast to Fraccarollo and Capart (2002), here the shear strength is not assumed to be the Coulomb yield criterion but an intrinsic property of the perfectly brittle (Issler, 2014) snow cover. $ {\tau}_{\rm c}$ hence represents the cohesion – or, possibly, the undrained shear strength if pore air pressure produces a mechanical feedback on the stress state – in a Tresca-type yield model.

2.2.3. Suspension model

A tangential jump entrainment model is here proposed to compute the suspension rate from the dense core into the PSC. The PSC exerts the shear stress $\Vert {{\boldsymbol \tau }_{21\parallel }} \Vert = k_{12}\rho _2\Vert {{\boldsymbol u}_2-{\boldsymbol u}_1} \Vert ^2$ on the dense core. The dense core reacts with at most a strength τ _su (Fig. 2a), which represents the resistance of snow grains to suspension. If $\Vert {{\boldsymbol \tau }_{21\parallel }} \Vert$ exceeds $ { \tau}_{\rm su}$, snow grains within the dense core become eroded and suspended into the PSC, upon which their speed changes from u₁ to u₂ (Fig. 2b). The suspension rate is therefore given by:

(25)$$S_{12} = \displaystyle{{\{ {\Vert {{\boldsymbol \tau }_{21\parallel }} \Vert -\tau_{{\rm su}}} \} } \over {\rho _1\Vert {{\boldsymbol u}_2-{\boldsymbol u}_1} \Vert }}.$$

Figure 2. Schematic of suspension process: (a) vertical section of the powder snow cloud travelling over the dense core; (b) momentum change of the just-suspended layer; (c) velocity-dependent dense core disaggregation process and consequent onset of suspension.

The resistance of snow grains to suspension, $ { \tau}_{\rm su}$, is the only free parameter in the suspension model. We assume that $ { \tau}_{\rm su}$ evolves during the motion of the snow avalanche (Fig. 2c). At ‖u₁‖ = 0, the snow grains are sintered and therefore not easily suspended, that is, $ {\tau}_{\rm su}$ is large. At higher velocities, following the initial disaggregation of the slab due to internal shearing and collisions, larger snow clods and smaller interstitial snow grains form. At large enough velocities of the basal layer, the so-formed small snow grains can be transported upward by escaping air and be more easily suspended into the PSC by turbulent eddies. All these processes result in $ {\tau}_{\rm su}$ being small at large velocities of the dense core. When the dense core decelerates, new bonds may form between the snow clods of the dense layer through sintering, which hinders suspension, that is, $ {\tau}_{\rm su}$ becomes large again. We here assume the speed of the dense core as a proxy for the formation of small snow grains and their resistance to be suspended, through the following heuristic formulation:

(26)$$\tau _{{\rm su}} = \tau _{\rm s}{\rm e}^{-\gamma _{\rm s}\Vert {{\boldsymbol u}_1} \Vert }, \;$$

where γ _s is a decay coefficient (m⁻¹s). In the modelling, we typically set γ _s = 1 m⁻¹s, which, for plausible values of k ₁₂ ≈ 0.04 and internal slab shear strength of $ {\tau}_{\rm s}$ = 5 kPa, initiates suspension at u ₁ ≈ 8 m s⁻¹. An approximate threshold velocity of 10 m s⁻¹ for initiating the breakage of snow bonds and forming a PSC was inferred by Voellmy (1955) and Hopfinger (1983) from observations. The dense core speed is here assumed to only influence the disaggregation process. However, the increase of the flow speed, and more specifically of the shear rate, will also dilate or fluidize the dense core (Issler and Gauer, 2008), and hence decrease its density. This process is not yet included in our model but may be implemented in the future within the variable-density conservation equations of layer 1.

When ‖u₁‖ ≫ 1/γ _s, the top part of the dense layer becomes completely disaggregated ($ { \tau}_{\rm su}$ ≈ 0) and hence is easily suspended into the PSC. In this situation, Eqn (25) simplifies to:

(27)$$S_{12} = \displaystyle{{\rho _2} \over {\rho _1}}k_{12}\Vert {{\boldsymbol u}_2-{\boldsymbol u}_1} \Vert .$$

This equation has a similar structure as the suspension model used by Nazarov (1991), which reads:

(28)$$S_{12} = m_{12}\displaystyle{{\sqrt {\rho _1\rho _2} } \over {\rho _1 + \rho _2}}\Vert {{\boldsymbol u}_2-{\boldsymbol u}_1} \Vert , \;$$

where m ₁₂ is an empirical coefficient. In typical situations, ρ ₂ ≪ ρ ₁, and hence Eqn (28) can be approximated as $S_{12} = m_{12}\sqrt {\rho _2/\rho _1} \Vert {{\boldsymbol u}_2-{\boldsymbol u}_1} \Vert$. By equating this last expression to Eqn (27), one finds

(29)$$k_{12}\approx m_{12}\sqrt {\displaystyle{{\rho _1} \over {\rho _2}}} .$$

The ratio of the depth-averaged densities, ρ ₁/ρ ₂, is typically 20–150, and Nazarov (1991) reported values of m ₁₂ ≈ 0.01–0.10. Hence, one may expect k ₁₂ ≈ 0.04–1. While in Nazarov's model the suspension coefficient m ₁₂ and the drag coefficient k ₁₂ are set independently, a unique coefficient k ₁₂ controls both suspension and drag in the proposed suspension model. If one increases k ₁₂, the suspension rate will also increase. At the same time, however, the interfacial drag increases as well and reduces the velocity difference ‖u₂ − u₁‖, which will reduce the suspension rate.

2.2.4. Dense-core deposition model

To model deposition of the dense core, we assume a reverse entrainment model, inspired by the work of Nikooei and Choi (2022). The dense core exerts the basal shear stress ‖$ {\boldsymbol \tau}_{10}$‖ on the substrate, whose internal strength is $ { \tau}_{\rm du}$ (Fig. 3a). If $ { \tau}_{\rm du}$ > ‖$ {\boldsymbol \tau}_{10}$‖, an infinitesimal layer decelerates from the flow velocity ‖u₁‖ to rest (Fig. 3b), that is, it deposits. The deposition rate is therefore given by:

(30)$$D_{10} = \displaystyle{{\{ {\tau_{{\rm du}}-\Vert {{\boldsymbol \tau }_{10}} \Vert } \} } \over {\rho_{\rm d} \Vert {{\boldsymbol u}_1} \Vert }}.$$

Figure 3. Schematic of the deposition process (0 = snow cover, 1 = dense core): (a) vertical section of the dense core travelling over the snow cover; (b) momentum change of the deposited layer; (c) schematics of the velocity-dependent dense core sintering process and onset of deposition.

Similarly to Eqn (26), we assume that the internal strength of the deposited snow, $ { \tau}_{\rm du}$, is controlled by a sintering process (Fig. 3c), where the decrease of flow velocity allows bonding between snow clods and hence the recovery of shear strength:

(31)$$\tau _{{\rm du}} = \tau _{\rm d}{\rm e}^{-\gamma _{\rm d}\Vert {{\boldsymbol u}_1} \Vert }, \;$$

where $ { \tau}_{\rm d}$ is the final internal shear strength of the deposited snow and γ _d is a decay coefficient. A simplified version of Eqn (30) may be obtained from a first-order linearization of Eqn (31):

(32)$$D_{10} = \displaystyle{{\{ {\tau_{\rm d}( {1-\gamma_{\rm d}\Vert {{\boldsymbol u}_1} \Vert } ) -\Vert {{\boldsymbol \tau }_{10}} \Vert } \} } \over {\rho_{\rm d} \Vert {{\boldsymbol u}_1} \Vert }}.$$

Recently, Rauter and Köhler (2019) proposed an empirical deposition model, which is here rewritten for a mass block (i.e. neglecting the pressure gradient term):

(33)$$D_{10( {{\rm RK}} ) } = \left\{{1-\displaystyle{{\Vert {{\boldsymbol u}_1} \Vert } \over {u_{{\rm dep}}}}} \right\}\displaystyle{{\{ {\Vert {{\boldsymbol \tau }_{10}} \Vert -\rho_1\Vert {{\boldsymbol g}_\parallel } \Vert h_1} \} } \over {\rho _1\Vert {{\boldsymbol u}_1} \Vert }}.$$

We compare the two models considering a flow of frictional material (‖$ {\boldsymbol \tau}_{10}$‖ = μρ ₁g _zh ₁) on a horizontal plane (${\boldsymbol g}_\parallel{ = } \, {\bf 0}, \;\;g_z = g$). In this case, by equating Eqns (32) and (33), we get:

(34)$$\matrix{ {\tau _{\rm d} = 2\mu \rho _1gh_1, \;} \hfill \cr {\gamma _{\rm d} = \displaystyle{1 \over {2u_{{\rm dep}}}}, \;} \hfill \cr } $$

which defines the order of magnitude of the two parameters $ { \tau}_{\rm d}$ and γ _d. The influence of the deposition model on the block height and velocity is evaluated in Figure 4 for a block with initial height h ₁(0) = 1 m and initial velocity u ₁(0) = 5 m s⁻¹ decelerating on a horizontal plane.

Figure 4. Flow depth (continuous lines, scale on the left) and flow velocity (dashed lines, scale on the right) as a function of the runout distance for a frictional mass-block model using different deposition models and parameters.

The evolution of depth of each block model is obtained by time integration of the corresponding volumetric conservation equations, d_th ₁ = −D ₁₀. The deposition rate D ₁₀ is calculated using Eqns (30) and (31) for the tangential jump deposition model and Eqn (33) for the model of Rauter and Köhler (2019). The flow velocity is then calculated after integrating the momentum conservation equation in time. For the proposed tangential jump deposition model, the momentum conservation equation is:

(35)$${ {\rm d}_{\rm t} (h_1 u_1) = g_x h_1 - {{ \max\ ({\tau}_{\rm du}, \Vert {\boldsymbol \tau}_{10}\Vert )}\over{\rho_1}} = g_x h_1 - {{\Vert {\boldsymbol \tau}_{10}\Vert}\over{\rho_1}} - u_1 D_{10},}$$

where $ { \tau}_{\rm du}$ arises in the deposition situation (D ₁₀ > 0) to account for the jump in shear stress from $ {\boldsymbol \tau}_{10}$ to arrest the material. Instead, Rauter and Köhler (2019) and Nikooei and Choi (2022) do not implement such jump condition in their momentum conservation equation, d_t(h ₁u ₁) = g _xh ₁ − ‖$ {\boldsymbol \tau}_{10}$‖/ρ ₁, that is, they neglect the term −u ₁D ₁₀ on the right-hand side of the equation, which is however crucial in the deposition situation (Hungr, 1990; Erlichson, 1991).

The material properties of the block are as follows: ρ ₁ = 200 kg m⁻³, μ = 0.3. For Rauter and Köhler's (2019) model, u _dep = 3 m s⁻¹ is used, while Eqn (30) was first evaluated with $ { \tau}_{\rm d}$ = 1.2 kPa and γ _d = 1/6 m⁻¹ s (obtained from Eqn (34)), and then with a realistic value of the deposited snow shear strength, $ { \tau}_{\rm d}$ = 5.0 kPa, and using γ _d = 1 m⁻¹ s. The latter values seem to provide a reasonable deposition curve, with deposition starting when the dense core velocity drops below 2 m s⁻¹. In the proposed model, larger shear resistance acts at the base of the block in the deposition situation, which causes faster deposition and deceleration compared to Rauter and Köhler's model. Note that combining the momentum and volumetric conservation equations of the proposed model leads to the equation of motion d_tu ₁ = g _x − ‖$ {\boldsymbol \tau}_{10}$‖/(ρ ₁h ₁) = g _x − μg _z: the flow mobility only depends on the slope and on the friction parameter and is independent of deposition (although the block may arrest earlier if it fully deposits). In other terms, the mobility of a decelerating flow cannot increase if deposition occurs concomitantly. The main practical advantage of using the physics-based deposition model of Eqn (30) is that the deposition rate is negatively correlated with the flow velocity, resulting in an increase of deposition depth at low speeds.

2.2.5. Particle settling model

The volumetric particle settling rate (per unit area) is given by the product of the settling velocity w _s of the snow particles composing the PSC and the snow concentration:

(36)$$D_{2i} = w_{\rm s}\;c_{\rm b}\cos \theta \;, \;$$

where i = 0, 1. The term cos θ is introduced because particles settle along the direction of gravity. The concentration of snow is evaluated at the bottom of the PSC layer: c _b = c(ζ ₂ = 0) = c ₀(ρ ₂ − ρ _a)/(ρ _s − ρ _a). For the settling onto the dense core, Nazarov (1991) assumes ρ _s = ρ ₁. We also make this assumption, which implies that both the snow particles and part of their surrounding air are lost by the PSC. Hence, here, c _b is the concentration within the PSC of the snow particles and of part of their surrounding interstitial air (Fig. 5a). They settle onto the dense core at the constant density ρ ₁ (Fig. 5b). Since layers 1 and 2 are contiguous, and the ambient air is on top of them, such surrounding air must come from the PSC, which further implies that the volume lost by the PSC (D ₂₁ in Eqn (9)) should be equal to the volume gained by the dense layer (ρ _s/ρ ₁ D ₂₁ in Eqn (5)). Hence, ρ _s = ρ ₁ is necessary if the assumption ρ ₁ = const. is made (i.e. when using conservation Eqns (5) and (7)). Instead, relaxing the hypothesis ρ ₁ = const. (i.e. when using conservation Eqns (3), (4) and (6)), one may assume the density of snow particles to be in the range ρ ₁ ≤ ρ _s < ρ _ice, ρ _ice being the density of snow crystals. If ρ _s = ρ _ice, only the snow crystals settle (Fig. 5c). In our simulations, we stick to the first case (ρ ₁ = const.) and hence use ρ _s = ρ ₁. For simplicity and consistency, we extend the hypothesis ρ _s = ρ ₀ to evaluate the term D ₂₀ as well. In the case where the PSC is flowing directly on top of the snow cover, particles may either directly settle onto layer 0 and come to rest or, if the speed of layer 2 is large enough (in the code, we set a threshold speed of 10 m s⁻¹), settle to reform a moving layer 1.

Figure 5. Sketch of the mass and volume changes, in the time interval Δt, of layers 1 (dashed black line) and 2 (continuous black line) due to settling of snow particles at a speed w _s. The dotted red line is pure air volume which might be ejected from the PSC into the ambient air. The continuous red line is the solid volume effectively settling onto the dense core. (a) Volume components within the PSC. (b) Settling model with a constant density of the dense core, where snow particles and their interstitial air settle from the PSC. (c) Settling model with variable density of the dense core, where only snow crystals settle from the PSC. (d) Nazarov (1991) model with a constant density of the dense core, where snow particles and interstitial air settle from the PSC, and the excess air is ejected out of the PSC. The fate of such excess air is not clear yet.

Nazarov (1991) assumes ρ _s = ρ ₁ to evaluate the volume gained by the dense core due to settling in the unit of time and area:

(37)$$D_{21( 1 ) } = w_{\rm s}\;c_0\displaystyle{{\rho _2-\rho _{\rm a}} \over {\rho _1-\rho _{\rm a}}}\cos \theta .$$

This corresponds to Eqn (36) although in his original formulation, he uses c ₀ = 1. The corresponding mass gained by the dense core in his model is ρ ₁D ₂₁₍₁₎, consistent with our model. Instead, regarding the volume lost by layer 2 in the unit of time and area, Nazarov (1991) assumed

(38)$$D_{21( 2 ) } = w_{\rm s} \cos \theta .$$

The corresponding mass which is lost in his model is ρ ₂D ₂₁₍₂₎. Hence, in Nazarov (1991), the volumes exchanged between the two layers are not the same (Fig. 5d). In other terms, Nazarov assumed that, during settling, layer 1 increases its volume by an amount corresponding to the snow particles and of part of the surrounding air (ρ _s = ρ ₁), while layer 2 loses the volume associated to both the solid particles and their whole surrounding air (ρ _s = ρ ₂). Hence, part of the air lost by layer 2 in Nazarov (1991) settles onto layer 1, while some (excess) air is, implicitly in his model, ejected vertically into the ambient air. Such ejected air may form a density stratification within the PSC (e.g. Turnbull and others, 2007), which might influence the concentration profile function, with concentration vanishing at the top (i.e. f _c(ζ ₂ = 1) = 0). In our work, we assume that such excess surrounding air remains conserved within the PSC (layer 2) without any change of the profile functions.

2.2.6. Air entrainment model

Ellison and Turner (1959) modelled the entrainment of lighter fluid on top of a heavier turbulent layer assuming that the entrainment rate is proportional to the flow velocity and a function of the bulk Richardson number,

(39)$$E_{{\rm a}2} = f( {{\rm Ri}} ) \Vert {{\boldsymbol u}_2} \Vert , \;$$

where the bulk Richardson number is defined as

(40)$${\rm Ri} = \displaystyle{{( {\rho_2-\rho_{\rm a}} ) g h_2 \cos \theta} \over {\rho _{\rm a}{\Vert {{\boldsymbol u}_2} \Vert }^2}}.$$

Different parametrizations of the normalized ambient-fluid entrainment speed (the function f(Ri)) have been proposed in the literature (Ellison and Turner, 1959; Ancey, 2004; Dellino and others, 2019). Within our model, we use the empirical equation proposed by Parker and others (1987), which is a fit of the experimental data of Ellison and Turner (1959), Lofquist (1960) and Fukuoka and others (1980) over a wide range of Richardson numbers (10⁻² to 10²):

(41)$$f( {{\rm Ri}} ) = \displaystyle{{0.075} \over {{( {1 + 718\;{\rm R}{\rm i}^{2.4}} ) }^{0.5}}}.$$

Hence, air entrainment becomes significant for low values of the Richardson number. Using typical values of the powder-snow avalanche parameters, for example, $\rho_2 \approx 2\,{\rm kg}\,{\rm m}^{-3}$, h ₂ ≈ 20 m, ‖u₂‖ ≈ 50 m s^–1, g cosθ ≈ 7 m s^–2, the Richardson number is 0.04, and from Eqn (41) we calculate E _a2/‖u₂‖ = 0.07, which is within the range of observed growth rates of the PSC (Issler and others, 2020).

2.3. Numerical implementation

The partial differential Eqns (2), (5), (7) and (9–11) are solved using a simplified version of the Method of Transport, which is based on (Fey and Jeltsch, 1992). Volume, mass and momentum are advected to the neighbouring cells (including diagonal neighbours) with the flow layer speed u_i(t _n), similarly to an upwind scheme. In contrast to the original scheme, the flow variables are not decomposed into component waves propagating at relative speed $\sqrt {g_zh_i\;}$ with respect to the flow layer speed, but the pressure gradient is explicitly included in the momentum balance equations. This is an acceptable approximation at large Froude numbers typical of dry-snow avalanches. The time step Δt_n at time t_n is chosen according to the Courant–Friedrichs–Lewy condition; if negative flow depths occur, the time step is repeated with reduced Δt_n. Numerical instabilities are, however, observed in some simulations at the initiation of motion – in particular for large values of k ₁₂, which would lead to a large mass of snow being suspended from the dense layer – and during the final runout stages. The numerical instabilities were partly mitigated in the current work by using a small maximum time step (0.05 s), which however resulted in quite long computation times for the large-scale simulations on a 3D topography, O(10 min) on a single core. A version of the code accounting for wave propagation will be implemented in the future, which may reduce the numerical instabilities. The pressure gradient term is evaluated as a central difference at the cell faces, where the pressure at each face is evaluated as a geometric average of the pressure at neighbour cells. The source terms are evaluated at the time step t _n. The volume, mass and momentum in each cell are therefore calculated at the time step t _n+1 using a forward Eulerian integration. More details on the numerical implementation are provided in (Issler, in preparation).

3. Calibration of the numerical model

The two-layer model requires specifying 18 parameters, in addition to the snow cover properties, to simulate both the dense flow, the PSC and their interactions with each other and with the snow cover. To better constrain the influence of each model parameter on the flow mobility and PSC formation, we first tested the model on two simplified 2-D parabolic topographies, where the mobility, maximum speed (Eqn (A3)) and stagnation pressure (Eqn (A2)) of the PSC are computed (Appendix A). The possible ranges of values of the model parameters are also described in the appendix. Most of the model parameters influence the PSC dynamics to some degree, but a few of them have the strongest effect. The most mobile and destructive PSCs are obtained with small values of k ₀₁, which create a fast dense core, large values of k ₁₂ and small values of w _s, both of which increase the effective density of the PSC. Entrainment of the snow cover by the dense core also increases the mobility of the PSC and its impact pressure. As expected, the simulations showed that formation of a powerful PSC requires a substantial drop height. In summary, the simulations on the 2-D parabolic tracks showed large sensitivity of the simulation results on the model parameters, highlighting the need to carefully select the input parameters. The simplified tracks allow to efficiently investigate the PSC dynamics for a wide range of parameter combinations, but they do not allow to select the model parameters precisely for a given snow avalanche event, where 3D topographical effects are important. Hence, a mixed snow avalanche that occurred in Norway is back-calculated. A powder-snow avalanche was released spontaneously in the Knutstugrovi gully (Lom municipality, Norway) on 27 February 2020. A dashboard camera mounted in a car recorded the PSC moving across the ice-covered lake (Fig. 6a). The car was stopped 48 s afterwards by a traffic light connected to an early-warning system (Fig. 6b). The runout of the dense and fluidized components of the mixed snow avalanche extended ~200 m below the road (Fig. 6c), while the PSC travelled almost 1 km farther on the lake. The dense core left deposits of 1 m or more along the road, while an average deposit of 0.5 m – mostly associated with the fluidized layer – was measured below the road. The PSC deposits on the lake were only 1–3 cm thick. The gully was partly forested. Birches were broken by the pressure exerted by the dense core and PSC and transported downslope with the avalanche (Fig. 6d). Beyond the stopping point of the dense and fluidized layers, the PSC did not break any large trees and shrubs, but some small trees. The PSC also reached a few houses on the south-eastern side of its trajectory, without causing damage.

Figure 6. The powder-snow avalanche in the Knutstugrovi gully on 27 February 2020. (a) The powder snow cloud travelling on the lake, as recorded by a dashboard camera. (b) East-side limit of the powder snow cloud recorded by the dashboard camera 48 s later and 800 m farther than the photo in (a). (c) Extent of the deposits. (d) Trees broken and transported by the dense flow, and trees not damaged by the powder snow cloud. Images (a) and (b) courtesy J. M. Kveum, images (c) and (d) by Henrik Langeland, NGI.

The simulations were run on a digital terrain model with 5 m resolution, which was interpolated from an original resolution of 10 m. The release area and the fracture depth of 0.8 m were identified after a field survey, from which an approximate release volume of 30 000 m³ is obtained. Based on the available information from the field survey and the weather data, the erodible snow cover depth is set to vary linearly from 0.1 m at 350 m a.s.l. to 0.3 m at 1150 m a.s.l. Similarly, the snow cover shear strength is set to vary linearly from 1 kPa at 350 m a.s.l., corresponding to denser snow, to 0.5 kPa at 1150 m a.s.l., corresponding to a loose wind-drifted snow.

The model parameters used for the back-calculation of the avalanche were determined through trial-and-error simulations and are shown in Table 1, ‘20200227 back-calculated’. The simulation results are shown in Figure 7 and compared to field observations. Low friction parameters of the dense/fluidized layer are used to capture its large mobility ($\alpha _1 = 26.5^\circ$). The simulated extent of the deposit of the dense layer was observed to be mostly dependent on the friction coefficient μ, while the deposit shear strength $ { \tau}_{\rm d}$ controls the onset of deposition upslope, in agreement with the observations made for the depositing block model in Section 2.2.4. The ratio ρ _d/ρ ₁ influences the deposit thickness. Based on field observations, the deposit was observed to be quite compact, and hence $ { \tau}_{\rm d}$ = 2 kPa and ρ _d = 300 kg m⁻³ were used. The deposition depths from layer 1 compare quite well with the deposit thicknesses measured in the field (insert in Fig. 7b), with a root mean square deviation between the measured and modelled deposition depths of 0.35 m. In the simulation, the dense/fluidized layer comes to a halt at the base of the gully, leading to the detachment of the PSC, which then continues its runout on the lake.

Table 1. Material and model parameters used for the back-calculation of the Knutstugrovi powder-snow avalanche on 2020-02-27

Simulations are also run for a hypothetical wet avalanche and using other model assumptions (constant concentration and velocity profiles and Nazarov's (1991) settling model).

Figure 7. Back-calculated results of the 2020-02-27 Knutstugrovi powder-snow avalanche. (a) Computed dense core speed. (b) Computed deposition depths. The insert compares the simulation results with the measured deposition. (c) Computed PSC maximum speed. (d) Computed maximum pressure at the base of the PSC. Estimated pressures based on extent and damage are also indicated.

The PSC parameters were back-calculated to match the deposit thickness over the lake, the PSC extent and the forest damage. Loose wind-drifted snow characterized the release area, which should be easily suspended, and hence $ { \tau}_{\rm s}$ = 1 kPa and k ₁₂ = 0.05 were used. To model the PSC travelling to the opposite side of the lake, a mid-range value of the snow particle settling velocity was used, w _s = 0.15 m s⁻¹. Substantially larger values of w _s would cause the PSC to die too early, and much smaller values of w _s cause too high deposition depths from the PSC on the lake. The back-calculated PSC parameters allowed us to accurately model the deposit thickness beyond the dense core, 0.02 m, as observed in the field measurement (Fig. 7b). This implicitly validates the modelled amount of snow suspended from the dense core to the PSC.

The lateral extent of the PSC could be estimated from the video recording of the dashboard camera (Fig. 6b): it is indicated as ‘< 0.1 kPa’ in Figure 7d, corresponding to the estimated low PSC impact pressure, and is well captured by the model. The extent of the PSC on the south-eastern side of the path could be estimated based on the absence of damage at the cabins. McClung and Schaerer (2022) indicate that pressures above 0.5 kPa may break windows. The windows and the wooden structure did not sustain any damage, which hence allowed us to estimate a possible upper bound of the impact pressure at this location (point indicated as ‘< 0.5 kPa’), although larger pressures cannot be completely excluded.

Further constraints on the PSC pressure can be derived from damage to trees. The many trees broken along the gully do not provide a bound on the PSC pressure because they may have been broken by the dense core. Instead, trees beyond the extent of layer 1 can be used to directly validate the impact pressures of the PSC. The average pressure needed to break a tree may be calculated as (Feistl and others, 2015):

(42)$$p_{\rm t} = \displaystyle{{\pi D_{\rm t}^3 \sigma _{\rm t}} \over {16w_{\rm t}H_{\rm t}^2 }}, \;$$

where D _t is the diameter of the tree, H _t is its height (the hypothesis that the PSC is higher than the tree is here made), w _t is the effective crown width and σ _t is the tensile strength of the tree. Notice that to derive Eqn (42), the PSC pressure is considered to be uniform, ignoring the shape factors f _c and f _u: this assumption is valid if H _t ≪ h ₂. Smaller birch trees were observed to be broken by the PSC. For these, we roughly set D _t ≈ 0.2 m, H _t ≈ 8 m, w _t ≈ 2 m and σ _t ≈ 40 MPa to get a lower bound of the PSC maximum basal pressure, p _2b,max = 0.5 kPa. In contrast, bigger trees were not broken by the PSC (Fig. 6d). For these, we set approximate tree parameters D _t ≈ 0.4 m, H _t ≈ 12 m, w _t ≈ 3 m and σ _t ≈ 50 MPa to get an upper bound of the PSC maximum basal pressure, p _2b,max = 1.5 kPa. This provides an order of magnitude of the PSC maximum pressure (point indicated as ‘≈ 1 kPa’ in Fig. 7d), which seems to be well captured in the back-calculation. Furthermore, smaller shrubs were only bent by the PSC. For these, we use D _t ≈ 0.1 m, H _t ≈ 1 m, w _t ≈ 1 m and σ _t ≈ 30 MPa to get an upper bound of the PSC maximum pressure, p _2b,max = 6 kPa (point indicated as ‘< 6 kPa’ in Fig. 7d). The higher pressure needed to break shrubs is indicative of the ability of these tree species to grow and resist in areas with possible or frequent powder-snow avalanche activity, but it is less useful in constraining the back-calculation of the 2020-02-27 avalanche. In summary, the simulation results generated using the back-calculated parameters allow us to obtain a good match with the available field data from the 2020-02-27 Knutstugrovi avalanche.

Three additional simulations are carried out (Table 1) to explore the significance of some model parameters and assumptions. Figure 8 shows the extent of the PSCs, taken as the 0.5 kPa PSC isobar, for the four simulated avalanches.

Figure 8. 0.5 kPa PSC pressure contours for the four simulations listed in Table 1. The maximum basal pressure distribution for the simulation using Nazarov's (1991) settling model is also shown. The white dashed line is used in Figure 9 to plot longitudinal profiles of relevant parameters of the four simulations.

A wet snow avalanche is simulated using high friction parameters and higher values of $ { \tau}_{\rm s}$ and $ { \tau}_{\rm d}$, and higher settling speed of the PSC (corresponding to bigger suspended particles). A very weak PSC is generated, the extent of which is significantly smaller than the one of the 2020-02-27 powder-snow avalanche and limited to the initial steep section of the gully. Uniform concentration and velocity profiles produce smaller PSC pressures than non-uniform ones. Nazarov's (1991) and our settling models lead to similar runout areas and PSC pressures (compare the colour rendering of Figs 7d, 8).

Figure 9 shows the computed maximum densities and normalized maximum velocities of the PSC along the profile line of Figure 8 for the four simulations. The maximum density of the PSC (Fig. 9a) for the 20200227 back-calculated simulation is ~2.5 times the air density near the release area, where snow starts to get suspended. The maximum density then decreases along the gully to an almost stationary value of 1.7ρ _a due to much air being entrained. As the PSC reaches the flat runout area, the density progressively decreases to the air density, as settling dominates over suspension. The depth-averaged density of layer 2 is of a similar order of magnitude as that indicated by Sovilla and others (2015), which suggests that the avalanche formation processes are adequately captured by the model. For the model assumptions and parameters used in our simulations, the PSC does not entrain the snow cover (for $ { \tau}_{\rm c}$ = 0.5 kPa, k ₀₂ = 0.025 and ρ ₂ = 3 kg m⁻³, the PSC would start entraining only at speeds above 82 m s⁻¹). Instead, the entrainment model and parameters used in Nazarov (1991) and Eglit (1998) predict fairly large entrainment also by the PSC, which exceeds the entrainment rates by PSCs inferred by Issler and others (2020) by an order of magnitude. Hence, this may explain why quite high PSC densities are modelled by Nazarov (1991) compared to our numerical predictions (Fig. 9a). The normalized speed of the PSC (Eqn (A3), Fig. 9b) reaches a maximum value of 0.7 (i.e. a maximum speed of 43 m s⁻¹) ‒ similar to that measured for other major snow avalanches (Gauer, 2014) ‒ to then decrease to near-zero values at the end of the flat runout area. As expected, both the density and velocity of the PSC generated in the wet-snow avalanche simulation are lower compared to the simulation of the 2020-02-27 powder-snow avalanche.

Figure 9. Simulation results along the profile line indicated in Figure 8. (a) The computed maximum density of the PSC normalized by the air density (ρ _a = 1.225 kg m⁻³). (b) The computed maximum normalized velocity of the PSC (Eqn (A3)).

It is of interest to compare simulations with the back-calculated (1) and the Nazarov (1991) model (4), which generated similar pressures (Figs 7d, 8 respectively). The peak velocity profiles of simulations 1 and 4 are similar; however, their peak density profiles are quite different in magnitude. The pressure similarity can be explained by considering the asynchronous distributions of the computed flow variables of the PSC: while peak velocities are towards the flow front, peak densities typically lag behind. The body and tail are characterized by higher Richardson numbers and hence less air is here entrained: the peak density is affected by D _2i and is therefore different in the two models. Instead, at the flow front, in virtue of larger velocities, much air is entrained. The term E _a2 becomes dominant over the term D _2i, and hence the density at the flow front becomes almost independent of which settling model is used. As the maximum pressure is also observed to be at the flow front, the maximum pressures calculated using the two models are therefore similar. Note also that concentration and velocity profiles that both decrease with ζ ₂ (Fig. 11) imply f _ρu > 1 and therefore advection of denser flow from behind towards the flow front. Consequently, the density behind the front – and thus the maximum density in the PSC – decreases more rapidly than in simulation 3 with constant profiles (f _ρu = 1).

4. Discussion and conclusion

The principal area of application envisaged for MoT-PSA is hazard mapping, both as a supporting tool for experts assessing small areas in detail and as a key element in an automated chain of tools producing hazard indication maps over large areas. Thus, the following questions arise: (i) Is the modelling concept suitable for both intended application areas? (ii) Can MoT-PSA simulate mixed snow avalanches adequately? (iii) What needs to be improved before the model can be applied confidently by avalanche experts?

The first question can be answered in the affirmative: MoT-PSA can be used in the same way as MoT-Voellmy and similar modelling tools, except that the user must specify the values of additional parameters. While MoT-PSA is markedly slower than MoT-Voellmy due to more than double the number of differential equations, it is still fast enough to be used in the large-scale hazard mapping system NAKSIN (Issler and others, 2023). As it is implemented, MoT-PSA can also be used as a dense avalanche model by deactivating suspension (i.e. setting k ₁₂ = 0), which would reduce it to a similar mathematical model as MoT-Voellmy (Issler, in preparation), upon which MoT-PSA was built. Preliminary simulations furthermore suggest that the runout of the dense core with the PSC activated does not significantly differ from simulations conducted with the PSC deactivated. Compared to other operational avalanche models, for example, RAMMS::Avalanche (Bartelt and others, 2017), the dense avalanche model implemented in MoT-PSA includes entrainment and deposition. Finally, MoT-PSA is capable of modelling a detached PSC, initialized with height, density and speed distributions; in this mode, laboratory experiments on density currents (Lofquist, 1960; Beghin and Olagne, 1991) and suspension flows (Keller, 1995; Dellino and others, 2019) can be simulated to test the suspension-layer component of MoT-PSA in detail.

A definitive answer to the second question requires critically reviewing published measurements and reported observations of a wide variety of PSAs and then back-calculating them. This task remains to be done, but the good agreement of the simulation of the 2020 Knutstugrovi event with the observations suggests that the main features of PSAs relevant for hazard mapping are captured adequately with reasonable parameter values.

The model describes the complex flow of PSAs in terms of simple models for two distinct layers interacting with each other and with the snow cover. Ideally, each of these process models should be validated separately against dedicated experiments. In some cases, like the entrainment of air along the upper surface of the PSC, laboratory experiments covering the entire range of the relevant non-dimensional numbers have been carried out and are used in MoT-PSA. The approximations for the turbulent and pressure drag can be tested and possibly improved by comparing to 3D simulations of the Navier–Stokes equations with a suitable turbulence model. Direct snow entrainment from the snow cover into the PSC should be compatible with measurements of snow-cover erosion by blowing snow. Measurements with frequency-modulated radar provide data on the entrainment and deposition of snow into and from the dense/fluidized layer, but interpretation of the data is not straightforward and neither the shear stresses exerted by the avalanche nor the shear strength of the substrate are well constrained. For other processes – for example, the suspension of fine snow grains from the dense or fluidized layer – it remains to be investigated whether dedicated experiments are feasible. It is generally accepted that the Voellmy friction law, which has been used provisionally in MoT-PSA, describes the dissipative processes in real avalanches poorly and that the density in the dense/fluidized layer varies strongly in space and time. These deficiencies in turn affect the process models for entrainment and deposition of the dense layer, and for suspension of snow into the PSC. Therefore, only partial validation at the process level is possible at present and the adequacy of the model must be assessed mainly by back-calculating the main features of observed events.

A crucial step for making MoT-PSA usable in practice is to develop recommendations for choosing the large number of model parameters in a given situation. In particular, the snow-cover shear strength, $ {\tau}_{\rm c}$, and erodible snow depth, h ₀, will depend on the target return period of the simulated event and on the local climate. The same will apply for the mean settling velocity of snow grains and for the shear strength $ { \tau}_{\rm su}$ at the upper surface of the dense/fluidized layer. At present, a reliable calibration of a Voellmy-type model with entrainment and deposition is also lacking. One may, however, reasonably expect ongoing work to soon provide reasonable procedures for determining the input data and selecting the model parameters.