Detrainment and braking of snow avalanches interacting with forests

Mountain forests provide natural protection against avalanches. They can both prevent avalanche formation in release zones and reduce avalanche mobility in runout areas. Although the braking effect of forests has been previously explored through global statistical analyses on documented avalanches, little is known about the mechanism of snow detrainment in forests for small and medium avalanches. In this study, we investigate the detrainment and braking of snow avalanches in forested terrain, by performing three-dimensional simulations using the material point method (MPM) and a large-strain elastoplastic snow constitutive model based on critical state soil mechanics. First, the snow internal friction is evaluated using existing field measurements based on the detrainment mass, showing the feasibility of the numerical framework and offering a reference case for further exploration of different snow types. Then, we systematically investigate the influence of snow properties and forest parameters on avalanche characteristics. Our results suggest that for both the cold and warm snow parameterized in our simulations, the detrainment mass decreases with the square of the avalanche front velocity before it reaches a plateau value. Furthermore, the detrainment mass significantly depends on snow properties. It can be as much as 10 times larger for warm snow compared to cold snow. By examining the effect of forest configurations, it is found that forest density and tree diameter have cubic and square relations with the detrainment mass, respectively. The outcomes of this study may contribute to the development of improved formulations of avalanche–forest interaction models in popular operational simulation tools and thus improve hazard assessment for alpine geophysical mass flows in forested terrain.

1 Introduction

The expansion of human activity in mountains has increased the risk associated with snow avalanches (De Biagi et al., 2012), which threaten infrastructures and human lives. Forests can mitigate snow avalanche hazards without expensive construction of heavy protective measures like concrete dams. Hence, it is crucial to understand how forests affect the avalanche dynamic behaviour. A forest can have two protective effects against avalanches. It stabilizes the snowpack in the release area (Viglietti et al., 2009) and reduces avalanche mobility in the runout zone. The second aspect has often been neglected because of the field observation that large avalanches can destroy forests with no significant deceleration (Bartelt and Stöckli, 2001). However, it has been shown that for small to medium avalanches (e.g. <10 000 m³), forests are normally not destroyed and can notably reduce the avalanche run-out distance (Teich et al., 2012a; Perzl et al., 2021).

Predominant numerical tools for modelling snow avalanches use two-dimensional approaches, like depth-averaged avalanche dynamics models (Christen et al., 2010). These methods are computationally efficient and a powerful tool for hazard mapping but cannot directly capture the interaction between individual trees and an avalanche. Therefore, many studies have worked on modelling the effect of a forest indirectly. There are mainly two approaches to simulate the protective effect of a forest against avalanches: the frictional approach and the detrainment approach (Feistl et al., 2012). The first accounts for the breakage of trees and debris entrainment by increasing the friction parameters of the Voellmy law (Voellmy, 1955) compared to open unforested terrain (Gruber and Bartelt, 2007; Takeuchi et al., 2018). The second approach is based on field observations and has a more solid physical meaning than the friction model. The mass and momentum of the snow stopped behind the trees are directly removed from the flow, which naturally leads to deceleration and run-out shortening of the avalanche. Assuming that no mass is entrained in the forest area, the detrainment rate ${\dot{Q}}_{\mathrm{d}}=-\frac{\mathrm{1}}{\mathit{\rho }}{\dot{M}}_{\mathrm{d}}$ is added to the mass balance, with M_d the detrainment mass, which corresponds to the mean mass stopped in the forest per unit of area. The rate of detrainment is quantified with a detrainment coefficient K (Feistl et al., 2014), which links the temporal derivative of the detrainment mass with V, the depth-averaged velocity of the avalanche as follows: $\frac{\mathrm{d}{M}_{\mathrm{d}}}{\mathrm{d}t}=-\frac{K}{\left|\right|\mathit{V}\left|\right|}$.

The two methods can both empirically recover the reduction in mass and momentum of an avalanche in an efficient way but may encounter difficulties and give unsatisfactory predictions in practical applications. According to Teich et al. (2012b), the friction approach is not effective for modelling small and medium avalanches in forested terrain, as the run-out distance was overestimated even by using the smallest turbulent friction coefficient. The detrainment approach highly relies on the data from field avalanches, with which the detrainment coefficient K can be determined to apply the detrainment function in numerical modelling (Teich et al., 2014). Based on existing studies, both the forest type and the snow properties have a crucial effect on the detrainment coefficient K (Teich et al., 2012b, 2014; Brožová et al., 2020). However, there is no systematic investigation on the controlling factors of the snow detrainment by forests under well-controlled conditions.

This study aims to quantify the amount of snow caught in forests and explore the mechanism of the detrainment, with comprehensive consideration of different avalanche features and forest configurations. We use a three-dimensional material point method (3D MPM), by which the fractures, collisions, and large deformations involved in snow avalanches can be well captured (Gaume et al., 2018, 2019; Li et al., 2020, 2021). More importantly, individual trees and their interactions with snow avalanches can be explicitly modelled without relying on empirical laws. Compared to the popular two-dimensional tools for modelling snow avalanches, the 3D MPM can fully resolve flow variations in all dimensions. Based on the simulations in this study, we can better understand the forest–avalanche interaction, reveal laws governing the effect of the key influencing factors on snow detrainment, and offer a basis for systematic calibration of the detrainment approach (Feistl et al., 2014) for operational purposes.

2 Methodology and setup

2.1 Numerical method

2.1.1 The material point method (MPM)

Our simulations are performed using the material point method (MPM), which is a hybrid Eulerian–Lagrangian method suitable to deal with large strains (Stomakhin et al., 2013). In MPM, Lagrangian particles are used to carry the information of position, velocity, and deformation gradient, while the Eulerian grid is used to compute the equations of motion and for updating the particle states. The governing equations are based on mass and momentum conservation as follows:

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

