**Research article**
12 Oct 2020

**Research article** | 12 Oct 2020

# The mechanical origin of snow avalanche dynamics and flow regime transitions

Xingyue Li Betty Sovilla Chenfanfu Jiang and Johan Gaume

^{1},

^{2},

^{3},

^{1,2}

**Xingyue Li et al.**Xingyue Li Betty Sovilla Chenfanfu Jiang and Johan Gaume

^{1},

^{2},

^{3},

^{1,2}

^{1}School of Architecture, Civil and Environmental Engineering, Swiss Federal Institute of Technology, Lausanne, Switzerland^{2}WSL Institute for Snow and Avalanche Research, SLF, Davos, Switzerland^{3}Computer and Information Science Department, University of Pennsylvania, Philadelphia, Pennsylvania, USA

^{1}School of Architecture, Civil and Environmental Engineering, Swiss Federal Institute of Technology, Lausanne, Switzerland^{2}WSL Institute for Snow and Avalanche Research, SLF, Davos, Switzerland^{3}Computer and Information Science Department, University of Pennsylvania, Philadelphia, Pennsylvania, USA

**Correspondence**: Johan Gaume (johan.gaume@epfl.ch)

**Correspondence**: Johan Gaume (johan.gaume@epfl.ch)

Received: 24 Mar 2020 – Discussion started: 28 Apr 2020 – Revised: 04 Aug 2020 – Accepted: 25 Aug 2020 – Published: 12 Oct 2020

Snow avalanches cause fatalities and economic damage. Key to their mitigation is the understanding of snow avalanche dynamics. This study investigates the dynamic behavior of snow avalanches, using the material point method (MPM) and an elastoplastic constitutive law for porous cohesive materials. By virtue of the hybrid Eulerian–Lagrangian nature of the MPM, we can handle processes involving large deformations, collisions and fractures. Meanwhile, the elastoplastic model enables us to capture the mixed-mode failure of snow, including tensile, shear and compressive failure. Using the proposed numerical approach, distinct behaviors of snow avalanches, from fluid-like to solid-like, are examined with varied snow mechanical properties. In particular, four flow regimes reported from real observations are identified, namely, cold dense, warm shear, warm plug and sliding slab regimes. Moreover, notable surges and roll waves are observed peculiarly for flows in transition from cold dense to warm shear regimes. Each of the flow regimes shows unique flow characteristics in terms of the evolution of the avalanche front, the free-surface shape, and the vertical velocity profile. We further explore the influence of slope geometry on the behavior of snow avalanches, including the effect of slope angle and path length on the maximum flow velocity, the runout angle and the deposit height. Unified trends are obtained between the normalized maximum flow velocity and the scaled runout angle as well as the scaled deposit height, reflecting analogous rules with different geometry conditions of the slope. It is found that the maximum flow velocity is mainly controlled by the friction between the bed and the flow, the geometry of the slope, and the snow properties. We reveal the crucial effect of both flow and deposition behaviors on the runout angle. Furthermore, our MPM modeling is calibrated and tested with simulations of real snow avalanches. The evolution of the avalanche front position and velocity from the MPM modeling shows reasonable agreement with the measurement data from the literature. The MPM approach serves as a novel and promising tool to offer systematic and quantitative analysis for mitigation of gravitational hazards like snow avalanches.

Snow avalanches have long been threatening infrastructures and human lives. Buildings, roads and railways can be severely damaged, causing profound economic losses. Moreover, the fatalities induced by snow avalanches are significant and are about 100 annually in the European Alps during the last 4 decades (Techel et al., 2016). Due to climate change, the frequency and risk of snow avalanches are still increasing (Choubin et al., 2019). It is thus of great importance to mitigate snow avalanche hazards, which highly relies on the understanding of their complex dynamic behavior.

Snow can behave as a fluid or as a solid under different conditions, leading to distinct behaviors of snow avalanches in reality (Gaume et al., 2011; Ancey, 2016). Characterizing different flow regimes of snow avalanches has played a significant role in hazard mapping and the design of mitigation measures (Gauer et al., 2008). Traditionally, two flow regimes were considered for snow avalanches, namely, dense snow avalanches and powder snow avalanches. A recent study by Köhler et al. (2018) highlighted the role of snow temperature in classifying the flow regimes and extended the traditional classification. The starting, flowing and stopping signatures of snow avalanches were used to distinguish seven flow regimes, including four dense flow regimes, two powder flow regimes and a snowball flow regime. Although the flow regimes have been identified based on macro flow behavior from the real measurements, their underlying physics remains unclear. Numerical and theoretical models can provide efficient and comprehensive analysis to shed light on the internal flow characteristics underpinning the macro flow behavior but are extremely challenging (Faug et al., 2018). To date, there has been no recognized model capable of capturing and analyzing the diverse behaviors of snow avalanches in a systematic and well-controlled way. Furthermore, the crucial effects of snow mechanical property and terrain geometry on snow avalanche dynamics have been widely recognized but are sparsely investigated due to practical challenges. Only limited numerical and real-measurement studies have been reported (Keshari et al., 2010; Fischer et al., 2012, 2015; Steinkogler et al., 2014).

Popular classical numerical tools for modeling snow avalanches primarily apply two-dimensional (2D) depth-averaged methods based on shallow-water theory (Naaim et al., 2013; Rauter et al., 2018), which fail to capture important flow characteristics along the surface-normal direction such as velocity distribution (Eglit et al., 2020). Nevertheless, 2D models are computationally efficient and provide acceptable accuracy, which serve as a powerful tool in many applications like hazard mapping. In comparison, three-dimensional (3D) simulations can fully resolve flow variations in all dimensions, which consequently require longer computation time. In recent years, particle-based continuum methods, including the smoothed particle hydrodynamics (SPH) method, the particle finite element method (PFEM) and the material point method (MPM), have gained increasing popularity in avalanche modeling, as they are able to easily handle large deformations and discontinuities (Abdelrazek et al., 2014; Salazar et al., 2016; Gaume et al., 2018). In particular, the MPM has proven to be an effective and efficient tool in investigating snow (Stomakhin et al., 2013; Gaume et al., 2018, 2019). Compared with the SPH method where boundary conditions are challenging to generalize, the MPM can readily address complex boundaries (Raymond et al., 2018). Moreover, the MPM does not suffer from the time-consuming neighbor searching that is inevitable in many mesh-free approaches like the SPH method (Mast et al., 2014). Both the PFEM and MPM use a set of Lagrangian particles and a background mesh to solve the mass and momentum conservation of a system. In contrast to the PFEM, each particle in the MPM has fixed mass, which naturally guarantees mass conservation. However, the fixed mass meanwhile leads to difficulty in adding or removing particles from the system (Larsson et al., 2020). The computational cost of the MPM is lower than that of the PFEM according to simulations with the same formulation (Papakrivopoulos, 2018).

This study applies the MPM in 2D (slope parallel and slope normal) to explore the distinct behaviors of snow avalanches and the key controlling factors of snow avalanche dynamics. To facilitate efficient computation and capture important flow features along the surface-normal direction, our 2D MPM modeling neglects variations along the flow width. The MPM is a hybrid Eulerian–Lagrangian approach, which uses Lagrangian particles to track mass, momentum and deformation gradient, and adopts an Eulerian background mesh to solve and update the motion of the particles. By virtue of the hybrid Eulerian–Lagrangian nature of the MPM, processes with large deformations, fractures, collisions and impacts can be well simulated (Mast et al., 2014; Gaume et al., 2018, 2019). In addition, continuous solid–fluid phase transition and the coexistence of solid-like and fluid-like behaviors can be captured with the implementation of proper constitutive models (Stomakhin et al., 2013; Gaume et al., 2018). The MPM has been increasingly adopted to investigate gravity-driven flows like landslides, debris flows and avalanches (Soga et al., 2016; Abe and Konagai, 2016; Gaume et al., 2018). This study will highlight the capability of the MPM in capturing different flow regimes of snow avalanches from fluid-like shear flow to solid-like sliding slab, by adopting a finite-strain elastoplastic constitutive model. Furthermore, it will be demonstrated that the proposed numerical approach serves as a promising tool to systematically study the key influencing factors of snow avalanche dynamics, including snow mechanical properties and slope geometry.

