The mechanical origin of snow avalanche dynamics and flow regime transitions

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.

1 Introduction

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 Methodology

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:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{D}\mathit{\rho }}{\mathrm{D}t}}+\mathit{\rho }\mathrm{\nabla }\cdot \mathit{v}=\mathrm{0},\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\rho }{\displaystyle \frac{\mathrm{D}\mathit{v}}{\mathrm{D}t}}=\mathrm{\nabla }\cdot \mathit{\sigma }+\mathit{\rho }\mathit{g},\end{array}$$

where ρ is density, t is time, v is velocity, σ is the Cauchy stress and g 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.

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:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\mathit{\sigma }={\displaystyle \frac{\mathrm{1}}{J}}{\displaystyle \frac{\partial \mathrm{\Psi }}{\partial {\mathbf{F}}_{\mathrm{E}}}}{\mathbf{F}}_{\mathrm{E}}^T,\end{array}$$

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:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }y(p,q)=(\mathrm{1}+\mathrm{2}\mathit{\beta })q^{\mathrm{2}}+M^{\mathrm{2}}(p+\mathit{\beta }{p}_{\mathrm{0}})(p-p_{\mathrm{0}}).\end{array}$$

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}\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₀ is the consolidation pressure and denotes the isotropic compressive strength. βp₀ 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{ϵ}}_{\mathrm{v}}^{\mathrm{p}}$, hardening or softening is implemented by expanding or shrinking the yield surface according to the following hardening law:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{p}_{\mathrm{0}}=K\text{sinh}\left(\mathit{\xi }\text{max}\left(-{\mathit{ϵ}}_{\mathrm{v}}^{\mathrm{p}},\mathrm{0}\right)\right),\end{array}$$

where K is the bulk modulus and ξ is the hardening factor. Under compression (${\dot{\mathit{ϵ}}}_{\mathrm{v}}^{\mathrm{p}}<\mathrm{0}$), p₀ increases, leading to hardening and promoting compaction. Under tension (${\dot{\mathit{ϵ}}}_{\mathrm{v}}^{\mathrm{p}}>\mathrm{0}$), p₀ 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).

Figure 1 Model setup for MPM modeling of snow avalanches on an ideal slope.

Table 1 Parameters adopted in the MPM simulations of snow avalanches on ideal slopes.

^* M values include 0.1, 0.5, 1.0 and 1.5. β values include 0.0, 0.3, 0.6 and 1.0. ξ values include 0.1, 0.5, 1.0 and 10.0. $p_{\mathrm{0}}^{\text{ini}}$ values include 3, 12, 21 and 30 kPa.

3 Snow avalanches on ideal slopes

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₀ 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₀. When θ is varied in groups I, II and III, the horizontal length L₀ is fixed and the drop height H₀ is adjusted as listed in Table 1. Alternatively, one could fix the drop height H₀ and change the horizontal length L₀. Note L₀∕h₀, L₀∕l₀ and L₀∕r are kept constant when L₀ is varied in groups II, IV and V by changing h₀, l₀ and r accordingly. An increased L₀ leads to the scale-up of the setup, resulting in the rise in the drop height H₀. 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.

Figure 2 (a) Flows in four typical flow regimes captured at t=8.3 s in the MPM simulations (Table 1, Group II). From top to bottom: cold dense, warm shear, warm plug and sliding slab. The cold dense flow is scaled up to be 3 times higher along the bed-normal direction for better visualization. (b) A flow with surges and small granules in transition from cold dense regime to warm shear regime at t=8.3 s. The color denotes velocity as in (a). (c) The early stage (t=5.5 s) of the warm shear flow in (a). Videos of the simulations are provided in the Supplement.

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.

Table 2 Snow properties adopted in the MPM modeling of the flows in the four typical flow regimes.

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=gsinθ 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.

Figure 3 Evolution of (a) front and (b) front velocity of the flows in the four typical regimes (Table 1, Group II).

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}\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.

Figure 4 Evolution of free-surface shape (a–d) and velocity profile at x = 50 m (e–h) for the flows in the different regimes (Table 1, Group II). (a, e) Cold dense; (b, f) warm shear; (c, g) warm plug; (d, h) sliding slab. The free surface of the cold dense flow in (a) is scaled up 15 times along the bed-normal direction for better visualization.

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₀, under the effects of the slope angle θ and the horizontal length L₀ (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₀).