where t is time, v is velocity, ρ is density, g is the gravitational acceleration, and σ is the Cauchy stress expressed by the elastoplastic constitutive relation:

$$\begin{array}{}\text{(3)}& \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).

In this study, the transfers between particles and grid use the affine particle-in-cell (APIC) method (Jiang et al., 2015). Compared to the particle-in-cell (PIC) and fluid implicit particle (FLIP) techniques, the APIC approach allows both the stable removal of the numerical dissipation and the preservation of angular momentum in addition to linear momentum.

Figure 1 Model setup of MPM simulations of snow avalanches impacting trees. The left and right are side and front views respectively. Five arrangements of trees in the right include (a) one tree, (b) three trees, (c) regularly staggered trees, (d) regularly aligned trees, and (e) random trees.

2.1.2 Large-strain elastoplastic model

For modelling snow, we use the modified Cam–Clay model with associative flow rule developed by Gaume et al. (2018). Given a stress τ, a mean effective stress and a deviatoric stress can be respectively obtained as $p=-\frac{\mathrm{1}}{\mathrm{3}}\mathrm{tr}\left(\mathit{\tau }\right)$ and $\mathbf{s}=\mathit{\tau }+p\mathbf{I}$, where I is the identity matrix. The Mises equivalent stress can then be derived as $q=\sqrt{\frac{\mathrm{3}}{\mathrm{2}}\mathbf{s}:\mathbf{s}}$. The yield surface in the space of the q–p is defined as

$$\begin{array}{}\text{(4)}& 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}$$

where β is the ratio between tensile and compressive strength and represents the cohesion, p₀ is the consolidation pressure, and M is the slope of the cohesionless critical state line, which denotes the internal friction of snow. For $y(p,q)\le \mathrm{0}$, the material follows the Hooke law (St Venant–Kirchhoff–Hencky strain), otherwise it behaves plastically. The hardening and softening of the material are modelled by respectively expanding and shrinking the yield surface according to the law

$$\begin{array}{}\text{(5)}& p_{\mathrm{0}}={\displaystyle \frac{E}{\mathrm{3}(\mathrm{1}-\mathrm{2}\mathit{\nu })}}\sinh \left(\mathit{\zeta }\max \left(-{\mathit{ϵ}}_{\mathrm{v}}^{\mathrm{p}},\mathrm{0}\right)\right),\end{array}$$

where ζ is the hardening factor, E the Young modulus, ν the Poisson ratio, and ${\mathit{ϵ}}_{\mathrm{v}}^{\mathrm{p}}=\log \left(\det \left({\mathbf{F}}^{\mathrm{P}}\right)\right)$ the volumetric plastic deformation.

2.2 Model setup

Different simplified setups are designed to test the protective role of different forest structures with various types of snow. As shown in the left of Fig. 1, the terrain is composed of an ideal slope with a constant inclination of $\mathit{\theta }=\mathrm{30}{}^{\circ }$ and a horizontal deposition zone. The bed friction coefficient of the terrain is assumed to be 0.5. The snow sample is initially placed at the top of the slope at a height of H₀=50.6 m with a prescribed initial velocity of v₀. We define h₀=1 m, l₀=20 m, and w₀=20 m as respectively the height, length, and width of the released snow. For all the snow types simulated in this study, they share the same Young modulus (E=3 MPa), Poisson ratio (ν=0.3), and initial density (ρ₀=200 kg m⁻³) (Li et al., 2020). As the material point method is a continuous approach, each material point corresponds to a piece of snow; therefore the snow density is the bulk density, including air and snow. Five types of tree arrangement are designed as shown in the right of Fig. 1, including (a) one tree, (b) three trees with an identical diameter and at the same elevation, (c) a regularly staggered forest composed of trees in a staggered arrangement, (d) a regularly aligned forest composed of trees in an aligned arrangement, and (e) a random forest in which trees are randomly located. All the trees are initially at the downstream of the released snow. In the setup of (c), (d), and (e), the width of the forests (along z) is larger than the maximum avalanche width, and the length of the forests along the flow direction is 40 m.

In addition to the different tree arrangements, the forest density ρ_forest (trees ha⁻¹) and the tree diameter d are also varied to study their effect on snow detrainment. The forest density in this study is defined as the number of trees per hectare. In the case of a regularly staggered forest, the forest density is controlled by the spacing between the trees e (Fig. 1). In this study, we assume that the trees are unbreakable and the trunk surface is rough. Therefore, the trees are modelled as rigid obstacles with a no-slip boundary condition. The height of the trees is set to 10 m so that the avalanches do not overtop the trees.

In order to both model the trees and shorten the computation time well, we use an Eulerian grid size of dx=0.15 m if it is not specified and eight particles per grid. The time step is constrained by the Courant–Friedrichs–Lewy (CFL) condition and the elastic wave speed to guarantee the simulation stability. The simulation data are exported every $\mathrm{1}/\mathrm{12}$ s.

3 Calibration with a documented case

To calibrate the friction of snow and verify other parameters (e.g. snow density, Young modulus, mesh size) adopted in our numerical modelling, we simulate the snow detrainment due to a group of three trees according to the field observation in Feistl et al. (2014). This case with three trees is selected for the calibration as the setup is simple, and the field measurement data are available for comparison. Based on the real condition (Feistl et al., 2014), the diameter of trees is set to 1 m, and the spacing between trees is 0.33 m in our simulation. The simulation results are quantitatively compared with the field data in terms of the volume and height of the deposited mass (i.e. detrainment) behind the trees.

By assuming the shape of the detrainment mass as a wedge, Feistl et al. (2014) calculated the detrainment volume as