## 2.1 The material point method (MPM)

Assuming a continuous material, the MPM discretizes it into Lagrangian particles (material points) to trace mass, momentum and deformation gradient and adopts Eulerian grids to solve the motion of the particles and update their states. In particular, the particle motion is governed by mass and momentum conservation as follows:

where *ρ* is density, *t* is time, ** v** is velocity,

**is the Cauchy stress and**

*σ***is the gravitational acceleration. As the mass carried by each particle does not vary, the balance of mass is satisfied naturally. The momentum balance is solved with a regular background Eulerian grid mesh and the discretization of the weak form of Eq. (2). The explicit MPM algorithm by Stomakhin et al. (2013) is applied with a symplectic Euler time integrator. Details of the adopted MPM time-stepping algorithm can be found in Stomakhin et al. (2013), Jiang et al. (2016) and Gaume et al. (2018). Note compared to Gaume et al. (2018), this study uses the affine particle-in-cell (APIC) method (Jiang et al., 2015), by which angular momentum is preserved in addition to linear momentum.**

*g*The MPM relies on a continuum description and requires a constitutive
model for the considered material. The Cauchy stress ** σ** in
Eq. (2) is related to the strain through
an elastoplastic constitutive law as follows:

where Ψ is the elastoplastic potential energy density,
**F**_{E} is the elastic part of the deformation gradient
**F** and *J*=det(**F**). Note that various constitutive
models can be implemented into the framework of the MPM to capture
different materials and their distinct behaviors. For example, a
nonassociated Mohr–Coulomb model was applied to model landslide and
dam failure (Zabala and Alonso, 2011; Soga et al., 2016) and a
nonassociated Drucker–Prager model was used to simulate sand
(Klár et al., 2016). In this study, we use the associated modified
cam clay model developed for snow (Gaume et al., 2018), which
reproduces mixed-mode snow fracture and compaction hardening.

## 2.2 Finite-strain elastoplastic model

The elastoplastic model in this study is borrowed from
Gaume et al. (2018), which consists of a mixed-mode
shear–compression yield surface, a hardening law and an associative
flow rule. We recall the main characteristics of the three key
components. On the basis of laboratory experiments
(Reiweger et al., 2015) and simulations based on X-ray computed
tomography (Hagenmuller et al., 2015; Chandel et al., 2015; Srivastava et al., 2017), the yield surface is defined in the space of the
*p*–*q* invariants of the stress tensor as follows:

*p* is the pressure calculated as $p=-\text{tr}\left(\mathit{\tau}\right)/d$,
where ** τ** is the Kirchhoff stress tensor and

*d*is the dimension.

*q*is the von Mises stress defined as $q=(\mathrm{3}/\mathrm{2}\phantom{\rule{0.25em}{0ex}}\mathbf{s}:\mathbf{s}{)}^{\mathrm{1}/\mathrm{2}}$, where $\mathbf{s}=\mathit{\tau}+p\mathbf{I}$ is the deviatoric stress tensor and

**I**is the identity matrix.

*p*

_{0}is the consolidation pressure and denotes the isotropic compressive strength.

*β*

*p*

_{0}is the isotropic tensile strength, where

*β*reflects the cohesion.

*M*is the slope of the critical state line, which characterizes the internal friction.

When the *p*–*q* state of the material is inside the yield surface
(i.e., $y(p,q)<\mathrm{0}$), the material behaves elastically and follows
Hooke's law (St Venant–Kirchhoff Hencky strain). Plastic behavior
happens if $y(p,q)=\mathrm{0}$. Depending on the volumetric plastic strain
${\mathit{\u03f5}}_{\mathrm{v}}^{\mathrm{p}}$, hardening or softening is implemented by expanding or
shrinking the yield surface according to the following hardening law:

where *K* is the bulk modulus and *ξ* is the hardening factor. Under
compression (${\dot{\mathit{\u03f5}}}_{\mathrm{v}}^{\mathrm{p}}<\mathrm{0}$), *p*_{0} increases, leading to
hardening and promoting compaction. Under tension (${\dot{\mathit{\u03f5}}}_{\mathrm{v}}^{\mathrm{p}}>\mathrm{0}$), *p*_{0} decreases, resulting in softening and allowing fracture.

A flow rule needs to be adopted when plastic behavior occurs. Referring to Gaume et al. (2018), this study uses an associative plastic flow rule reported by Simo (1992) and Simo and Meschke (1993). The applied flow rule follows the principle of maximum plastic dissipation, which maximizes the rate of plastic dissipation. It is worth noting that the second law of thermodynamics is fully satisfied by using the plastic model with the flow rule. More details can be found in Gaume et al. (2018).

## 3.1 Model setup

To examine the behavior of snow avalanches under well-controlled
conditions, the setup with a rectangular snow sample and an ideal slope
is adopted as illustrated in Fig. 1. The snow
sample is initially placed at the top of the slope and will flow down
under gravity. The inclined slope is connected to the horizontal
ground using a circular arc with a central angle equaling the slope
angle *θ*. The arc and the horizontal ground are named the connecting-arc zone and deposition zone, respectively, in the following discussion. To investigate different flow regimes of snow avalanches, the
properties of snow are systematically varied, including the friction
coefficient *M*, the tension ∕ compression ratio *β*, the
hardening factor *ξ* and the initial consolidation pressure
${p}_{\mathrm{0}}^{\text{ini}}$. In addition, the effect of slope angle *θ*
and horizontal length *L*_{0} in Fig. 1 is
studied with five groups of simulations as summarized in
Table 1. For each group, the snow properties are
changed within the prescribed ranges. Groups I, II and III are conducted
to study the influence of slope angle *θ*, while groups II, IV and
V are designed to examine the effect of the horizontal length
*L*_{0}. When *θ* is varied in groups I, II and III, the horizontal
length *L*_{0} is fixed and the drop height *H*_{0} is adjusted as listed
in Table 1. Alternatively, one could fix the drop
height *H*_{0} and change the horizontal length *L*_{0}. Note *L*_{0}∕*h*_{0},
*L*_{0}∕*l*_{0} and *L*_{0}∕*r* are kept constant when *L*_{0} is varied in
groups II, IV and V by changing *h*_{0}, *l*_{0} and *r* accordingly. An
increased *L*_{0} leads to the scale-up of the setup, resulting in the
rise in the drop height *H*_{0}. Detailed parameters adopted in the MPM
simulations are summarized in Table 1. The size
of the background Eulerian mesh in the MPM is selected to be small enough
to guarantee the simulation accuracy and resolution and meanwhile be
large enough to shorten the computation time. The time step is
constrained by the CFL condition and the elastic wave speed to secure
the simulation stability. The simulation data are exported every
1∕24 s.

## 3.2 Typical flow regimes

In each group of our MPM simulations, four typical flow regimes are captured with the changing mechanical properties of snow. Figure 2a shows four representative cases in Group II, where distinct flow regimes are observed with the different snow properties summarized in Table 2. From top to bottom, we observe Regime 1 to Regime 4. The flow in Regime 1 behaves as a fluid or a dry cohesionless granular flow, whose free surface is continuous. Since the height of the flow in Regime 1 is excessively small compared with the others, it is scaled up to be 3 times higher along the bed-normal direction for better visualization in Fig. 2a. Small surges are observed especially at the front of the flow in Regime 1. The flow in Regime 2 demonstrates a more fluctuated free surface and a discontinuous tail, due to the occurrence of a granulation process. The flow height of the granular flow is higher compared with that of the flow in Regime 1, since the granules can be notably accumulated in the connecting-arc and deposition zones. The flow in Regime 3 demonstrates ductile behavior and slides down the slope and reaches the horizontal deposition zone with no significant deformation and no cracks. In contrast, clear cracks and broken pieces are noticed in the flow in Regime 4. The initial snow sample in Regime 4 breaks into multiple blocks shortly after its release from the slope.

