1. Introduction
Mixed snow avalanches (henceforth also called powder snow avalanches) are characterized by three distinct flow regimes (Sovilla and others, Reference Sovilla, McElwaine and Louge2015): 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, Reference Bartelt2017) 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, Reference Kulikovskiy and Sveshnikova1977). Eglit (Reference Eglit and Shahinpoor1983) and Nazarov (Reference Nazarov1991) 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 (Reference Fukushima and Parker1990) combined the centre-of-mass approach of Kulikovskiy and Sveshnikova with the four-equation model of turbidity currents by Parker and others (Reference Parker, Fukushima and Pantin1986) to model non-Boussinesq PSCs. In addition to the volume, mass and momentum conservation equations already identified by Eglit (Reference Eglit and Shahinpoor1983), Fukushima and Parker (Reference Fukushima and Parker1990) 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 (Reference Naaim and Gurer1998) 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 (Reference Sampl and Zwinger2004) 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 (Reference Bartelt, Buser, Vera Valero and Bühler2016) proposed a two-layer depth-averaged model for powder snow avalanches, whose underlying physical assumptions are, however, debated (Issler and others, Reference Issler, Jenkins and McElwaine2018).
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, Reference Sovilla, McElwaine and Louge2015) and deposits (Issler and others, Reference Issler, Gauer, Schaer and Keller2020). 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, Reference Issler1998) 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 (Reference Eglit and Shahinpoor1983) and Nazarov's (Reference Nazarov1991) two-layer model to 3-D terrain (Issler, Reference Issler2023). 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, Reference Schaerer and Salway1980; Issler and others, Reference Issler, Gauer, Schaer and Keller1996, Reference Issler, Gauer, Schaer and Keller2020; Sovilla and others, Reference Sovilla, McElwaine and Louge2015) – 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.
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
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−3. Snow clods in the dense/fluidized layer have considerably smaller density, ρ s, in the range 200–600 kg m−3 (McClung and Schaerer, Reference McClung and Schaerer1985) but may be embedded in a matrix of snow grains (Issler and others, Reference Issler, Gauer, Schaer and Keller2020). 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, ρ 0 to be constant, and the dense/fluidized layer, ρ 1, 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 ρ 0 is eroded to a depth Δh 0 by the dense-flow layer with density ρ 1 and flow depth h 1, ρ 0 remains the same but ρ 1 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 (ρ 1 ≈ ρ 0), the density of layer 1 does not change after erosion ($ \rho_1^{\prime}$
≈ ρ 1), and mass and volume are hence exactly conserved. However, with typical values for the dense core, ρ 0 ≈ 0.5ρ 1, Δh 0 ≈ 0.3h 1, one obtains $ \rho_1^{\prime}$
≈ 0.85ρ 1; similarly, for a dilute, fluidized flow with ρ 0 ≈ 3ρ 1, $ \rho_1^{\prime}$
≈ 1.5ρ 1 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}$
= ρ 1), 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 ρ 0h 0(x, y, t), ρ 1h 1(x, y, t), ρ 2h 2(x, y, t), momentum hold-ups ρ 1h 1u1(x, y, t) = (ρ 1h 1u 1, ρ 1h 1v 1) and ρ 2h 2u2(x, y, t) = (ρ 2h 2u 2, ρ 2h 2v 2), and bulk density of the PSC, ρ 2(x, y, t). The density of the snow cover, ρ 0, and of the dense core, ρ 1, are instead assumed to be constant.
In the presence of an erodible bed, the evolution of the bed depth (h 0) can be written as an Exner equation:
where E 01 and E 02 represent the entrainment rates (m s–1) of the bed into layer 1 and layer 2, respectively, D 10 represents the deposition rate of the dense layer, and D 20 represents the settling rate of snow particles from the PSC onto the snow cover. If ρ 0 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
where S 12 is the suspension rate from the dense layer into the PSC and D 21 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 ρ 1. Hence, Eqn (3) is not considered. Instead, using the condition ρ 1 = const., the mass conservation Eqn (4) simplifies to
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 ρ 1 = const. When the snow cover is entrained, it loses a volume E 01 per unit time and area, while the dense core gains a volume ρ 0/ρ 1 E 01. For instance, in the case ρ 1 < ρ 0, 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 10 per unit time and area, while the flow loses a volume ρ d/ρ 1 D 10 per unit time and area. The situation where ρ d > ρ 1 physically corresponds to the dense core compressing and sintering to form the final deposit, as is observed in nature (Issler and others, Reference Issler, Gauer, Schaer and Keller2020).
The momentum conservation for the dense layer is written as
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, Reference Gray, Wieland and Hutter1999), but was for simplicity not implemented in the model at this stage.
Introducing the hypothesis ρ 1 = const., Eqn (6) reduces to
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:
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
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 pn 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 pn|| 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 ζ 1 = (z − h 0)/h 1, that is, ρ 1*(ζ 1) = ρ 1 and u1*(ζ 1) = u1 (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 ζ 2 = (z − h 0 − h 1)/h 2. The density profile of the PSC is assumed as
where c = (ρ 2 − ρ 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 0, c 1, c 2 may be chosen based on laboratory experiments (e.g. Hermann and Hutter, Reference Hermann and Hutter1991) 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
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 0, s 1, s 2 may be chosen based on experimental observations (e.g. Hermann and Hutter, Reference Hermann and Hutter1991).
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
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 01, E 02, D 10, D 20, S 12, D 21, 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:
where μ is the Coulomb friction coefficient and k 01 is the drag coefficient. The basal normal stress is defined as
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 ρ 1 ≫ ρ 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 > 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:
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:
where h = h 1 + h 2. The term (1 + (∂xh)2 + (∂yh)2) is introduced to provisionally model the speed-up of the air flow relative to the suspension-layer speed u2 as the air gets deflected at the PSC–air interface (i.e. the deflected speed is enhanced by a factor (1 + (∂xh)2 + (∂yh)2)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)2 + (∂yh)2]−1/2) which acts on an enhanced upper surface compared to the basal area by a factor [1 + (∂xh)2 + (∂yh)2]1/2.
The stagnation pressure is defined as (Fig. 1)
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)2 + (∂yh)2]1/2 is again introduced to account for the enhanced upper surface compared to the basal area. pn is therefore projected along the flow direction (u2/‖u2‖) to obtain
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 (Reference Fraccarollo and Capart2002) – 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):
where ‖$ {\boldsymbol \tau}_{10}$
‖ = μσ zz + k 01ρ 1‖u1‖2 and ‖$ {\boldsymbol \tau}_{20}$
‖ = k 02ρ 2‖u2‖2 are the shear stresses inside the flow just above the interface to the snow cover. If there is entrainment (i.e. when h 0 > 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 (Reference Fraccarollo and Capart2002), here the shear strength is not assumed to be the Coulomb yield criterion but an intrinsic property of the perfectly brittle (Issler, Reference Issler2014) 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 u1 to u2 (Fig. 2b). The suspension rate is therefore given by:
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 ‖u1‖ = 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:
where γ s is a decay coefficient (m−1s). In the modelling, we typically set γ s = 1 m−1s, which, for plausible values of k 12 ≈ 0.04 and internal slab shear strength of $ {\tau}_{\rm s}$
= 5 kPa, initiates suspension at u 1 ≈ 8 m s−1. An approximate threshold velocity of 10 m s−1 for initiating the breakage of snow bonds and forming a PSC was inferred by Voellmy (Reference Voellmy1955) and Hopfinger (Reference Hopfinger1983) 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, Reference Issler and Gauer2008), 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 ‖u1‖ ≫ 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:
This equation has a similar structure as the suspension model used by Nazarov (Reference Nazarov1991), which reads:
where m 12 is an empirical coefficient. In typical situations, ρ 2 ≪ ρ 1, 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
The ratio of the depth-averaged densities, ρ 1/ρ 2, is typically 20–150, and Nazarov (Reference Nazarov1991) reported values of m 12 ≈ 0.01–0.10. Hence, one may expect k 12 ≈ 0.04–1. While in Nazarov's model the suspension coefficient m 12 and the drag coefficient k 12 are set independently, a unique coefficient k 12 controls both suspension and drag in the proposed suspension model. If one increases k 12, the suspension rate will also increase. At the same time, however, the interfacial drag increases as well and reduces the velocity difference ‖u2 − u1‖, 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 (Reference Nikooei and Choi2022). 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 ‖u1‖ to rest (Fig. 3b), that is, it deposits. The deposition rate is therefore given by:
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:
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):
Recently, Rauter and Köhler (Reference Rauter and Köhler2019) proposed an empirical deposition model, which is here rewritten for a mass block (i.e. neglecting the pressure gradient term):
We compare the two models considering a flow of frictional material (‖$ {\boldsymbol \tau}_{10}$
‖ = μρ 1g zh 1) on a horizontal plane (${\boldsymbol g}_\parallel{ = } \, {\bf 0}, \;\;g_z = g$
). In this case, by equating Eqns (32) and (33), we get:
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 1(0) = 1 m and initial velocity u 1(0) = 5 m s−1 decelerating on a horizontal plane.
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, dth 1 = −D 10. The deposition rate D 10 is calculated using Eqns (30) and (31) for the tangential jump deposition model and Eqn (33) for the model of Rauter and Köhler (Reference Rauter and Köhler2019). 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:
where $ { \tau}_{\rm du}$
arises in the deposition situation (D 10 > 0) to account for the jump in shear stress from $ {\boldsymbol \tau}_{10}$
to arrest the material. Instead, Rauter and Köhler (Reference Rauter and Köhler2019) and Nikooei and Choi (Reference Nikooei and Choi2022) do not implement such jump condition in their momentum conservation equation, dt(h 1u 1) = g xh 1 − ‖$ {\boldsymbol \tau}_{10}$
‖/ρ 1, that is, they neglect the term −u 1D 10 on the right-hand side of the equation, which is however crucial in the deposition situation (Hungr, Reference Hungr1990; Erlichson, Reference Erlichson1991).
The material properties of the block are as follows: ρ 1 = 200 kg m−3, μ = 0.3. For Rauter and Köhler's (Reference Rauter and Köhler2019) model, u dep = 3 m s−1 is used, while Eqn (30) was first evaluated with $ { \tau}_{\rm d}$
= 1.2 kPa and γ d = 1/6 m−1 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−1 s. The latter values seem to provide a reasonable deposition curve, with deposition starting when the dense core velocity drops below 2 m s−1. 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 dtu 1 = g x − ‖$ {\boldsymbol \tau}_{10}$
‖/(ρ 1h 1) = 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:
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(ζ 2 = 0) = c 0(ρ 2 − ρ a)/(ρ s − ρ a). For the settling onto the dense core, Nazarov (Reference Nazarov1991) assumes ρ s = ρ 1. 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 ρ 1 (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 21 in Eqn (9)) should be equal to the volume gained by the dense layer (ρ s/ρ 1 D 21 in Eqn (5)). Hence, ρ s = ρ 1 is necessary if the assumption ρ 1 = const. is made (i.e. when using conservation Eqns (5) and (7)). Instead, relaxing the hypothesis ρ 1 = const. (i.e. when using conservation Eqns (3), (4) and (6)), one may assume the density of snow particles to be in the range ρ 1 ≤ ρ 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 (ρ 1 = const.) and hence use ρ s = ρ 1. For simplicity and consistency, we extend the hypothesis ρ s = ρ 0 to evaluate the term D 20 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−1), settle to reform a moving layer 1.
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 (Reference Nazarov1991) 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 (Reference Nazarov1991) assumes ρ s = ρ 1 to evaluate the volume gained by the dense core due to settling in the unit of time and area:
This corresponds to Eqn (36) although in his original formulation, he uses c 0 = 1. The corresponding mass gained by the dense core in his model is ρ 1D 21(1), consistent with our model. Instead, regarding the volume lost by layer 2 in the unit of time and area, Nazarov (Reference Nazarov1991) assumed
The corresponding mass which is lost in his model is ρ 2D 21(2). Hence, in Nazarov (Reference Nazarov1991), 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 = ρ 1), while layer 2 loses the volume associated to both the solid particles and their whole surrounding air (ρ s = ρ 2). Hence, part of the air lost by layer 2 in Nazarov (Reference Nazarov1991) 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, Reference Turnbull, McElwaine and Ancey2007), which might influence the concentration profile function, with concentration vanishing at the top (i.e. f c(ζ 2 = 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 (Reference Ellison and Turner1959) 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,
where the bulk Richardson number is defined as
Different parametrizations of the normalized ambient-fluid entrainment speed (the function f(Ri)) have been proposed in the literature (Ellison and Turner, Reference Ellison and Turner1959; Ancey, Reference Ancey2004; Dellino and others, Reference Dellino, Dioguardi, Doronzo and Mele2019). Within our model, we use the empirical equation proposed by Parker and others (Reference Parker, Garcia, Fukushima and Yu1987), which is a fit of the experimental data of Ellison and Turner (Reference Ellison and Turner1959), Lofquist (Reference Lofquist1960) and Fukuoka and others (Reference Fukuoka, Fukushima, Murata and Arai1980) over a wide range of Richardson numbers (10−2 to 102):
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 2 ≈ 20 m, ‖u2‖ ≈ 50 m s–1, g cosθ ≈ 7 m s–2, the Richardson number is 0.04, and from Eqn (41) we calculate E a2/‖u2‖ = 0.07, which is within the range of observed growth rates of the PSC (Issler and others, Reference Issler, Gauer, Schaer and Keller2020).
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, Reference Fey and Jeltsch1992). Volume, mass and momentum are advected to the neighbouring cells (including diagonal neighbours) with the flow layer speed ui(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 Δtn at time tn is chosen according to the Courant–Friedrichs–Lewy condition; if negative flow depths occur, the time step is repeated with reduced Δtn. Numerical instabilities are, however, observed in some simulations at the initiation of motion – in particular for large values of k 12, 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 01, which create a fast dense core, large values of k 12 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.
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 m3 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/ρ 1 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−3 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.
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 (Reference Nazarov1991) settling model).
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 12 = 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−1. 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 (Reference McClung and Schaerer2022) 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, Reference Feistl2015):
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 2. 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.
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 (Reference Nazarov1991) 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 (Reference Nazarov1991) 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 (Reference Sovilla, McElwaine and Louge2015), 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 02 = 0.025 and ρ 2 = 3 kg m−3, the PSC would start entraining only at speeds above 82 m s−1). Instead, the entrainment model and parameters used in Nazarov (Reference Nazarov1991) and Eglit (Reference Eglit1998) predict fairly large entrainment also by the PSC, which exceeds the entrainment rates by PSCs inferred by Issler and others (Reference Issler, Gauer, Schaer and Keller2020) by an order of magnitude. Hence, this may explain why quite high PSC densities are modelled by Nazarov (Reference Nazarov1991) 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−1) ‒ similar to that measured for other major snow avalanches (Gauer, Reference Gauer2014) ‒ 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.
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−3). (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 (Reference Nazarov1991) 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 ζ 2 (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, Reference Issler2023). As it is implemented, MoT-PSA can also be used as a dense avalanche model by deactivating suspension (i.e. setting k 12 = 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, Reference Bartelt2017), 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, Reference Lofquist1960; Beghin and Olagne, Reference Beghin and Olagne1991) and suspension flows (Keller, Reference Keller1995; Dellino and others, Reference Dellino, Dioguardi, Doronzo and Mele2019) 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 0, 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.
Acknowledgements
This work was supported by NGI's special grant for snow avalanche research from the Norwegian Petroleum and Energy Department, administrated by the Norwegian Directorate of Waterways and Energy (NVE), and NGI's general research fund from the Research Council of Norway. H. V. thanks Peter Gauer for discussions on the scaling behaviour of snow avalanches, and Betty Sovilla for discussions on mixed snow avalanches. D. I. is grateful to Margarita E. Eglit, Peter Gauer, Hansueli Gubler, Felix Hermann, Kolumban Hutter, Mark Schaer and Betty Sovilla for illuminating and enjoyable discussions on powder-snow avalanches over the course of many years. We are grateful to two anonymous reviewers for their careful reading of the manuscript and insightful comments, which helped to improve this paper.
Appendix A – Sensitivity analysis on a 2-D parabolic track
A sensitivity study is conducted on simplified parametrized 2-D topographies to evaluate the influence of the model parameters on the dynamics of the powder snow cloud. A typical avalanche path may be approximated as a parabola (Lied and Bakkehøi, Reference Lied and Bakkehøi1980). Gauer (Reference Gauer2018) proposed the following parametrization for a path with a constant-slope release area, parabolic track and horizontal runout (Fig. 10):
where x is the horizontal coordinate (with value 0 at the front of the release area) and z is the vertical coordinate (with the horizontal runout plane set at altitude 0). θ 0 is the inclination of the release area, and H sc,f is the total drop height from the front of the release area. Two different parabolic topographies are tested: P1 characterized by H sc,f = 1500 m and $\theta _0 = 50^\circ$
, which is a long and steep avalanche track and should favour the formation of large powder snow clouds; P2 characterized by H sc,f = 700 m and $\theta _0 = 45^\circ$
, which is a shorter avalanche track, where smaller powder snow clouds are expected. To keep the test configuration as simple as possible, 2-D simulations are run along a cross-section of the parabola. This setup is only representative of an avalanche that does not show significant transverse spreading.
Parabolic track and release area.
The initial avalanche volume plays an important role in the flow mobility. For the calibration, we simulate ‘major’ avalanche events characterized by large volumes. In most of our simulations, we deactivate the snow cover entrainment, to focus on the essential features of the powder snow cloud formation. For all the simulations, a (normal-to-slope) fracture depth (D r) of 1 m is used. Large slab avalanches are typically characterized by a width (cross-slope) to length (down-slope) ratio between 2 and 6 (McClung and Schaerer, Reference McClung and Schaerer2022). For our simulation, we select a ratio of 2. Hence, the longitudinal length of the release area is determined as $L = \sqrt {V_{\rm r}/( {2D_{\rm r}} ) }$
, which is distributed along a plane inclined at θ 0, on top of the parabolic track. Since 2-D simulations are performed, the actual volumes per unit width used in the simulations are equal to V r/(2L). A grid with a uniform horizontally projected cell length of 5 m is used.
A sensitivity study is carried out, computing quantities characterizing the flow dynamics for different values of the model parameters. The ranges of the input model parameters and their default values (in bold) are shown in Table 2. The concentration profile is heuristically defined based on experiments by Hermann and Hutter (Reference Hermann and Hutter1991): it varies linearly from 4/3 at the base of the PSC to 2/3 at the top of the PSC layer (Fig. 11). A nose-shaped velocity profile is assumed for the PSC, which has a value of 1.4 at the base and 0.13 at the top. The density of the dense core is constant and assumed in the range of 50–300 kg m−3, the lower value being used to model a fluidized flow. Typical literature values of the friction parameters of the dense core, μ and k 01, are assumed (Bartelt and others, Reference Bartelt2017). The drag coefficient between the dense core and the powder snow cloud, k 12, is assumed to be ~10–20 times the value of k 01 (depending on the simulation). This accounts for the fact that the powder snow cloud is governed by turbulence, while the dense core is typically dominated by Coulomb friction. The drag coefficient between the powder snow cloud and the snow cover, k 02, is fixed to 0.5 ⋅ k 12 to account for a smoother snow cover surface compared to the rougher dense core surface. The magnitude of k 02 in our study is similar to that reported by Fukushima and Parker (Reference Fukushima and Parker1990). In preliminary simulations, we furthermore observed that k 02 does not significantly affect the dynamics of the powder-snow avalanche, which is also in agreement with the results by Fukushima and Parker (Reference Fukushima and Parker1990), and hence k 02 was not included in the sensitivity study. Similarly, k a2 has negligible influence on the results and was not varied in the study. The settling velocity of snow particles is assumed between 0.05 and 0.50 m s−1. Note that our model does not explicitly consider turbulence keeping particles in suspension. To compensate for this, lower values of the settling speed probably must be used. In one of the simulations, entrainment of a 0.5 m thick, weak snow cover ($ { \tau}_{\rm c}$
= 700 Pa) is activated. In all the simulations, deposition of the dense core is kept active (with ρ d = ρ 1, $ { \tau}_{\rm d}$
= 5 kPa and γ d = 1 m−1 s). Deposition was anyway observed to have negligible effects on the simulation results. The internal shear strength of the dense core at rest was assumed equal to $ { \tau}_{\rm s}$
= 5 kPa (with γ s = 1 m−1 s). One additional simulation is run with $ { \tau}_{\rm s}$
= 0 to study the influence of snow shear strength on the suspension mechanism.
Summary of material and model parameters used in the sensitivity study on the 2-D parabolic tracks
Concentration and velocity profiles assumed for the simulations.
Different flow dynamics parameters are calculated in each simulation, as proxies for the flow mobility and PSC dynamics:
• The flow mobility is measured by the runout angle α (tan α = H/l, where H is the drop height from the top of the release area to the front of the avalanche, and l is the horizontal distance between the two points). An α angle is calculated for the PSC (α 2) which represents its mobility. α 2 is calculated from the distal limit where the stagnation pressure at the base of the PSC,
(A2)$$p_{2{\rm b}} = \displaystyle{1 \over 2}\rho _{2{\rm b}}u_{2{\rm b}}^2 \;$$drops below 0.5 kPa. The maximum basal stagnation pressure (p 2b,max) is also computed at the end of the parabola ($0^\circ$
).
• The normalized maximum flow velocity is calculated for the PSC as (Gauer, Reference Gauer2018, Reference Gauer2020)
(A3)$$\upsilon _{2, {\rm max}} = \displaystyle{{u_{2, {\rm max}\;}} \over {\sqrt {gH_{{\rm sc}}/2} }}.$$
The results of the sensitivity study are shown in Figures 12 and 13 for the long and steep track P1. As expected, the flow mobility of the PSC (Fig. 12) is in most cases larger than that of the dense core ($\alpha _1 \approx 37^\circ$
in the default simulation), and is sensitive to all the parameters of the model. High values of the Voellmy-friction parameters reduce the velocity of the dense layer, and they consequently also decrease the mobility of the PSC. This dependency is particularly significant for the turbulence coefficient k 01, which has the largest influence on the dense layer velocity and is hence observed to be negatively correlated to the maximum velocity of the PSC. k 12 is instead positively correlated to the mobility and maximum velocity of the PSC. High values of k 12 favour suspension of snow from the dense core: hence, the effective density of the PSC becomes larger, and so does the driving force due to gravity whereas the braking effect of air entrainment depends on ρ a but not on ρ 2. A high value of the dense core internal shear strength ($ { \tau}_{\rm s}$
= 5 kPa) retards the formation of the PSC and increases the shear resistance at the base of the powder snow cloud (cf. Eqn (20)), thereby reducing its mobility. Large values of the snow particle settling velocity (e.g. w s = 0.5 m s−1) reduce the effective density of the PSC and its acceleration along the steep part of the parabola and thus its mobility. Entrainment of snow (mostly by the dense core) nourishes the PSC and hence increases its mobility and velocity. An ignition effect due to snow cover entrainment is also reported by Fukushima and Parker (Reference Fukushima and Parker1990).
Results of the sensitivity analysis on the runout angle of the powder snow cloud, α 2, on the avalanche track P1.
Results of the sensitivity analysis of the maximum basal pressure of the powder snow cloud at the end of the parabola, $p_{2{\rm b, max}}( {0^\circ } )$
, on the avalanche track P1.
Figure 12 shows that the mobility of the PSC varies strongly with the dense core density ρ 1 at values below 100 kg m−3. In particular, ρ 1 = 50 kg m−3 represents a fluidized layer. The low value of ρ 1 causes the density of the PSC to be low, which in turn causes large air entrainment, reducing the effective driving force. The large entrainment of air for ρ 1 = 50 kg m−3 also produces a retarding effect on the PSC, whose mobility is hence reduced. However, this effect should be considered an artefact of the simplified description of the dense/fluidized layer: fluidization occurs in the head of the avalanche, accompanied by significantly larger flow depths than in the dense core and often with high entrainment rates. Thus, even though the fine snow particles are all the more easily suspended if the head is completely fluidized, the mass in the fluidized head will nevertheless not decrease, as was the case in the simulation with MoT-PSA. This shortcoming of the model will need to be addressed when the dense/fluidized layer is modelled more accurately with spatially and temporally variable density.
Figure 13 shows the maximum basal stagnation pressure (p 2b,max) computed at the end of the parabola ($0^\circ$
). As discussed above, a weak PSC (in terms of density and velocity) is formed for the simulations with a ‘slow’ dense core (k 01 = 0.005), for low suspension capability (k 12 = 0.02), for heavy snow particles characterized by high settling speed (w s = 0.5 m s−1) and for the simulation of the fluidized layer (ρ 1 = 50 kg m−3). For these simulations, the maximum basal stagnation pressure is lower than 0.5 kPa. Instead, when a mature PSC forms, it accelerates more significantly along the steep parts of the topography, and the stagnation pressures are consequently higher (up to 2 kPa at the end of the parabola in the simulations). In agreement with observations, destructive PSCs form from fast-moving avalanches of loose and dry snow (i.e. with high suspension and low settling).
The simulations for the shorter avalanche track P2 are presented in Figures 14 and 15. The track profile was chosen to be similar to the cross-section of the Knutstugrovi avalanche, which is back-calculated in Section 3. The trends for all the parameters are like the simulations on P1. However, on P2 it is observed that the simulated powder snow clouds are smaller and less mobile compared to the powder snow clouds in P1. This is explained by lower velocities reached on the shorter avalanche track P2, which lead to less snow being suspended from the dense core, causing a lower density of the powder snow cloud and hence lower driving force. On tracks with smaller drop heights, a mature powder snow cloud can only form from dry snow avalanches, that is, using low values of μ, k 01, w s, and high values of k 12 in the model. Like the simulations on P1, snow cover entrainment produces a more mobile PSC. On the shorter track, however, the snow cover needs to be significantly weaker than on larger avalanche tracks for entrainment to happen.
Results of the sensitivity analysis on the runout angle of the powder snow cloud, α 2, on the avalanche track P2.
Results of the sensitivity analysis on the maximum basal pressure of the powder snow cloud at the end of the parabola, $p_{2{\rm b, max}}( {0^\circ } )$
, on the avalanche track P2.
Open access