$$\begin{array}{}\text{(6)}& W={\displaystyle \frac{l{h}_{\mathrm{w}}^{\mathrm{2}}}{\mathrm{2}\tan \left(\mathit{\theta }\right)}},\end{array}$$

with l the width of the group of trees, h_w the wedge height, and θ the slope angle. In our simulations, the detrainment volume is obtained by summarizing the volume of all the particles in the detrainment mass without the assumption of the wedge shape. The height of the detrainment mass is calculated as the vertical distance from its free surface to the tree foot (as illustrated in Fig. 2). As summarized in Table 1, five simulation cases with different snow friction coefficient M were conducted, while all other parameters are fixed (e.g. the tension compression ratio β=0.3, the hardening factor ζ=1, the initial consolidation pressure $p_{\mathrm{0}}^{\mathrm{ini}}$ = 30 kPa). With the reduction of snow friction from 1.5 to 0.8, the volume and height of the detrainment mass decrease, and the case with M=0.8 has a good agreement with the field data. In addition, Fig. 2 shows that the shape of the detrainment mass does not differ much with the change of snow friction. The simulation with M=0.8 gives a similar profile of the detrainment mass as that observed from the field (Fig. 2d and e). This calibration case demonstrates the good feasibility of our numerical tool in capturing the key features of detrainment and serves as a reference case for the study on the effect of snow properties in the following section.

Table 1 Deposition height and volume of snow in simulations with different M and in the field observation.

Figure 2 Cross-sectional view (passing through the centre of the middle tree) of the snow stopped behind a group of three trees for different values of M: (a) M=1.2, (b) M=1, and (c) M=0.8. Perspective view of the stopped mass (d) in the simulation with M=0.8 and d=1 m and (e) from the field observation (Feistl et al., 2014).

4 Results

A qualitative study of the forest–avalanche interaction will be first presented where different avalanche characteristics will be introduced (Sect. 4.1). Then, the main forest factors influencing the detrainment mass in the case of a regular staggered forest will be discussed (Sect. 4.2), and a newly proposed unique law for estimation of the detrainment mass will be suggested (Sect. 4.3). Finally, we will show an energy investigation involving the braking process as well as the influence of snow properties and forest structure (Sect. 4.4).

4.1 Influence of avalanche features

4.1.1 Snow properties at the forest scale

In this study, different snow properties are adopted in the MPM modelling, to study their influence on the braking effect of the forest. The snow properties are defined by the parameters of the modified Cam–Clay model (Gaume et al., 2018): the friction coefficient M, the tension compression ratio β, the hardening factor ζ, and the initial consolidation pressure $p_{\mathrm{0}}^{\mathrm{ini}}$. These parameters have their physical basis and can be determined according to physical properties of snow and/or parametric study (Li et al., 2020). In particular $\mathit{\beta }{p}_{\mathrm{0}}^{\mathrm{ini}}$ denotes the isotropic tensile strength, M is linked to the internal friction angle of snow ϕ as follows $\mathit{\varphi }={\sin }^{-\mathrm{1}}\left(\frac{\mathrm{3}M}{\mathrm{6}+M}\right)$ (Sadrekarimi and Olson, 2011), and the hardening factor reflects how fast the load increases with the displacement in the plastic stage. According to our systematic study on the effect of snow properties on the avalanche behaviour (Li et al., 2020), it was found that the tensile strength $\mathit{\beta }{p}_{\mathrm{0}}^{\mathrm{ini}}$ and βM consistently increase from cold to warm avalanches, which suggests that these two terms control the different snow behaviours.

Figure 3 Flow profile for three different flow regimes: (a) Case 1, (b) Case 2, and (c) Case 3, for a regular staggered forest at t=10 s, e=8 m, d=1 m, and v₀=10 m s⁻¹.

Table 2 Snow properties adopted in the MPM modelling of the flows for three typical flow regimes.

Based on the calibration case reported in Sect. 3, $\mathit{\beta }{p}_{\mathrm{0}}^{\mathrm{ini}}=\mathrm{9}$ kPa and βM=0.24 (Case 2 in Table 2), we can obtain different snow types by changing $\mathit{\beta }{p}_{\mathrm{0}}^{\mathrm{ini}}$ and βM. Firstly, to reproduce the behaviour of a cold snow (Case 1 in Table 2), we decrease the tensile strength $\mathit{\beta }{p}_{\mathrm{0}}^{\mathrm{ini}}=\mathrm{0}$ kPa and βM=0, which leads to a cohesionless and low-friction snow. And to capture the behaviour of a wet and warm snow, $\mathit{\beta }{p}_{\mathrm{0}}^{\mathrm{ini}}=\mathrm{9}$ kPa and βM=0.36 (Case 3 in Table 2) are adopted, giving a relatively cohesive and frictional snow with granules and blocks.

As illustrated in Fig. 3, a regular staggered forest is used for the three cases in Table 2. It is apparent that the snow properties have a significant effect on the accumulated snow behind the trees. Indeed, the total amount of stopped snow increases with the internal friction and the cohesion of snow. According to our simulations, some arches due to the jamming effect (Feistl et al., 2014) appear in Case 3, especially upstream of the avalanche. This occurrence of the arches promotes the stabilization of snow on the slope and is tightly related to the high internal friction and cohesion of snow. The jamming effect and the formation of the arches observed in the avalanche–forest interaction are similar to the behaviour of a granular flow in a two-dimensional hopper (Lai et al., 2001).

Figure 4 Profile of accumulated snow for different flow regimes: (a) Case 1, (b) Case 2, and (c) Case 3, regular staggered forest at t=10 s, e=8 m, d=1 m, and v₀=10 m s⁻¹.