The four flow regimes in Fig. 2a show
similar features to the identified flow regimes by
Köhler et al. (2018) based on real measurements and
observations, namely, cold dense, warm shear, warm plug and sliding
slab regimes. The detailed macro flow characteristics and internal shearing behavior of the flows will be expatiated in the following
section. The information of the energy of the flows is provided in
Appendix A. In addition, flows in transition between the
flow regimes can be captured. For example,
Fig. 2b shows a flow in transition
from cold dense to warm shear flow regimes, as it demonstrates the
characteristics of both regimes. Significant surges and roll waves
occur in Fig. 2b, showing similarity
with the cold dense flow. On the other hand, small granules are
observed at the early stage of the flow (see Supplement Video S4),
presenting characteristics of the warm shear flow. This transitional
flow is modeled with *M**β*=0.15 (*M*=0.5; *β*=0.3), whose
properties are in between the cold dense and warm shear flows
listed in Table 2. Although
only one flow regime is characterized for each of the flows in
Fig. 2a, it is worth noting that a
single snow avalanche may display different flow regimes at its
different locations and at different instants
(Kern et al., 2009). Figure 2c shows
the warm shear flow in Fig. 2a at its
early stage (*t*=5.5 s), where the red and blue materials are in compression and tension, respectively, and the dark gray denotes the
initial state of the material. The red particles are hardly visible as
they are located at the bottom layer of the flow, which indicates that
the snow inside the core of the avalanche is mainly under tension or
at the initial state. At *t*=5.5 s, the flow demonstrates the
characteristics of sliding slab, as the dark-gray part slides along
the slope and is seldom sheared. The initial sliding slab in
Fig. 2c can indeed transform into the
warm shear flow in Fig. 2a after it
reaches the deposition zone.

All the demonstrated flows in the four typical regimes share identical
initial and boundary conditions except for the snow properties. From
the simulations, it is clear that higher *M* and *β* give a more
frictional and cohesive flow, since they reflect the internal friction
and cohesion, respectively. For instance, the *M* and *β* of the
flow in the warm shear regime are higher than those of the cold dense
regime, which facilitates the formation of the granules and the higher
flow height. The hardening coefficient *ξ* and the initial
consolidation pressure ${p}_{\mathrm{0}}^{\text{ini}}$ also affect the flow
behavior, whose influence depends on the *M* and *β* according to
our sensitivity study. As listed in
Table 2, the tensile strength
$\mathit{\beta}{p}_{\mathrm{0}}^{\text{ini}}$ and *M**β* consistently increase from the
cold dense to the sliding slab flow regimes, which gives indications of
the possible underpinning factors controlling the transition of the
flow regimes.

### 3.2.1 Front evolution

Figure 3 illustrates the evolution of the
front position and velocity for the flows in the four typical flow
regimes in Fig. 2a. In some of the
simulated flows, scattered particles are observed at the flow front
(i.e., the warm shear flow and sliding slab flow in Supplement
Video S2), which need to be excluded in the determination of the front
position as they separate from the main body of the flow. Hence, the
front position is determined by ruling out 1 % of the
particles at the front of the flow. The front in
Fig. 3a is calculated as the horizontal
distance between the current front position and the initial front
position. The gray band in Fig. 3a shows the
region where a flow front enters the connecting-arc zone, below and
above which the flow front is on the slope and in the horizontal
deposition area, respectively. The evolution of front velocity in
Fig. 3b is plotted with two constant
velocities ${v}_{\text{max}}^{\mathrm{b}}$ and ${v}_{\text{max}}^{\mathrm{f}}$, which are the
theoretical velocities of a sliding rigid block with and without, respectively,
consideration of bed friction. Referring to
Fig. 1, if a rigid block slides down the slope,
its path length prior to the connecting-arc zone is $l={L}_{\mathrm{0}}/\text{cos}\mathit{\theta}-r\text{tan}\mathit{\theta}-\mathrm{0.5}{l}_{\mathrm{0}}$. Its acceleration
along the flow direction is ${a}^{\mathrm{b}}=g(\text{sin}\mathit{\theta}-\mathit{\mu}\text{cos}\mathit{\theta})$ considering gravity and
friction or *a*^{f}=*g*sin*θ* with a frictionless bed, where
*μ* is the bed friction coefficient fixed to 0.5 as listed in
Table 1. Given *l*, *a*^{b} and *a*^{f}, the theoretical
velocities when the block goes to the end of the slope can be
calculated as ${v}_{\text{max}}^{\mathrm{b}}=\sqrt{\mathrm{2}{a}^{\mathrm{b}}l}$ and
${v}_{\text{max}}^{\mathrm{f}}=\sqrt{\mathrm{2}{a}^{\mathrm{f}}l}$ with and without, respectively, consideration of the bed friction.

When the flows are on the slope (front<53.82 m), the front of the cold dense flow is the fastest, followed by the warm shear and warm plug flows, and the sliding slab is the slowest. Indeed, the front velocity in Fig. 3b generally gives a consistent trend. The warm shear and sliding slab flows demonstrate fluctuations at around 1.2 and 2.2 s in Fig. 3b, due to the breakage of the snow sample. The tensile strength of the snow in the warm shear flow (9 kPa) is smaller than in the sliding slab flow (21 kPa), leading to the earlier breakage and earlier fluctuations in Fig. 3b. The fluctuations in the front velocity of the cold dense flow last longer, which might be due to the turbulence and surges as observed in Fig. 2a. After the flows enter the connecting-arc zone ($\mathrm{53.82}\phantom{\rule{0.125em}{0ex}}\mathrm{m}<\text{front}<\mathrm{60.25}$ m), the fronts of the cold dense, warm plug and sliding slab flows evolve smoothly in Fig. 3a, while the front of the warm shear flow is sharply slowed down at the end of the connecting-arc zone before it goes to the horizontal deposition zone. In the warm shear flow, discrete granules form a discontinuous flow front. After the scattered granules arrive at the end of the connecting-arc zone, they quickly stop, leading to the stagnancy in the increase in the front position until the arrival of continuously incoming granules. Indeed, the warm shear flow in Fig. 3b shows the sharp velocity reduction at around 7.5 s, after which the growth in velocity occurs thanks to the incoming flow. A significant drop is also observed from the front velocity of the cold dense flow at around 4.8 s in Fig. 3b, corresponding well to the moment that the flow front reaches the connecting-arc zone in Fig. 3a. Indeed, this reduction is mainly contributed to by the changed slope geometry in addition to the turbulence and surges in the cold dense flow. In Fig. 3b, the maximum front velocities of the warm shear, warm plug and sliding slab are close to ${v}_{\text{max}}^{\mathrm{b}}$, while the maximum front velocity of the cold dense flow reaches ${v}_{\text{max}}^{\mathrm{f}}$. This indicates different dominant factors of the maximum flow velocity. The maximum velocity of the flows in warm shear, warm plug and sliding slab regimes is chiefly governed by the frictional behavior between the flow and the bed. In contrast, the maximum velocity of the cold dense flow is mainly controlled by the snow properties, where the low friction and low cohesion facilitate a higher velocity. When the flow fronts enter the deposition zone (front>60.25 m), all the flows start to slow down gradually. It is noticed that the front velocity of the warm shear flow shows fluctuations from around 7.5 s, which are chiefly because of the discrete nature of the snow granules at the flow front (see Supplement Video S1). As the front velocity of the warm shear flow decreases at around 7.5 s, the warm plug flow exceeds the warm shear flow. Nevertheless, the final front of the warm shear flow goes further as it stops later. The final fronts of the four flows show a consistent relation like that obtained when they are on the slope. Before the flows stop, the decelerations of the fronts (slope of velocity in Fig. 3b) are similar, which might be governed by bed friction.

### 3.2.2 Evolution of free-surface shape and vertical velocity profile