Figure 5 Evolution of the maximum velocity with the normalized deposit height for varying slope angles θ.

Figure 6 Evolution of the maximum velocity with the normalized deposit height for different horizontal lengths L₀.

Figure 7 Unified relation between the scaled maximum velocity and the scaled deposit height.

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₀ (see Fig. 1), the drop height H₀ 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₀ 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₀, 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:

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{v}_{\text{max}}^{*}={\displaystyle \frac{v_{\text{max}}-v_{\text{max}}^{\mathrm{b}}\left[\mathrm{1}-e^{-M(\mathrm{1}+\mathit{\beta })/{\mathit{\kappa }}_{\mathrm{1}}}\right]}{v_{\text{max}}^{\mathrm{f}}}},\text{(7)}& {\displaystyle }& {\displaystyle }{h}^{*}={\displaystyle \frac{h_{\mathrm{d}}}{h_{\mathrm{0}}}}\left[\mathrm{1}-e^{-M(\mathrm{1}+\mathit{\beta }){\mathit{\kappa }}_{\mathrm{2}}/l}\right].\end{array}$$

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 κ₁ and κ₂ are used to account for other possible factors including snow compressive strength and brittleness, where κ₁ is dimensionless and κ₂ has a dimension of length. In this study, κ₁ and κ₂ 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₀, 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.

Figure 8 Evolution of the maximum velocity with α for varying slope angles θ.

Figure 9 Evolution of the maximum velocity with α for different horizontal lengths L₀.

From Figs. 8 and 9, the effects of θ and L₀ 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⁻¹ 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 Unified relation between the scaled maximum velocity and the scaled α.

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

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathit{\alpha }}^{*}={\displaystyle \frac{\mathit{\alpha }}{{\mathit{\alpha }}^{\mathrm{b}}}},\end{array}$$

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.

4 Snow avalanches on irregular terrains

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⁻³ for new snow and 250 kg m⁻³ 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₀. The red bands in Figs. 11–14 denote measurement error.

Table 3 Adopted parameters in the five MPM simulations of snow avalanches on real terrains.

^* Calibrated parameter.

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 11 Front velocity distribution along the flow path for Case I: Weissfluh north ridge, 12 March 1982, a1 (Davos, Switzerland). Drop height H₀=236 m.

Figure 12 Front velocity distribution along the flow path for Case II: Weissfluh north ridge, 12 March 1982, a2 (Davos, Switzerland). Drop height H₀=177 m.

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.

Figure 13 Front velocity distribution along the flow path for Case III: Himmelegg, 14 February 1990 (Arlberg, Austria). Drop height H₀=352 m.

Figure 14 Front velocity distribution along the flow path for Case IV: Ryggfonn, 2 May 2006 (Stryn, Norway). Drop height H₀=303 m.

Figure 15 Front velocity distribution along the flow path for Case V: Vallée de la Sion (VdlS), 31 January 2003 (Sion, Switzerland). Drop height H₀=1246 m.

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.

5 Conclusions and discussion

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⁻¹. 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.

Appendix A: Energy evolution of the flows in the four typical flow regimes

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

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }{\dot{w}}^P=-\mathit{\tau }:{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathcal{L}}_v{\mathbf{b}}^{\mathrm{E}}\right){\left({\mathbf{b}}^{\mathrm{E}}\right)}^{-\mathrm{1}},\end{array}$$

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

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }{\dot{w}}^P=\mathit{\tau }:\dot{\mathit{\gamma }}{\displaystyle \frac{\partial y}{\partial \mathit{\tau }}}=\dot{\mathit{\gamma }}\widehat{\mathit{\tau }}\cdot {\displaystyle \frac{\partial y}{\partial \widehat{\mathit{\tau }}}}.\end{array}$$

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.

Figure A1 Evolution of potential and kinetic energy of the flows in the four typical flow regimes. (a) Cold dense; (b) warm shear; (c) warm plug; (d) sliding slab.

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.

Figure A2 Evolution of dissipated energy of the flows in the four typical flow regimes.

Figure A3 Energy dissipation inside the flow and through the boundary bed in the flows with the different flow regimes.

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).