Independent of the snow properties, there are a lot of defections as the flow front is split into many branches (Fig. 3), in contrast to the flow without forest obtained from our simulation. Figure 3 demonstrates that the more cohesionless the snow is, the more the snow is laterally spread in the forest and gets around the trees easily. As shown by Luong et al. (2020), for the cohesionless flow regime (Case 1), the snow on the side of the flow moves very fast. After being initially deflected, this lateral snow moves in a diagonal channel without collision, even if a part of this effect is due to the arrangement of trees like a Galton board.

In addition to the lateral motion of the avalanches, the length and height of the stopped snow are further analysed with the side view as illustrated in Fig. 4. The insets at the upper corner of Fig. 4 show the stopped mass (in red) at the uppermost trees. It is interesting to observe that the height of the wedge of stopped mass does not vary much with the different snow properties. Meanwhile, the length of the snow wedge in the flow direction changes notably. With the more cohesive and frictional snow from Case 1 to Case 3, the length of the wedge increases, and its free surface becomes more horizontal (inset c). According to Fig. 4, the snow properties also affect the evolution of stopped mass along the flow direction. While the shapes of the snow wedges appear to be similar along the slope in Cases 1 and 2, the profile of the snow changes significantly in Case 3. As blocks and granules exist in Case 3, there is more snow stopped by the trees upstream. It is noticed that the uppermost trees in Case 3 have less stopped mass than the second uppermost trees, mainly due to the severe initial impact of the avalanche on the uppermost trees.

4.1.2 Avalanche velocity

To define the detrainment mass stopped by the trees, a criteria based on the evolution of particle velocity is used. If a snow particle has a velocity smaller than 0.5 m s⁻¹ throughout the flowing process, it is considered stopped and detrained. This definition borrows the concept by Feistl et al. (2014), where the snow is directly removed from the avalanche mass once it is stopped.

Figure 5 Evolution of the detrainment mass per unit of area of snow with the velocity (regular staggered forest, e=8 m).

The effect of avalanche velocity on detrainment snow is studied by imposing an initial velocity v₀ ranging from 0 to 24 m s⁻¹. Figure 5 shows the evolution of the detrainment mass with the front velocity at the entrance of the forest for the different snow properties (in Sect. 4.1.1). The detrainment mass per unit area, M_d, is calculated by dividing the total detrainment mass by the total area of influence. The latter is obtained by multiplying the number of trees with snow accumulation by the area of influence of each tree (e²).

Firstly, a strong effect of snow properties on the stopped snow is observed from Fig. 5. For instance, the total amount of snow stopped in Case 3 can be 10 times larger than that in Case 1 when the velocity is small. However, for each type of snow the same behaviour is observed, and the detrainment mass decreases when the velocity increases and tends to a constant value for high velocities. The initial negative correlation with the low velocities is in agreement with the field observation that a lower velocity leads to more stopped snow (Feistl et al., 2014). Therefore, we propose that the detrainment mass decreases with velocity square according to

$$\begin{array}{}\text{(7)}& M_{\mathrm{d}}\left(V_{\mathrm{f}}\right)={\displaystyle \frac{p_{\mathrm{1}}+p_{\mathrm{2}}{V}_{\mathrm{f}}^{\mathrm{2}}}{V_{\mathrm{f}}^{\mathrm{2}}}},\end{array}$$

with $V_{\mathrm{f}}=\sqrt{v_{\mathrm{xf}}^{\mathrm{2}}+v_{\mathrm{yf}}^{\mathrm{2}}+v_{\mathrm{zf}}^{\mathrm{2}}}$ the front velocity of the avalanche when it enters the forest. The front velocity is calculated based on the velocity of the first 15 % of the particle flow to eliminate the effect of scattering particles at the avalanche front. The parameters of this model (p₁(kg s⁻²), p₂(kg m⁻²)) (Fig. 5) are identified based on least-squares method (Table 3), and both depend on the snow properties. p₂ represents the detrainment mass reached for very high velocities (Eq. 7); it is the mass stopped independently of the velocity.

Similar to the effect of velocity on the detrainment mass, Kyburz et al. (2020) showed that the impact pressure increases with the velocity square. This increase in the impact pressure causes the breakage of snow during its collision with the trees, which leads to a lower volume of snow stuck behind the trees. The detrainment mass for the highest velocity is rather due to the type of snow; this part of snow is initially compacted against the trunk at the front of the avalanche and is not impacted by the velocity.

4.2 Influence of forest parameters

4.2.1 Tree diameter

The influence of tree diameter is studied for the one-tree configuration with a mesh size of 0.1 m for a better modelling of thin trees. In the case of a single tree, we use the mass stopped m_d (kg), not normalized with the area around the tree as the area is difficult to define. For all the snow properties, both the maximum and the final masses stopped m_d have the same evolution trend with the tree diameter (Fig. 6). Feistl et al. (2014) proposed the following cubic model for the detrainment volume W induced by single-tree interaction:

$$\begin{array}{}\text{(8)}& W={\displaystyle \frac{d^{\mathrm{3}}}{\mathrm{12}\tan \left(\mathit{\theta }\right)\tan (\mathit{\delta }/\mathrm{2}{)}^{\mathrm{2}}}},\end{array}$$

with d the diameter of the tree, θ the slope angle, and δ the top wedge angle of the pyramid formed by stopped snow.

Figure 6 Evolution of the detrainment mass observed for the one-tree arrangement and the fits with the tree diameter for three flow regimes: Case 1, Case 2, and Case 3 with a front velocity of respectively 12.5, 10.9, and 10.75 m s⁻¹. “Maximum” refers to the maximum mass stored behind the tree and “final” refers to the final mass stored. A slope of 30^∘ and a top wedge angle of 60^∘ (from measurements Feistl et al., 2014) are used for the theoretical model (Eq. 8). The removed point denotes a special case not considered in the proposed square relation in Eq. (9) since the entire avalanche in this case is stopped due to the low flow velocity and the high tree diameter. Please note that it is a coincidence that Case 1 maximum and Case 2 final agree well.