Figure 4 shows the evolution of the free surface
and velocity profile of the flows in the different flow regimes. The
velocity profile is obtained at *x*=50 m. The free surface
of the cold dense flow in Fig. 4a is scaled up
15 times along the bed-normal direction to visualize the fluctuations
at the surface. The height of the cold dense flow is much smaller than
the initial flow height, since it is highly sheared throughout its
flow depth as shown in Fig. 4e. The velocity
profiles in Fig. 4e are smooth, indicating the
continuous shearing along the flow depth. The shape of the velocity
profile in Fig. 4e does not change much with
time, while the flow speed and the shear rate decrease as the flow
tends to stop. Generally speaking, the cold dense flow behaves as a
fluid or a noncohesive granular flow, in agreement with the
characterization by Köhler et al. (2018). The warm shear flow in
Fig. 4b demonstrates a fluctuated free surface
because of the granules formed. Correspondingly, its velocity profile
shows fluctuations as well. As illustrated in
Fig. 4f, the warm shear flow is fully sheared
along the flow depth direction before it stops. Moreover, its flow
depth can exceed the initial flow height due to the piling up and
accumulation of snow granules. The shear behavior and the piling up
feature are indeed consistent with the warm shear regime identified by
Köhler et al. (2018). Nevertheless, instead of a noncohesive
granular flow characterized by Köhler et al. (2018), the flow in our
MPM modeling has cohesion (see
Table 2), which helps the
formation of the granules. The warm plug flow remains a block and is
seldom sheared when it slides on the slope. Upon its arrival at the
connecting-arc zone, significant shearing occurs due to the changed
shape of the connecting-arc zone. As shown in
Fig. 4c and g, the front of the warm plug flow
is notably sheared at *t*=8.0 s, the flow body is only
sheared at the bottom layer at *t*=8.3 s and the flow tail
is seldom sheared at *t*=8.8 s. The sliding slab in
Fig. 4d shows the sliding down of the slabs
from *t*=9.0 to 9.2 s and the accumulated slabs in the
connecting-arc and deposition zones at *t*=10.5 s. As there
are particles stopping on the slope, the tail of the free surface
collapses onto the slope. The shearing behavior inside the slabs is
extremely limited as shown in Fig. 4h. Both the
warm plug and the sliding slab behave as solid-like objects, while the
snow of the sliding slab flow is more brittle and produces slab
fractures.

## 3.3 Effect of slope angle and path length on flow dynamics

### 3.3.1 Maximum velocity and deposition height

Figures 5 and
6 demonstrate the evolution of the
maximum avalanche velocity on the slope *v*_{max} with the
normalized avalanche deposit height *h*_{d}∕*h*_{0}, under the effects of
the slope angle *θ* and the horizontal length *L*_{0} (reflecting
the path length), respectively. The deposit height *h*_{d} is defined as
the maximum avalanche height along the bed-normal direction after snow
avalanches stop. The maximum velocity of snow avalanches is usually
obtained before their arrival at the deposition zone. To analyze our
MPM data with theoretical predictions, the maximum velocity
*v*_{max} in Figs. 5 and
6 is determined when the flow is on
the slope. For all the simulated groups in
Figs. 5 and
6, a similar trend is observed, with
three zones indicating different flow characteristics and flow
regimes. Zone A shows a maximum velocity which tends to be negatively
correlated with the deposit height. A typical flow regime is the cold
dense regime, in which a higher maximum flow velocity leads to a
longer runout distance and a smaller deposit height. The black data
in Zone A reflect small snow friction and cohesion, which agrees with
the snow properties of the cold dense flow. Zone B has a deposit
height close to the initial height of the snow sample. This
characteristic is normally captured from the warm plug and sliding
slab flow regimes. Note that, in addition to the typical case in the
sliding slab regime in Fig. 4d which
demonstrates accumulated snow in the deposition zone, there are other
cases with slabs sliding down the slope and stopping in the deposition
zone without piling up and snow accumulation. These cases in the
sliding slab regime give a final deposit height close to the initial
flow height. The high snow friction and cohesion reflected by the
light color in Zone B indeed indicate the snow properties of the warm
plug and sliding slab flows. In Zone C, the deposit height is notably
larger than the initial height. In this case, representative flow
regimes are the warm shear flow and the sliding slab flow, where the
accumulation of snow can be significant after the flows deposit. It is
found that the flowing and deposition behaviors of snow avalanches are
primarily controlled by the snow friction and snow cohesion (*M* and
*β*), as we observe the clear transition of colors denoting
*M**β* in the different zones in
Figs. 5 and
6. The scattered colors of some
points, such as the dark points in Zone C, indicate the additional
effects of snow brittleness (reflected by *ξ*) and snow compressive
strength (*p*_{0}).

Slope angle is a key factor in evaluating the trigger, flow and
deposition of snow avalanches
(Gaume, 2012; Sovilla et al., 2010). Figure 5
shows the positive correlation between the slope angle *θ* and
the maximum velocity on the slope *v*_{max}. When *θ* is
varied with a fixed *L*_{0} (see Fig. 1), the
drop height *H*_{0} is increased accordingly, which gives a larger
initial potential energy of the flow and consequently a higher
*v*_{max}. The effect of increased path length reflected by
*L*_{0} is similar to the outcome of the growth of *θ*, as shown in
Fig. 6. It is interesting to observe
the similar trend for the different groups of simulations with varying
*θ* and *L*_{0}, which hints at an analogous physical rule behind the
trend. Indeed, a unified relation can be obtained as shown in
Fig. 7, by scaling *v*_{max}
and *h*_{d} as follows:

The normalization of *v*_{max} takes into account
${v}_{\text{max}}^{\mathrm{b}}$ and ${v}_{\text{max}}^{\mathrm{f}}$, which
are the theoretical predictions with a frictional bed and a
frictionless bed, respectively. The consideration of
${v}_{\text{max}}^{\mathrm{b}}$ and ${v}_{\text{max}}^{\mathrm{f}}$
reflects the influence of bed friction, slope angle and path
length. In addition, the effect of snow friction (*M*) and cohesion
(*β*) is also considered. The deposit height *h*_{d} is
scaled with the initial height of the snow slab, the snow properties
and the path length. Note two constant coefficients *κ*_{1} and *κ*_{2} are used to
account for other possible factors including snow compressive strength
and brittleness, where *κ*_{1} is dimensionless and *κ*_{2} has
a dimension of length. In this study, *κ*_{1} and *κ*_{2} are
0.2 and 200 m, respectively. According to
Eq. (6), when the friction *M* and cohesion *β* of
snow are high, the numerator is close to
${v}_{\text{max}}-{v}_{\text{max}}^{\mathrm{b}}$, and a zero numerator
indicates a maximum velocity close to the theoretical prediction
considering a rigid block sliding on a frictional bed. As shown in
Fig. 7, the data around the zero
line hint that the maximum velocity of the flows is chiefly
controlled by the friction between the flow and the bed. On the other
hand, when *M* and *β* tend to zero, ${v}_{\text{max}}^{*}$ approaches
${v}_{\text{max}}/{v}_{\text{max}}^{\mathrm{f}}$ in Eq. (6),
reflecting how close the maximum flow velocity is to the theoretical
prediction with a rigid block sliding on a frictionless
bed. Correspondingly, the cases with small *M* and *β* in
Fig. 7 reflect a maximum flow
velocity primarily governed by snow properties, instead of by the
frictional behavior between the flow and the bed. A representative
case is the cold dense flow in Fig. 3b,
where its maximum velocity is close to ${v}_{\text{max}}^{\mathrm{f}}$
as the flow is highly sheared. Furthermore, data below the zero line
are observed in Fig. 7,
corresponding to the cases where the snow box either stays on the
slope with limited displacement or slides down the slope with a
velocity that sometimes decreases (i.e., not a constant acceleration as
assumed in the calculation of ${v}_{\text{max}}^{\mathrm{b}}$).

### 3.3.2 Maximum velocity and *α*

The runout angle *α* is defined as $\mathit{\alpha}=\text{arctan}(H/L)$. *H* and *L* are total vertical drop and total
horizontal reach, respectively, determined based on the top point of
the release zone and the front of the final deposit
(Lied and Bakkehøi, 1980). Figures 8
and 9 show the relation between
*v*_{max} and *α*, including MPM data and real-measurement
data collected from McClung and Gauer (2018). For the five groups
of MPM simulations varying the slope angle *θ* and horizontal
length *L*_{0}, all of them largely follow a two-stage relation between
*v*_{max} and *α*: an initially decreasing
*v*_{max} and a subsequently constant *v*_{max} with the
increase in *α*. As demonstrated in
Figs. 8 and
9, the first stage mainly consists
of cases with low friction and cohesion, while the second stage is
chiefly composed of cases with high friction and cohesion. At the
first stage, a higher *v*_{max} leads to a longer runout
distance and thus a smaller *α*. This indicates the dominant
effect of *v*_{max} in controlling the runout distance, which
might be due to the positive correlation between the velocity and the
kinetic energy of a snow avalanche. Indeed, it has been recognized
that the runout distance is tightly related to the kinetic energy of
the flow upon its arrival at the deposition zone
(Sovilla et al., 2006). Obviously, the dominance of *v*_{max}
disappears at the second stage, as a similar *v*_{max} gives a
notably different *α*. The runout distance at this stage is
mainly affected by the deposition behavior, instead of by the flowing
behavior. For example, assuming a warm plug flow and a warm shear flow
sharing an identical *v*_{max} before they reach the deposition
zone, their runout distances can differ significantly, since the warm plug flow
may stop abruptly while the warm shear flow may gradually become
steady and have a relatively longer runout distance.

From Figs. 8 and
9, the effects of *θ* and
*L*_{0} on *v*_{max} are similar, as both of them have a positive
correlation with *v*_{max}. In addition, the slope angle
*θ* also influences the maximum runout angle as shown in
Fig. 8. The larger the slope
angle, the larger the maximum runout angle. This is due to the
definition of the runout angle *α*, which gives a maximum runout
angle close to the slope angle *θ*. The maximum runout angle is
reached when a flow stops on the slope with the modeled
configuration. With $\mathit{\theta}=\mathrm{30}{}^{\circ}$, several flows stay on the
slope and have *α*≈*θ*. All the flows with *θ*=40 and 50^{∘} go to the connecting-arc and deposition
zones, giving *α*<*θ*. Note that the runout angle *α* has
been correlated to the mean slope angle *β* in exploring the
runout distance of snow avalanches
(Lied and Bakkehøi, 1980; Barbolini et al., 2000; Delparte et al., 2008). As
ideal slopes (see Fig. 1) are adopted here, the
mean slope angle *β* is close to the slope angle *θ*. Indeed,
the positive correlation between the maximum runout angle and the
slope angle in Fig. 8 agrees with
the *α*−*β* model (Lied and Bakkehøi, 1980). The MPM results
in Figs. 8 and
9 generally fall into the range of
the real-measurement data from McClung and Gauer (2018). In
particular, the lower bound of *v*_{max} from the real
measurements is recovered with the MPM simulation. Note the case with
*v*_{max}=70 m s^{−1} serving as the upper bound of
the field data was a powder snow avalanche, whose behavior differs
significantly from the simulated dense snow avalanches. In addition, the path
length of the upper-bound case is significantly higher than the
adopted ones in the MPM simulations (McClung and Gauer, 2018). This
upper-bound case can indeed be captured with our MPM modeling by
varying the model setup but is not the focus here. It was reported by
McClung and Gauer (2018) that the runout angle has a negative
correlation with the maximum front velocity but with wide scatter as
observed from the blue dots in
Figs. 8 or
9. According to our sensitivity
study, the scatter might be a result of different terrain conditions
(e.g., slope angle), release positions (e.g., path length) and snow
properties. In addition, some data might be on the plateau stage where
the runout distance is governed by the deposition behavior instead of
the maximum front velocity.

Figure 10 demonstrates a unified trend with the dimensionless velocity
${v}_{\text{max}}^{*}$ in Eq. (6) and *α*^{*} as follows:

where *α*^{b} is calculated by assuming a sliding rigid
block. Referring to Fig. 1, the velocity of the
block increases from 0 to $\sqrt{\mathrm{2}{a}^{\mathrm{b}}l}$ as it slides down
from upstream to the end of the frictional slope. With an
assumption that the velocity of the block does not change before and
after it goes across the connecting-arc zone, its runout distance on
the deposition zone can be calculated, with an initial velocity of
$\sqrt{\mathrm{2}{a}^{\mathrm{b}}l}$, a constant acceleration of −*μ**g* and a
final velocity of 0. *α*^{b} can then be derived as
*α*^{b}=arctan*μ*. It is interesting to obtain
*α*^{b} that is solely dependent on the bed friction coefficient
*μ*. The scaled runout angle ${\mathit{\alpha}}^{*}=\mathrm{1}$ means a runout distance
fully consistent with the prediction using the sliding rigid block
theory, while ${\mathit{\alpha}}^{*}<\mathrm{1}$ and ${\mathit{\alpha}}^{*}>\mathrm{1}$ denote a runout distance
which is longer and shorter, respectively, than the predicted one with the sliding
rigid block theory. Indeed, the data with ${\mathit{\alpha}}^{*}<\mathrm{1}$
in Fig. 10 generally have low
friction and cohesion, which reasonably produce the longer runout
distances. On the contrary, the cases with ${\mathit{\alpha}}^{*}>\mathrm{1}$ are typically
more frictional and cohesive, leading to the shorter runout
distances. Note that the data close to the zero line in
Fig. 10 correspond to the cases
at the plateau stage in Figs. 8
and 9. As discussed in
Fig. 7, when the friction *M* and
cohesion *β* are high, a zero ${v}_{\text{max}}^{*}$ comes from a
maximum velocity *v*_{max} that approaches the theoretical
prediction ${v}_{\text{max}}^{\mathrm{b}}$ with consideration of a rigid
block sliding on a frictional bed, indicating the maximum velocity
*v*_{max} is dominated by the frictional behavior between the
flow and the bed. On the other hand, the *v*_{max} of the cases
far above and below the zero line in
Fig. 10 are governed by the snow
properties.

To testify the capability of the MPM modeling in capturing key dynamic
features (i.e., front velocity and position) of snow avalanches, five
reported real avalanches with different complex terrains are
simulated. All the cases are modeled in 2D, neglecting the variation
along the flow width direction. The adopted geometry of the terrains
is borrowed from the corresponding literatures. As no detailed snow properties of
the avalanches were measured and reported, the applied snow properties
in the MPM refer to the description of the snow type and snow
condition. In particular, three of the avalanches mainly consisted of
dry, loose and new snow, while the other two were chiefly composed of
wind-packed and settled old snow. Correspondingly, two groups of snow
properties are adopted as summarized in
Table 3. Based on the determined snow
densities (150 kg m^{−3} for new snow and 250 kg m^{−3} for old snow),
Young's modulus and the tensile strength can be estimated using the
relations from Gaume (2012) and Scapozza (2004). The
friction of the slope is calibrated according to the existing data of
the real
avalanches. Figures 11–15
show the MPM simulation results in comparison with the reported data
from the literatures. Particularly, the evolution of the scaled front
velocity is examined along the flow path (Gauer, 2014). The front
velocity from the field was obtained by means of Doppler radar devices
and photo analyses. Different measurement approaches may give
different velocities but are generally consistent with one another
(Rammer et al., 2007). The comparison basis between velocities from
numerical modeling and real measurements needs to be carefully
checked, as it is sometimes questionable
(Fischer et al., 2014; Rauter and Köhler, 2020). For example,
depth-averaged velocities from numerical modeling cannot be directly
compared to peak intensity velocities from Doppler radar measurements
(Rauter and Köhler, 2020). In
Figs. 11–15,
the front velocity from the MPM is determined as the approach velocity
(Rauter and Köhler, 2020), which is calculated from the evolution of
the front position with time and is assumed to be comparable with the
data from the different measurement techniques. Note that this
approach velocity has a different definition from the velocity of the
particles at the front of a flow, although their values are almost
identical in our simulations. The geometry of the terrain is denoted
by the dashed gray curves in
Figs. 11–15,
where the coordinates *x* and *y* are normalized with the vertical
drop height *H*_{0}. The red bands in
Figs. 11–14
denote measurement error.