Comparing the simulated final mass with the previous analytical proposal of Feistl et al. (2014), we can notice that this model is close to the flow in Case 2 calibrated for a tree diameter of 1 m. It proves that the calibration is consistent. However, the numerical result shows that the cubic model in Eq. (8) is not verified. This observation is in accordance with the analysis of the shape of the snow wedge behind the tree; indeed the height of the wedge is not linearly influenced by the tree diameter, whereas the diameter has a large effect on the length and the width of the wedge. Thus, the evolution of the stopped mass with the tree diameter follows a square law for each type of snow.

Figure 7 Distribution of the maximum detrainment mass behind each tree during the avalanche (regular staggered forest, d=1 m, e=5.5 m, snow properties: Case 2): (a) v₀=10 m s⁻¹; (b) v₀=22 m s⁻¹.

While we dispose of the evolution of the detrainment mass with the tree diameter for the one-tree arrangement, in the case of an arrangement of several trees, the spatial distribution of the detrainment mass shows that trees interact with each other (Fig. 7). Comparing to the one-tree result, more snow is stopped for the first row than for the one-tree case, and for the lower rows the detrainment mass can be either underestimated or overestimated. This is not only a geometric effect, because the pattern of the spatial distribution of detrainment mass behind each tree depends on the velocity (Fig. 7). Thus, it reveals the difficulty to predict the distribution of detrainment mass. For a better understanding of the phenomenon, the study must be carried out at the forest scale, and study of the forest density would help to consider the collective behaviour and interaction of the trees.

Figure 8 Evolution of the detrainment mass per unit of area with the tree diameter for different types of snow and front velocity (regular staggered forest, e=8 m).

For the study of a forest, the influence of tree diameter is examined by keeping the same forest density and changing the basal surface area per hectare $\mathrm{\Phi }=\frac{A_{\mathrm{trunk}}}{A_{\mathrm{forest}}}$, where A_trunk is the total area occupied by the trunks and A_forest the total forest area.

Figure 8 depicts that the same behaviour is observed for different velocities and snow properties, and the trend is similar to the case of one tree. These results suggest that the detrainment mass evolves according to a square law with the tree diameter (Eq. 9), with p₃ and p₄ (1 m⁻²) evaluated based on the least-squares method (Table 3).

$$\begin{array}{}\text{(9)}& M_{\mathrm{d}}\left(V_{\mathrm{f},\mathrm{d}}\right)=M_{\mathrm{d}}\left(V_{\mathrm{f}}\right)\left(p_{\mathrm{3}}+p_{\mathrm{4}}{d}^{\mathrm{2}}\right)\end{array}$$

4.2.2 Forest density

The effect of forest density is studied with a constant tree diameter of 0.6 m in Fig. 9a and a fixed stand density index (SDI) (Reineke, 1933) in Fig. 9b. The range of the adopted SDI from 400 to 1347 trees ha⁻¹ in Fig. 9a corresponds to normal to very dense forests (Abegg et al., 2020). Figure 9a shows the evolution of the detrainment mass under the effect of the forest density, for different front velocities. In previous models (Feistl, 2015), it was assumed that the increase in the forest density only increases the number of trees/obstacles and consequently the number of wedges, which led to a linear law between the mass stopped and the forest density. However, our results suggest that the forest density has an influence on the stopped mass behind each individual tree. Indeed, compared to a forest with a lower density, a denser forest leads to more detrainment not only because there are more trees but also due to more stopped mass for each tree as snow jamming occurs. This increased detrainment due to the collective behaviour of the trees is observed for all the different velocities simulated. Based on the obtained simulation data, the cubic law is the best-fitting model; therefore the proposed law in Eq. (9) can be improved with consideration of the effect of forest density according to Eq. (10) below:

$$\begin{array}{}\text{(10)}& M_{\mathrm{d}}(V_{\mathrm{f},\mathrm{d}},{\mathit{\rho }}_{\mathrm{forest}})=M_{\mathrm{d}}\left(V_{\mathrm{f},\mathrm{d}}\right)\left(p_{\mathrm{5}}+p_{\mathrm{6}}{\mathit{\rho }}_{\mathrm{forest}}^{\mathrm{3}}\right),\end{array}$$

with p₅ and p₆ ((trees ha⁻¹)⁻³) evaluated based on the least-squares method (Table 3).

Figure 9 (a) Evolution of the detrainment mass per unit of area with the forest density for different front velocities (regular staggered forest, snow properties: Case 2, d=0.6 m). (b) Evolution of the detrainment mass per unit of area with the tree diameter for a constant stand density index of 925 trees ha⁻¹ (regular staggered forest, v₀=10 m s⁻¹, snow properties: Case 2).

The evolution of the detrainment mass with the tree diameter is shown in Fig. 9b, where the stand density index is fixed at 925 trees ha⁻¹, and the forest density is changing. When the tree diameter is small (d<0.6 m), the detrainment mass decreases with the tree diameter, due to the reduction of the forest density and the number of trees. Indeed, fewer trees give a weaker collective effect of the forest and thus stop less snow. A different trend is observed when the tree diameter becomes larger (d>0.6 m), indicating the competing roles played by the collective effect of the trees and the individual tree effect. The larger the tree diameter, the larger the stopped mass behind each individual tree, which changes the decreasing trend when the tree diameter is small. In addition to the simulation results, the prediction of the detrainment mass with the proposed model in Eq. (10) is shown in Fig. 9b. As observed, the model prediction agrees well with the simulation data, showing the good performance of the model and the possibility to use the stand density index in the model.

4.3 Model for detrainment mass