Case I and II in Figs. 11 and 12 are two avalanches successively released at the northwest flank of the Weissfluh north ridge (Gubler et al., 1986; Gauer, 2014), whose velocity was measured with continuous-wave (CW) Doppler radar. The snow forming the first avalanche was dry and mostly loose new snow, which produced a powder cloud. In comparison, the second avalanche consisted of wind-packed snow, which led to blocky slab-type release. It is noticed that the consistency between the MPM results and the measured data is better for the second avalanche in Fig. 12. The underestimated maximum front velocity in Fig. 11 might be due to the challenge of capturing the powder cloud of the first avalanche with the MPM. The front velocity of a powder snow avalanche is normally obtained from the frontal dilute part, whose velocity can be higher than the dense core of the avalanche (Sovilla et al., 2015). In addition, the neglection of entrainment in the simplified MPM simulation may also contribute to the discrepancy in Fig. 11. It is suspected that the first release induced much more entrainment than the second one, considering the availability of the snow to be entrained.

Figure 13 shows the avalanche of Case III, which happened after strong snowfall (Gauer, 2014). There were no field observations due to the stormy weather. The velocity was measured by pulsed Doppler radar. The snow was conjectured to be dry or only slightly moist. The adopted snow properties in the MPM modeling refer to those of new snow, which are assumed to be identical to the snow in Case I as listed in Table 3. Figure 13 illustrates reasonable agreement between the MPM and the measured data, in terms of both the final front position and the maximum front velocity. During the flowing process, the MPM tends to underestimate the front velocity, which might be related to the dry nature of the snow as discussed in Case I (Fig. 11). Compared with cases I and II, the front velocity evolution of Case III fluctuates more, as the terrain is more irregular.

Case IV in Fig. 14 is a snow avalanche composed of snow cornices at the release position and settled old snow in the track (Gauer, 2014). The consistency between the MPM data and the estimated velocity from a series of timed photographs (Gauer, 2014) is satisfactory, except for the overestimated velocity at the beginning of the flow. The overestimated front velocity from the MPM is tightly related to the abruptly steepened slope at $x/{H}_{\mathrm{0}}\approx \mathrm{0.2}.$ The increase in front velocity in reality was not as sharp as the MPM result, which might be due to the effect of more entrainment especially at the steep part of the slope. Indeed, it was reported that the maximum velocity of a simulated snow avalanche without entrainment is higher than that with entrainment, given the same runout distance (Sovilla and Bartelt, 2002; Sovilla et al., 2007). Moreover, the measurement data are based on photo series with intervals of 1 s, while the time gaps of the MPM data are 1∕24 s. The maximum front velocity in reality might be lost within the 1 s interval of the measurement.

Figure 15 demonstrates the data of Case V, including those from real measurements and RAMMS (Christen et al., 2010) as well as from our MPM modeling. The front velocity from the field was calculated from timed photographs. The dry-snow avalanche in Case V was artificially triggered, and the development of the powder part was reported (Christen et al., 2010). To be consistent, the adopted properties for the dry snow in cases I and III are used for Case V. As shown in Fig. 15, both the RAMMS and MPM results show reasonable consistency with the real-measurement data. The calculated final front positions from RAMMS and the MPM are similar, while the maximum front velocity is underestimated and overestimated by RAMMS and the MPM, respectively. As discussed in Case IV, the MPM modeling does not take entrainment into account, which might be the reason for the overestimated front velocity.

This study explores the dynamics of snow avalanches with the material
point method (MPM) and an elastoplastic constitutive model. By virtue
of the capability of the MPM in simulating processes with large
deformations, fractures and collisions and coexistence of solid- and
fluid-like behaviors, a wide range of distinct snow avalanches with
diverse flow behaviors have been investigated. The reported four flow
regimes for dense snow avalanches from real observations have all been
captured from our MPM simulations, including cold shear, warm shear,
warm plug and slab sliding regimes. Moreover, in transition from cold
shear to warm shear flow regimes, flows with surges and small granules
are observed. The evolution of the avalanche front, the free-surface
shape and the vertical velocity profile shows distinct
characteristics for the different flow regimes, underpinning the
identification of flow regimes. In addition to the flow surface and the
shear behavior presented in this study, other features of the flow may
also be used to pinpoint the flow regimes, such as snow temperature
and liquid water content (Köhler et al., 2018). Although they are not
explicitly taken into account in this study, the changing snow
properties in our MPM modeling are capable of capturing the
characteristics of the different regimes. Furthermore, distinct
stopping mechanisms and maximum velocities were reported for the four
regimes (Köhler et al., 2018). For example, the cold dense regime was
identified by a starving stopping mechanism, where the flow deposits and
stops firstly from its tail and then towards its front. And the velocity of the
cold dense regime was reported to be smaller than
30 m s^{−1}. It is noticed that the simulated flow with the MPM
does not fully follow these descriptions, which might be due to the
idealized MPM setup and different terrain conditions.

We have systematically examined the effects of snow properties, slope
angle and path length on the flow and deposition behaviors of snow
avalanches, including the maximum flow velocity on the slope, the
runout angle and the avalanche deposit height. It is found that snow
friction and cohesion are closely related to the behavior of snow
avalanches. Low snow friction and cohesion give fluid-like behavior
and highly sheared flows, while high snow friction and cohesion lead
to solid-like flow with limited shearing. Both slope angle and path
length have a positive correlation with the maximum flow velocity on
the slope, while their effects on the deposit height are
trivial. Furthermore, unified trends have been obtained with
normalization of the maximum flow velocity, the deposit height and the
runout angle, revealing analogous physical rules under the different
conditions. Key controlling factors of *v*_{max} have been
identified, including the friction between the bed and the flow and the
geometry of the slope, as well as the snow properties. Depending on
snow properties, the runout angle is either controlled by the flow
behavior of a snow avalanche before its arrival at the deposition
zone or governed by its deposition behavior. It should be noted that
a wide range of material parameters have been adopted for the
sensitivity study. The combination of extreme flow properties leading
to very high velocity might not be realistic for snow avalanches. The
material parameters need to be carefully calibrated for investigation
of real snow avalanches.

The MPM modeling has been calibrated and tested through simulations of real snow avalanches on irregular terrains. The calculated avalanche front position and velocity from the MPM show reasonable agreement with the measurement data from the literature. The behavior of dense snow avalanches has been well recovered by the MPM. A discrepancy was observed particularly for avalanches which developed a powder cloud above the dense core, as the powder cloud has not been modeled here.

The presented research focuses on examining the flow regimes and flow dynamics of snow avalanches with idealized conditions, which is a preliminary study serving as the basis for investigating more realistic and complex snow avalanches. The 2D ideal slope with a constant inclination could be further changed to other shapes to be more realistic, such as a parabolic track. Although the 2D setups were used to efficiently conduct the systematic study including more than 1000 cases, it is fully possible to explore interesting cases with 3D MPM simulations (Gaume et al., 2019). Future studies will take into account real topography in 3D and recover the natural boundary conditions of snow avalanches. In addition, a new framework will need to be developed for investigating snow avalanches with a powder cloud, by considering a new constitutive law for the cloud and its interaction with both the dense core of snow avalanches and the air around the cloud (e.g., air friction). To further consider entrainment, the snow cover could be explicitly simulated with our model. This would however significantly increase the computational time. Alternatively, one could add a mass flux rate term to the mass balance equation, which considers the snow cover as a rigid boundary and estimates the entrained mass based on empirical and theoretical relations (Naaim et al., 2004; Issler and Pérez, 2011). Meanwhile, the momentum conservation needs to be adjusted to account for the momentum change in snow avalanches due to entrainment. Despite the assumptions and idealization applied in this study, it is demonstrated that the MPM model provides a promising pathway towards systematic and quantitative investigations of snow avalanche dynamics and flow regime transitions under the effects of snow mechanical properties and terrain geometries, which can improve our understanding of wet snow avalanches and offer analysis for avalanche dynamics with the influence of climate change.

The constitutive model adopted in this study perfectly satisfies the second law of thermodynamics. Following the derivation in Gaume et al. (2018), proving energy does not increase is equivalent to proving the plastic dissipation rate ${\dot{w}}^{P}(\mathit{X},t)$ is nonnegative. ${\dot{w}}^{P}$ can be computed as

where ** τ** is the Kirchhoff stress tensor,
𝓛

_{v}is the Lie derivative and

**b**

^{E}is the elastic right Cauchy–Green strain tensor. Since we use an associative flow rule, ${\mathcal{L}}_{v}{\mathbf{b}}^{\mathrm{E}}=-\mathrm{2}\dot{\mathit{\gamma}}\frac{\partial y}{\partial \mathit{\tau}}{\mathbf{b}}^{\mathrm{E}}$ (see Eq. 10 in Gaume et al., 2018), ${\dot{w}}^{P}$ can be expressed as

Recall Eq. (11) in Gaume et al. (2018) that $\dot{\mathit{\gamma}}\ge \mathrm{0}$. Furthermore, $\widehat{\mathit{\tau}}\cdot \frac{\partial y}{\partial \widehat{\mathit{\tau}}}\ge \mathrm{0}$ because the yield surface is a convex function of $\widehat{\mathit{\tau}}$ which includes the origin. Therefore ${\dot{w}}^{P}\ge \mathrm{0}$. Note this result holds for any isotropic plasticity model that has a convex yield function and an associative flow rule.

The evolution of kinetic and potential energy of the flows in the four typical flow regimes (i.e., cold dense, warm shear, warm plug, sliding slab) is demonstrated in Fig. A1. As expected, the potential energy of the flows initially decreases as the flows move down from the slope and then becomes steady after the flows stop. The kinetic energy of the flows firstly increases and then reduces until it vanishes. It is noticed that the kinetic energy of the sliding slab shows fluctuations in the deceleration phase, due to the interactions between the separating slabs in the flow after they reach the connecting-arc zone (see Supplement Video S1).

Figure A2 shows the dissipated energy of the flows in the four cases. The dissipated energy increases before it reaches the final steady state. The growth rate of the dissipated energy varies for the different flows as they have distinct flow behaviors. Nevertheless, the final energy dissipation does not show much difference for the different flows. This is because of the identical initial potential energy and the similar final potential energy of the flows.

The energy dissipation is contributed to by (1) the internal force of the material and (2) the external force on the material from the boundary slope. As illustrated in Fig. A3, in all the four cases, the dissipated energy from the boundary is much higher than that dissipated inside the material. This is consistent with the results in Gracia et al. (2019).

All the relevant data are available on Zenodo at https://doi.org/10.5281/zenodo.3965795 (Li et al., 2020a).

Videos of the avalanches presented in Fig. 2 are available on Zenodo at https://doi.org/10.5281/zenodo.3944698 (Li et al., 2020b).

The supplement related to this article is available online at: https://doi.org/10.5194/tc-14-3381-2020-supplement.

JG designed the study and obtained the funding. XL conducted the study, performed the simulations and wrote the paper under the supervision of JG. BS commented on the paper and provided guidance on the flow regime transitions. CJ developed the MPM tool, which has been used in this study.

The authors declare that they have no conflict of interests.

The first author acknowledges Lars Blatny for his support with the language editing of this paper.

This research has been supported by the Swiss National Science Foundation (grant no. PCEFP2_181227).

This paper was edited by Florent Dominé and reviewed by Frédéric Dufour and two anonymous referees.

Abdelrazek, A. M., Kimura, I., and Shimizu, Y.: Numerical simulation of a small-scale snow avalanche tests using non-Newtonian SPH model, Journal of Japan Society of Civil Engineers, 70, I_681–I_690, 2014. a

Abe, K. and Konagai, K.: Numerical simulation for runout process of debris flow using depth-averaged material point method, Soils Found., 56, 869–888, 2016. a

Ancey, C.: Snow avalanches, in: Oxford Research Encyclopedia of Natural Hazard Science, 2016. a

Barbolini, M., Gruber, U., Keylock, C., Naaim, M., and Savi, F.: Application of statistical and hydraulic-continuum dense-snow avalanche models to five real European sites, Cold Reg. Sci. Technol., 31, 133–149, 2000. a

Chandel, C., Srivastava, P. K., and Mahajan, P.: Determination of failure envelope for faceted snow through numerical simulations, Cold Reg. Sci. Technol., 116, 56–64, 2015. a

Choubin, B., Borji, M., Mosavi, A., Sajedi-Hosseini, F., Singh, V. P., and Shamshirband, S.: Snow avalanche hazard prediction using machine learning methods, J. Hydrol., 577, 123929, https://doi.org/10.1016/j.jhydrol.2019.123929, 2019. a

Christen, M., Kowalski, J., and Bartelt, P.: RAMMS: Numerical simulation of dense snow avalanches in three-dimensional terrain, Cold Reg. Sci. Technol., 63, 1–14, 2010. a, b

Delparte, D., Jamieson, B., and Waters, N.: Statistical runout modeling of snow avalanches using GIS in Glacier National Park, Canada, Cold Reg. Sci. Technol., 54, 183–192, 2008. a

Eglit, M., Yakubenko, A., and Zayko, J.: A review of Russian snow avalanche models – From analytical solutions to novel 3D models, Geosciences, 10, 77, https://doi.org/10.3390/geosciences10020077, 2020. a

Faug, T., Turnbull, B., and Gauer, P.: Looking beyond the powder/dense flow avalanche dichotomy, J. Geophys. Res.-Earth, 123, 1183–1186, 2018. a

Fischer, J.-T., Kowalski, J., and Pudasaini, S. P.: Topographic curvature effects in applied avalanche modeling, Cold Reg. Sci. Technol., 74, 21–30, 2012. a

Fischer, J.-T., Fromm, R., Gauer, P., and Sovilla, B.: Evaluation of probabilistic snow avalanche simulation ensembles with Doppler radar observations, Cold Reg. Sci. Technol., 97, 151–158, 2014. a

Fischer, J.-T., Kofler, A., Fellin, W., Granig, M., and Kleemayr, K.: Multivariate parameter optimization for computational snow avalanche simulation, J. Glaciol., 61, 875–888, 2015. a

Gauer, P.: Comparison of avalanche front velocity measurements and implications for avalanche models, Cold Reg. Sci. Technol., 97, 132–150, 2014. a, b, c, d, e

Gauer, P., Issler, D., Lied, K., Kristensen, K., and Sandersen, F.: On snow avalanche flow regimes: Inferences from observations and measurements, in: Proceedings Whistler 2008 International Snow Science Workshop, 21–27 September 2008. a

Gaume, J.: Prédétermination des hauteurs de départ d'avalanches. Modélisation combinée statistique-mécanique, PhD thesis, Grenoble, 2012. a, b

Gaume, J., Chambon, G., and Naaim, M.: Quasistatic to inertial transition in granular materials and the role of fluctuations, Phys. Rev. E, 84, 051304, https://doi.org/10.1103/PhysRevE.84.051304, 2011. a

Gaume, J., Gast, T., Teran, J., van Herwijnen, A., and Jiang, C.: Dynamic anticrack propagation in snow, Nat. Commun., 9, 1–10, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Gaume, J., van Herwijnen, A., Gast, T., Teran, J., and Jiang, C.: Investigating the release and flow of snow avalanches at the slope-scale using a unified model based on the material point method, Cold Reg. Sci. Technol., 168, 102847, https://doi.org/10.1016/j.coldregions.2019.102847, 2019. a, b, c

Gracia, F., Villard, P., and Richefeu, V.: Comparison of two numerical approaches (DEM and MPM) applied to unsteady flow, Computational Particle Mechanics, 6, 591–609, 2019. a

Gubler, H., Hiller, M., Klausegger, G., and Suter, U.: Messungen an Fliesslawinen. Zwischenbericht, in: Internal Report 41, Swiss Federal Institute for Snow and Avalanche Research, 1986. a

Hagenmuller, P., Chambon, G., and Naaim, M.: Microstructure-based modeling of snow mechanics: a discrete element approach, The Cryosphere, 9, 1969–1982, https://doi.org/10.5194/tc-9-1969-2015, 2015. a

Issler, D. and Pérez, M. P.: Interplay of entrainment and rheology in snow avalanches: a numerical study, Ann. Glaciol., 52, 143–147, 2011. a