Based on the systematic study of the influence of the avalanche features and forest parameters, a model for the detrainment mass per unit of area can be proposed for a regular staggered forest (Eq. 11).

Table 3 Parameter value of the law (11) for different snow properties (Table 2).

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{M}_{\mathrm{d}}& (V_{\mathrm{f}},d,{\mathit{\rho }}_{\mathrm{forest}})=\\ & {\displaystyle \frac{p_{\mathrm{1}}+p_{\mathrm{2}}{V}_{\mathrm{f}}^{\mathrm{2}}}{V_{\mathrm{f}}^{\mathrm{2}}}}\left(p_{\mathrm{3}}+p_{\mathrm{4}}{d}^{\mathrm{2}}\right)\left(p_{\mathrm{5}}+p_{\mathrm{6}}{\mathit{\rho }}_{\mathrm{forest}}^{\mathrm{3}}\right),\end{array}\end{array}$$

with V_f the frontal velocity of the avalanche when it reaches the forest, ρ_forest (trees ha⁻¹) the number of trees per hectare, d the tree diameter, and $p_{[\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5},\mathrm{6}]}$ the model parameters given in Table 3.

Figure 10 depicts the evolution of the detrainment mass predicted from Eq. (11) and that from MPM simulations, showing that the developed law in Eq. (11) predicts the detrainment mass well. As expected, Case 1 and Case 3 have the smallest and largest detrainment mass respectively.

Figure 10 Evolution of the detrainment mass predicted with the model (Eq. 11) and with the observation; the coefficient of determination r² for the model prediction on the detrainment mass is 0.9912.

4.4 Energy analysis of the avalanche–forest interaction

4.4.1 Evolution of energy and energy dissipation

After the evaluation of the detrainment mass, it is needed to quantify the effect of the detrainment process in the braking effect of forest against avalanches. To understand the physical process of the braking effect, the analysis of the evolution of the kinetic and potential energy with and without forest is made (Fig. 11). As expected, we observe when the avalanche enters the forest that the kinetic energy decreases while without forest the kinetic energy continues to increase steadily. An important difference of kinetic energy between the simulation with and without forest when the avalanche front reaches the bottom of the slope is observed, which confirms the significant braking effect of forest. In addition, in terms of potential energy, whereas the potential energy without forest decreases linearly as a function of time, with forest, due to the mass which stays on the slope, the potential energy decreases more slowly with time. Consequently at the end of the simulation, the potential energy does not vanish due to the detrainment mass stopped by the forest.

Figure 11 Temporal evolution of the kinetic and potential energy without forest and with a regular staggered forest (Case 2, v₀=6 m s⁻¹).

To understand the kinetic energy loss with the forest, two physical processes are considered, namely detrainment and dissipation. First, a detrainment behaviour accounts for the stopped snow mass behind the trees, whose potential energy is not converted to kinetic energy. Considering the avalanche potential energy with and without forest, E_pf and E_p0 respectively, the change of the potential energy due to the detrainment with the presence of the forest E_pd can be defined as

$$\begin{array}{}\text{(12)}& E_{\mathrm{pd}}=E_{\mathrm{pf}}-E_{\mathrm{p}\mathrm{0}}.\end{array}$$

Secondly, the dissipation process is associated with avalanche–tree interactions including collision and friction, which lead to plastic dissipation and random fluctuation. This dissipation due to the trees is defined as ${\stackrel{\mathrm{̃}}{\mathrm{E}}}_{\mathrm{f}}$, whose calculation is detailed in Appendix A.

Figure 12 Temporal evolution of the detrainment energy and dissipation due to the forest (Case 2, v₀=6 m s⁻¹).

Figure 12 shows the energy change due to the detrainment and dissipation processes (i.e. E_pd and ${\stackrel{\mathrm{̃}}{\mathrm{E}}}_{\mathrm{f}}$), which contribute to the reduction of the avalanche kinetic energy with the forest. Taking the avalanche with forest as a reference, when its front enters the forest, the energy loss due to detrainment and dissipation starts to grow (at around 2.3 s in Fig. 12). The increasing rate of ${\stackrel{\mathrm{̃}}{\mathrm{E}}}_{\mathrm{f}}$ tends to decrease with time, as the velocity of the flow decreases (see the kinetic energy in Fig. 11) and the dissipation by ground friction reduces. In contrast, the growth rate of E_pd increases with time. Indeed, after small wedges are formed behind and around the trees, the spacing between trees decreases, which leads to faster blockage of incoming snow particles. It should be noted that the forest enhances energy loss due to detrainment and dissipation but meanwhile reduces the avalanche velocity. Hence, the presence of the forest leads to lower energy dissipation due to ground friction, as it is positively correlated with the avalanche velocity.

4.4.2 Influence of the forest structure and the snow properties on the energy loss

The same experiment is carried out for three types of snow (Table 2) and three types of forest (regular aligned, regular staggered, and random) with the same forest density and tree diameter, in order to analyse the separate processes of detrainment and dissipation. Figure 13 depicts the detrainment and dissipation energies when the front of the avalanche reaches the bottom of the slope, where the front is defined as the first 1 % of the avalanche. It reveals that the detrainment represents a larger part of the energy loss compared to the dissipation. The detrainment increases with the cohesion and the internal friction of snow, whereas the dissipation due to forest changes slightly with the snow properties. For the three types of snow, the dissipation is higher for the regular staggered forest, due to the Galton board arrangement for which collisions and deflections are more frequent. The most notable difference for the different forest arrangements in Fig. 13 is that much lower detrainment is observed for the aligned forest, indicating that much less snow is stopped. Indeed, when trees are regularly aligned, a large part of the flow passes through the space between trees without collisions (Fig. 14). These results suggest that in the case of a regular forest, the staggered arrangement should be privileged for a better protective effect. The detrainment energy when the avalanche reaches the bottom of the slope is similar between the random and the regular staggered arrangements.