Jiang, C., Schroeder, C., Selle, A., Teran, J., and Stomakhin, A.: The affine particle-in-cell method, ACM T. Graphic., 34, 1–10, 2015. a

Jiang, C., Schroeder, C., Teran, J., Stomakhin, A., and Selle, A.: The material point method for simulating continuum materials, in: ACM SIGGRAPH 2016 Courses, 1–52, 2016. a

Kern, M., Bartelt, P., Sovilla, B., and Buser, O.: Measured shear rates in large dry and wet snow avalanches, J. Glaciol., 55, 327–338, 2009. a

Keshari, A. K., Satapathy, D. P., and Kumar, A.: The influence of vertical density and velocity distributions on snow avalanche runout, Ann. Glaciol., 51, 200–206, 2010. a

Klár, G., Gast, T., Pradhana, A., Fu, C., Schroeder, C., Jiang, C., and Teran, J.: Drucker-prager elastoplasticity for sand animation, ACM T. Graphic., 35, 1–12, 2016. a

Köhler, A., McElwaine, J., and Sovilla, B.: GEODAR Data and the flow regimes of snow avalanches, J. Geophys. Res.-Earth, 123, 1272–1294, 2018. a, b, c, d, e, f, g

Larsson, S., Prieto, J. M. R., Gustafsson, G., Häggblad, H.-Å., and Jonsén, P.: The particle finite element method for transient granular material flow: modelling and validation, Computational Particle Mechanics, 1–21, https://doi.org/10.1007/s40571-020-00317-6, 2020. a

Li, X., Sovilla, B., Jiang, C., and Gaume, J.: Supplementary data for “The mechanical origin of snow avalanche dynamics and flow regime transitions”, Zenodo, https://doi.org/10.5281/zenodo.3965795, 2020a. a

Li, X., Sovilla, B., Jiang, C., and Gaume, J.: Supplementary videos for “The mechanical origin of snow avalanche dynamics and flow regime transitions”, Zenodo, https://doi.org/10.5281/zenodo.3944698, 2020b. a

Lied, K. and Bakkehøi, K.: Empirical calculations of snow–avalanche run–out distance based on topographic parameters, J. Glaciol., 26, 165–177, 1980. a, b, c

Mast, C. M., Arduino, P., Miller, G. R., and Mackenzie-Helnwein, P.: Avalanche and landslide simulation using the material point method: flow dynamics and force interaction with structures, Computat. Geosc., 18, 817–830, 2014. a, b

McClung, D. and Gauer, P.: Maximum frontal speeds, alpha angles and deposit volumes of flowing snow avalanches, Cold Reg. Sci. Technol., 153, 78–85, 2018. a, b, c, d

Naaim, M., Naaim-Bouvet, F., Faug, T., and Bouchet, A.: Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects, Cold Reg. Sci. Technol., 39, 193–204, 2004. a

Naaim, M., Durand, Y., Eckert, N., and Chambon, G.: Dense avalanche friction coefficients: influence of physical properties of snow, J. Glaciol., 59, 771–782, 2013. a

Papakrivopoulos, V.: Development and preliminary evaluation of the main features of the Particle Finite Element Method (PFEM) for solid mechanics, Master's thesis, Delft University of Technology, 2018. a

Rammer, L., Kern, M., Gruber, U., and Tiefenbacher, F.: Comparison of avalanche-velocity measurements by means of pulsed Doppler radar, continuous wave radar and optical methods, Cold Reg. Sci. Technol., 50, 35–54, 2007. a

Rauter, M. and Köhler, A.: Constraints on Entrainment and Deposition Models in Avalanche Simulations from High-Resolution Radar Data, Geosciences, 10, 9, https://doi.org/10.3390/geosciences10010009, 2020. a, b, c

Rauter, M., Kofler, A., Huber, A., and Fellin, W.: faSavageHutterFOAM 1.0: depth-integrated simulation of dense snow avalanches on natural terrain with OpenFOAM, Geosci. Model Dev., 11, 2923–2939, https://doi.org/10.5194/gmd-11-2923-2018, 2018. a

Raymond, S. J., Jones, B., and Williams, J. R.: A strategy to couple the material point method (MPM) and smoothed particle hydrodynamics (SPH) computational techniques, Computational Particle Mechanics, 5, 49–58, 2018. a

Reiweger, I., Gaume, J., and Schweizer, J.: A new mixed-mode failure criterion for weak snowpack layers, Geophys. Res. Lett., 42, 1427–1432, 2015. a

Salazar, F., Irazábal, J., Larese, A., and Oñate, E.: Numerical modelling of landslide-generated waves with the particle finite element method (PFEM) and a non-Newtonian flow model, Int. J. Numer. Anal. Meth., 40, 809–826, 2016. a

Scapozza, C.: Entwicklung eines dichte-und temperaturabhängigen Stoffgesetzes zur Beschreibung des visko-elastischen Verhaltens von Schnee, PhD thesis, ETH Zürich, 2004. a

Simo, J. and Meschke, G.: A new class of algorithms for classical plasticity extended to finite strains. Application to geomaterials, Comput. Mech., 11, 253–278, 1993. a

Simo, J. C.: Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory, Comput. Method. Appl. M., 99, 61–112, 1992. a

Soga, K., Alonso, E., Yerro, A., Kumar, K., and Bandara, S.: Trends in large-deformation analysis of landslide mass movements with particular emphasis on the material point method, Géotechnique, 66, 248–273, 2016. a, b

Sovilla, B. and Bartelt, P.: Observations and modelling of snow avalanche entrainment, Nat. Hazards Earth Syst. Sci., 2, 169–179, https://doi.org/10.5194/nhess-2-169-2002, 2002. a

Sovilla, B., Burlando, P., and Bartelt, P.: Field experiments and numerical modeling of mass entrainment in snow avalanches, J. Geophys. Res.-Earth, 111, F03007, https://doi.org/10.1029/2005JF000391, 2006. a

Sovilla, B., Margreth, S., and Bartelt, P.: On snow entrainment in avalanche dynamics calculations, Cold Reg. Sci. Technol., 47, 69–79, 2007. a

Sovilla, B., McElwaine, J. N., Schaer, M., and Vallet, J.: Variation of deposition depth with slope angle in snow avalanches: Measurements from Vallée de la Sionne, J. Geophys. Res.-Earth, 115, F02016, https://doi.org/10.1029/2009JF001390, 2010. a

Sovilla, B., McElwaine, J. N., and Louge, M. Y.: The structure of powder snow avalanches, C. R. Phys., 16, 97–104, 2015. a

Srivastava, P., Chandel, C., and Mahajan, P.: Micromechanical modeling of elastic and strength properties of snow, SLAM3-Slab Avalanche Multiscale Mechanical Modeling, 3–5, 2017. a

Steinkogler, W., Sovilla, B., and Lehning, M.: Influence of snow cover properties on avalanche dynamics, Cold Reg. Sci. Technol., 97, 121–131, 2014. a

Stomakhin, A., Schroeder, C., Chai, L., Teran, J., and Selle, A.: A material point method for snow simulation, ACM T. Graphic., 32, 1–10, 2013. a, b, c, d

Techel, F., Jarry, F., Kronthaler, G., Mitterer, S., Nairz, P., Pavšek, M., Valt, M., and Darms, G.: Avalanche fatalities in the European Alps: long-term trends and statistics, Geogr. Helv., 71, 147–159, https://doi.org/10.5194/gh-71-147-2016, 2016. a

Zabala, F. and Alonso, E.: Progressive failure of Aznalcóllar dam using the material point method, Géotechnique, 61, 795–808, 2011. a

- Abstract
- Introduction
- Methodology
- Snow avalanches on ideal slopes
- Snow avalanches on irregular terrains
- Conclusions and discussion
- Appendix A: Energy evolution of the flows in the four typical flow regimes
- Data availability
- Video supplement
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Methodology
- Snow avalanches on ideal slopes
- Snow avalanches on irregular terrains
- Conclusions and discussion
- Appendix A: Energy evolution of the flows in the four typical flow regimes
- Data availability
- Video supplement
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References