Figure 13 Energy detrainment and dissipation for different snow properties (1: Case 1; 2: Case 2; 3: Case 3 with M=1) and for three types of forest structure (regular aligned, regular staggered, and random).

Given the same snow properties (e.g. Case 2), the final stopped mass with the aligned, staggered, and random arrangements is 2147, 3002, and 13061 kg respectively. The same trend is observed for other considered snow properties in this study (Cases 1 and 3). Compared to the regular aligned forest, the regular staggered forest stops more snow due to a smaller spacing (along the diagonal direction) between the trees (Fig. 14a and b). The random arrangement gives the highest stopped mass, as each tree cluster serves as a dam and collectively blocks a large amount of snow (Fig. 14c). Consequently, in terms of the final run-out distance, the avalanche in the random forest travels a shorter distance than that in the regular staggered forest, which suggests that the random one has a higher protective effect.

Figure 14 Flow height of the avalanche in a forest with (a) staggered arrangement, (b) aligned arrangement, and (c) random arrangement, (Case 2, v₀=6 m s⁻¹).

5 Discussion and conclusions

In this study, we have explored the avalanche–forest interactions with the material point method (MPM) and an elastoplastic constitutive model. The aim is to highlight the dependencies of the detrainment mass and understand the physical process of the braking effect of forest with respect to avalanche dynamics. We have systematically examined the effect of snow properties, by studying both cold and warm snow. It is found that snow friction and cohesion have a large effect on detrainment mass behind the trees. Low snow friction and cohesion give fluid-like behaviour and highly sheared flows, and consequently little snow is stopped. When snow becomes more frictional and cohesive, large wedges and arches appear behind the trees, which contribute to the large amount of snow stopped. Due to the scarcity of the snow triaxial test (Scapozza and Bartelt, 2003), the snow parameters in our study are determined and calibrated according to their physical basis and previously conducted parametric study. Nevertheless in the future, these parameters could be evaluated based on extensive triaxial tests for different snow types.

Interestingly, our analysis of the effect of avalanche velocity reveals that the detrainment mass decreases with the velocity and reaches a plateau value for high velocities (typically higher than 15 m s⁻¹). While the detrainment mass linearly decreases with increasing velocity in Feistl et al. (2014), our numerical results suggest a decrease with velocity square. Knowing the relation between the detrainment mass per unit of area M_d and the detrainment coefficient $K=\frac{M_{\mathrm{d}}.V_{\mathrm{f}}^{\mathrm{2}}.\cos \left(\mathit{\theta }\right)}{l_{\mathrm{avalanche}}}$ (Feistl, 2015), the detrainment rate can be assessed, with V_f the front velocity, θ the slope angle, and l_avalanche the avalanche length. For an open forest or old forests, Feistl et al. (2014) proposed a value of K=10 Pa independently of the snow properties and the front velocity. With our simulation results, for Case 1, we obtain a value of K=5.31 for a velocity of 7.4 m s⁻¹ and K=32.4 for a velocity of 24 m s⁻¹, while for Case 2, we obtain K=10.1 for a velocity of 7.7 m s⁻¹ and K=86.6 for a velocity of 24.7 m s⁻¹. This underlines a high dependency of the detrainment coefficient on the snow properties and on the front velocity of the avalanche. The comparison with the previous value (Feistl et al., 2014) of K for an open forest suggests that it is only valid for low velocities. When the avalanche velocity is high, the plateau stage of the detrainment mass obtained from our study is higher than the decreasing stage by Feistl et al. (2014). Consequently, an implementation of our proposed velocity-based model would lead to higher detrainment mass, which results in lower runout distance and smaller velocity.

In addition, we have examined different forest arrangements, by varying the forest density and the tree diameter, which both can easily be identified from a forest inventory or remote sensing (Bebi et al., 2016). It is found that the detrainment mass increases with the cube of forest density and the square of tree diameter. Contrary to the geometric model proposed by Feistl et al. (2014), we found that the height of the wedge does not vary linearly with the tree diameter. This comprehensive numerical investigation allowed us to propose a unified law for the detrainment mass (Eq. 11). Therefore, it enables an estimation of the K value as shown before, which may serve as an input to 2D depth-averaged methods to better predict the detrainment mass.

From an energy point of view, our results suggest that the detrainment has the largest influence on the energy loss, which claims the use of the detrainment approach. However, the additional dissipation due to forest (${\stackrel{\mathrm{̃}}{\mathrm{E}}}_{\mathrm{f}}$) is not negligible, although the detrainment approach removes not only the potential energy of the removed particles but also their kinetic energy. In view of the kinetic energy levels compared to the dissipation added with forest, the kinetic energy of the particles removed seems too low to represent the dissipation added with forest. This result implies that an increase in the bed friction in combination with the detrainment approach could be justified to take into account the increase in energy dissipation linked to the effect of forest.

The presented research focuses on examining the interaction between forest and avalanches with idealized conditions, which allow new outlooks for further studies. Here, trees are modelled as rigid obstacles and with a no-slip boundary condition. Hence, this study is only valid for small to medium avalanches where the forest is not destroyed, and the trees act as obstacles. Secondly taking into account the true friction coefficient of the tree would allow us, as suggested by Teich et al. (2014), to define a law for different species. In the future, adding a parameter to account for forest type, crown cover and surface roughness based on remote sensing data (Brožová et al., 2020) could help the operation of the proposed model and improve the evaluation of detrainment in Bayesian networks for risk assessment (Stritih, 2021).

Furthermore, slope angle can affect the detrainment mass, as it can significantly vary with the front velocity V_f. However, its effect on the geometry of the wedges is closely related to other factors like avalanche velocity and bed friction and therefore is difficult to be quantified independently. Future studies on the relation between the effect of slope angle and that of other factors will need to be conducted to introduce slope angle in our proposed model for the detrainment mass. In addition, the consideration of a slope with a constant inclination is an ideal condition. This could be further changed to any other shape to be more realistic, and other tree arrangements could be explored as well. It would also be interesting to study longer forest with larger released volumes, to obtain longer permanent regimes; this would however significantly increase the computational time.

Moreover, this study focuses only on the avalanche dynamics and not on the release in forest area; therefore the conclusions on the effect of forest arrangement should be carefully interpreted in practical forest design for avalanche mitigation. As the avalanche release could be greatly affected by the forest as well, it is interesting to investigate the stabilization of snowpack under the effect of trees using MPM simulations in the future.

Despite the assumptions and idealization applied in this study, it highlights the main factors influencing the detrainment mass in a unique law depending on snow mechanical properties, front velocity, and forest parameters. This outcome of this work can be applied not only for calibration of depth-averaged models used operationally (e.g. Christen et al., 2010), but also for forest management. Moreover, according to energy analysis, a good calibration of the detrainment approach seems to be not enough to model all processes driving avalanche–forest interaction. This suggests that an approach combining detrainment and friction increase could be more appropriate. Finally, although this study focused on snow avalanches, the methods used and the general relationships found are also relevant to other geophysical mass flows interacting with forests and can pose strong impacts on hazard assessment and risk management.

Appendix A: Calculation of energy dissipation due to forest

By defining Δe the additional dissipation at each time step, it is possible to separate it into two parts (Li et al., 2020): the dissipation inside the flow w_int and through the boundary bed w_b, and in the case of forest, the dissipation ${\stackrel{\mathrm{̃}}{\mathrm{e}}}_{\mathrm{f}}$ is added. Hence, the total dissipation in the case without forest is $\mathrm{\Delta }{e}_{\mathrm{0}}=w_{{\mathrm{b}}_{\mathrm{0}}}+w_{{\mathrm{int}}_{\mathrm{0}}}\approx w_{{\mathrm{b}}_{\mathrm{0}}}$, as the internal dissipation can be negligible compared to the bed friction dissipation (Li et al., 2020). Similarly, in the forest case $\mathrm{\Delta }{e}_{\mathrm{f}}=w_{{\mathrm{b}}_{\mathrm{f}}}+w_{{\mathrm{int}}_{\mathrm{f}}}+{\stackrel{\mathrm{̃}}{\mathrm{e}}}_{\mathrm{f}}\approx w_{{\mathrm{b}}_{\mathrm{f}}}+{\stackrel{\mathrm{̃}}{\mathrm{e}}}_{\mathrm{f}}$.

As the velocity of the avalanche is reduced with forest, the frictional dissipation with the ground decreases. To evaluate the dissipation added with forest, we can not simply subtract the dissipation obtained with and without forest. Therefore we use Eq. (A1), which allows us to calculate ${\stackrel{\mathrm{̃}}{\mathrm{E}}}_{\mathrm{f}}$ based on our simulations. Indeed, by replacing the expression of Δe₀ and Δe_f (Eq. A2), and by assuming that the bed dissipation w_b depends linearly on the velocity (i.e $\frac{w_{{\mathrm{b}}_{\mathrm{f}}}(i.\mathrm{\Delta }t)}{v_{\mathrm{f}}(i.\mathrm{\Delta }t)}\approx \frac{w_{{\mathrm{b}}_{\mathrm{0}}}(i.\mathrm{\Delta }t)}{v_{\mathrm{0}}(i.\mathrm{\Delta }t)}$) (Eq. A2), the definition of the dissipation added due to forest is obtained (Eq. A3). So at the time step n, with Δt the time step, it is possible to obtain the total dissipation added by the forest ${\stackrel{\mathrm{̃}}{\mathrm{E}}}_{\mathrm{f}}$.

$$\begin{array}{}\text{(A1)}& {\displaystyle }\sum _{i=\mathrm{0}}^{n-\mathrm{1}}\left({\displaystyle \frac{\mathrm{\Delta }{\mathrm{e}}_{\mathrm{f}}(i.\mathrm{\Delta }t)}{v_{\mathrm{f}}(i.\mathrm{\Delta }t)}}-{\displaystyle \frac{\mathrm{\Delta }{e}_{\mathrm{0}}(i.\mathrm{\Delta }t)}{v_{\mathrm{0}}(i.\mathrm{\Delta }t)}}\right).v_{\mathrm{f}}(i.\mathrm{\Delta }t)\text{(A2)}& {\displaystyle }=\sum _{i=\mathrm{0}}^{n-\mathrm{1}}({\displaystyle \frac{w_{{\mathrm{b}}_{\mathrm{f}}}(i.\mathrm{\Delta }t)+{\stackrel{\mathrm{̃}}{\mathrm{e}}}_{\mathrm{f}}(i.\mathrm{\Delta }t)}{v_{\mathrm{f}}(i.\mathrm{\Delta }t)}}-{\displaystyle \frac{w_{{\mathrm{b}}_{\mathrm{0}}}(i.\mathrm{\Delta }t)}{v_{\mathrm{0}}(i.\mathrm{\Delta }t)}}).v_{\mathrm{f}}(i.\mathrm{\Delta }t)\text{(A3)}& {\displaystyle }\approx \sum _{i=\mathrm{0}}^{n-\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathrm{e}}}_{\mathrm{f}}(i.\mathrm{\Delta }t)\text{(A4)}& {\displaystyle }={\stackrel{\mathrm{̃}}{\mathrm{E}}}_{\mathrm{f}}(n.\mathrm{\Delta }t)\end{array}$$