Potential and challenges of depth-resolved three-dimensional MPM simulations: a case study of the 2019 ‘Salezer’ snow avalanche in Davos

Abstract

Avalanche modeling is an essential tool to assess snow avalanche hazard. Today, most popular numerical approaches adopt depth-averaged equations. These methods are computationally efficient but limited in capturing processes occurring in the flow depth direction, e.g. erosion or deposition, which are often considered using ad hoc parameterizations or neglected completely. However, processes such as snow erosion, can crucially influence flow dynamics and run-out and are often not negligible. We address these issues by using a new three-dimensional (3-D) model, based on the material point method and finite strain elastoplasticity. To assess the possibilities and challenges associated with these highly detailed but computationally expensive calculations, we simulated the ‘Salezer’ snow avalanche that released in Davos, Switzerland in 2019. To reproduce the event in our simulations, we use the release areas mapped in a photogrammetric drone survey and estimate the snow conditions on the day of the event. We compare macroscopic features, such as flow outline and snow deposition of the simulated avalanche to field observations. An in-depth analysis of transient 3-D flow structures at the avalanche head not only demonstrates the degree of physical detail in the model, but also highlights challenges which still need to be addressed.

1. Introduction

The goal to understand and predict the dynamic behavior of snow avalanches is often to mitigate avalanche danger by estimating e.g. the avalanche flow velocity and run-out to plan suitable countermeasures. Models that are widely used today are based on the analogy between avalanches and floods, implementing a set of depth-averaged equations derived from the Navier–Stokes equations. Due to the depth integration these so-called Saint-Venant models involve a number of complex assumptions about the flow dynamics, as well as ad hoc parametrizations and conceptual or empirical models of dynamic processes. The parameters involved in these parametrizations and empirical laws need to be calibrated from historical events (e.g. McDougall and Hungr, 2004; Zugliani and Rosatti, 2021), which implies strong limits in their predictive capacity. Especially the Coulomb and turbulent friction parameters in the widely used Voellmy rheological model play an important role governing the run-out distance of the avalanche, but are not comparable to a physically measurable mechanical property of snow.

Furthermore, Saint-Venant-like models suffer from shortcomings when simulating flows on steep or high-curvature terrain due to the depth integration of the flow equations. While improvements were made to resolve this issue (e.g. Gray and others, 1999; Pudasaini and Hutter, 2003), all depth-averaged models used today inherently suffer from this limitation to some degree. This limitation has special importance for snow avalanches, because snow avalanches mostly occur in steep alpine terrain. Moreover, in times of a warming climate, the frequency and characteristic of the snow avalanche hazard is transforming as well (Lazar and Williams, 2008; Castebrunet and others, 2014; Naaim and others, 2016). This creates the need for more physics-based models with better predictive capacity compared to the most commonly used depth-averaged models, which are often calibrated using historic data.

In the past few decades, novel high-resolution measurement technologies (e.g. Thibert and others, 2008; Kern and others, 2009; Sovilla and others, 2015; Gauer and Kristensen, 2016; Köhler and others, 2018) were used to improve the physical understanding of the processes governing snow avalanche dynamics. The interpretation of the measurements and the development of numerical models, which consider the analogy between snow avalanches and granular flows (e.g. Sampl and Granig, 2009; Li and others, 2021; Ligneau and others, 2022), and reproduce the experimental observations in ever greater physical detail, allow for an even deeper insight into the dynamic flow processes.

One particular modeling approach, namely the material point method (MPM), received increased attention because it performs well in simulating the large material deformations, as well as aggregation and fracturing processes that materials undergo in geophysical mass flows, including snow avalanches. A recent implementation of MPM, which has been developed to simulate crack initiation and propagation for snow avalanche release (Stomakhin and others, 2013; Gaume and others, 2018), proved also to perform well in simulating the dynamics of hazardous geophysical mass movements in general (e.g. Gaume and Puzrin, 2021; Wolper and others, 2021; Cicoira and others, 2022; Li and others, 2022b). The respective studies demonstrate that this MPM model is able to reproduce dynamic flow processes such as snow entertainment, surges, flow regime transitions and snow granulation (Li and others, 2020, 2021, 2022b), which are important to study the dynamics of snow avalanches.

In the present study, we further push the boundaries of the mesh resolution and the physical detail, which can be achieved in fully three-dimensional (3-D) simulations of snow avalanches over an explicitly simulated erodible bed, and thus, exploring the possibilities and limitations of this up-to-date 3-D MPM snow avalanche model. Fully 3-D simulations come with a considerably higher computational cost, which has to be balanced with an increase in the physical relevance of the results of the 3-D model compared to other methods. In order to test the validity and relevance of the model in a quantitative way, in this paper we apply the MPM to a test case scenario of a relatively well-documented real avalanche event in Davos, Switzerland. The ‘Salezer’ snow avalanche event is described in Section 2. Details about how we simulate the avalanche are presented in Section 3. In Section 4, we show the results of the MPM simulation including comparisons with measurements of a drone survey after the event. In Sections 5 and 6, we discuss the simulation results and the comparisons with the measurements and draw conclusions about the applicability and future potential of the MPM model.

2. Avalanche field observations and data

The ‘Salezer’ avalanche occurred on 15 January 2019 in Davos, Switzerland, following a heavy snowfall that deposited 90 cm of snow in 3 d (SLF, 2019), resulting in a total snow depth of ~250 cm, in the avalanche release zone. To ensure the safety of the road and heliport below, the Salezer Horn slope is regularly triggered with explosive charges to cause controlled avalanche release. On 15 January 2019, a large avalanche was released on the ridge near the summit of Salezer Horn. Due to the large amount of erodible snow available along the path, the powder-snow avalanche reached a very large size, which over-passed the tunnel protecting the main road, crossed a car park and finally flew onto the ice surface of Lake Davos.

In the following sections, we describe the observation and measurement data from the avalanche event.

2.1 Release area and flow outline

A drone survey with a sensefly eBee RTK was carried out on 15 January 2019 after the event including photogrammetric measurement of the surface elevation and an orthophoto of the whole avalanche path. The mean flight altitude above ground was 195 m resulting in 679 images with a mean spatial resolution of 4 cm covering an area of 2.11 km² in total. Based on the orthophoto it was possible to identify three release areas and approximately map the outline of the dense flow of the avalanche. The avalanche control crew in the helicopter reported, that after the avalanche started in the primary release area, the flow of the avalanche itself led to the destabilization of the two secondary release areas. With high probability, the secondary releases were triggered sequentially by the disturbance induced in the snow cover, caused by the main body of the avalanche flowing by. The crown of the primary avalanche release was located at an elevation of 2456 m above sea level (a.s.l.) close to the summit of Salezer Horn, while the lowest point of the run-out was at 1556 m a.s.l. Thus, overall the avalanche covered a height difference of 900 m and a path length of ~2.5 km.

Furthermore, experts of the WSL Institute for Snow and Avalanche Research SLF extracted the approximate outline of the avalanche dense flow based on the snow surface texture visible in the orthophoto. The inset in Figure 1b visualizes that the distinction of the dense flow and the powder part of the avalanche was not always obvious from the orthophoto of the drone mapping and, therefore, can only be considered as approximative.

Figure 1. Panels a and b show an overview map and the orthophoto mapped from the drone survey of the avalanche track with release areas (red, orange and yellow shaded areas), and dense flow outline (purple), respectively. The inset in panel b shows a close-up of the granulation patterns in the dense part and the powder part of the avalanche, as well as the undisturbed snow cover. Map source: Swiss Federal Office of Topography.

2.2 Erosion and deposition

From the drone survey images, we calculated a digital surface model (DSM) with a spatial resolution of 10 cm. Due to the boundary conditions with the high avalanche danger and the start of the World Economic Forum with the corresponding closure of the airspace, we could not distribute ground control points. However, due to the eBee RTK capability the geolocation accuracy in the range of centimeters is possible. The intrinsically calculated values for the surface models and the orthophotos are x = 2.58 cm, y = 2.68 cm and z = 3.69 cm. These values agree with previous campaigns, where we achieved similar geolocation accuracies applying check points measured with differential GNSS. A second, snow-free flight was performed on 24 July 2019 (e.g. Eberhard and others, 2021).

We calculate the snow height distribution on the terrain after the event shown in Figure 2 by subtracting the snow-free DSMfrom the snow covered DSM. Although the snow height distribution before the event is not available in this case, the data provide good indication of where snow was eroded or deposited by the avalanche. However, in the absence of accurate snow height distribution data before the event and with a potentially considerable amount of deposits in the lake, it is not possible to make an accurate mass balance for this event.

Figure 2. Snow height distribution calculated from the photogrammetric drone survey. The inset shows a close-up of the primary release area marked with the red-dotted outline. Map source: Swiss Federal Office of Topography.

From the snow height distribution after the event, we can infer that the release height at the crown of the three avalanche release areas is highly variable, because of locally large deposits of wind-drifted snow. The release heights in all three release areas vary from ~0.5 up to 2.0 m. In Figure 2 the outline of the release areas is visible from the distinct drop in snow height, e.g. shown in the inset on the left for the primary release area.

Because the water level of the lake is reduced in winter, the deposition height is not accurate in the area of the lake. Moreover, snow which is cleared from the roads of Davos is deposited by the local authorities at the south-western tip of the lake. In Figure 2 this is visible from the dark red triangle in this region, which is therefore not relevant for our analysis.

2.3 Front velocity

At the south-western tip of the lake several persons were present during the event, recording a part of the avalanche with their mobile phone devices, while the avalanche was approaching. For the analysis, we use a private video, which is available online (Youtube, 2019). We extract the approach velocity of the avalanche front from the video by defining four control points along the flow path, which are shown in Figure 5b. Thereby we calculate the velocity from the elapsed time in the video and the distance between the control points. The control points are at the entrance (point 1 in Fig. 5b), near the middle (point 2) and at the exit (point 3) of the ‘Salezer Tobel’ gully, as well as at the edge of the avalanche tunnel roof protecting the main road (point 4). In the first two sections (points 1–3) the avalanche flows inside the gully with an average slope of 46$^\circ$. In the last section, between the exit of the gully and the tunnel roof (points 3 and 4), the avalanche flows on a wide open slope with an average inclination of 26$^\circ$. The average velocities between points 1 and 2, 2 and 3 and 3 and 4 are 42, 47 and 28 m s⁻¹, respectively. By combining extreme values of the ranges of the position and time span extracted from the video for the velocity calculation, we obtain an estimate of an uncertainty of up to ±5 m s⁻¹ of the approach velocity. The inaccuracy of the velocity estimate mainly arises due to the perspective view and the temporal resolution in the video.

2.4 Avalanche flow on the lake

In the run-out zone, the avalanche was interacting with Lake Davos. This bears the risk of generating an impulse wave, which could potentially endanger further infrastructure beyond the run-out of the avalanche. However, in the present case study this was not observed, as on the day of the event, the lake surface was 5 m below the maximum capacity and was covered with an ice sheet. The blasting crew in the helicopter reported that the ice at the side of the avalanche was only starting to crack ~10 s after the avalanche head stopped at the other side of the lake as shown in Figure 3. Considering that the terrain close to the impact point is almost flat, we assume that the avalanche flew almost parallel to the ice surface, and the normal forces exerted by the avalanche on the ice surface were low compared to a steeper impact. This makes it less likely for the ice to break and an impulse wave to be generated due to the impact.

Figure 3. Photographs of the avalanche flowing into the lake taken from the helicopter crew. The image in panel a is taken at the time when the avalanche reached the other side of the lake. The images in panels b and c are taken 12 and 30 s after the image in panel a, respectively. The blue arrows and dots mark the north direction and the location of the south-western tip of the lake, respectively. Photographs: V. Meier.

3. Numerical modeling of the event with MPM

In this study, we aim to test the possibilities and limitations of a novel fully 3-D numerical MPM model to simulate snow avalanches. In the numerical model, we distinguish two main components: the MPM solver and a constitutive material model described in Sections 3.1 and 3.2, respectively.

In Sections 3.3 and 3.4, we describe how we represent the snowpack on the day of the event and the original topography in the numerical model, respectively.

3.1 The material point method

The MPM solves the conservation of mass and momentum equations in a hybrid Lagrangian and Eulerian way. On the one hand, the Lagrangian particles (material points) are a discretized representation of the continuum material and advect material properties such as mass, velocity and momentum. On the other, an Eulerian background mesh is used to compute forces, solve the equation of motion and apply boundary conditions. To map the material properties between the Lagrangian particles and the Eulerian grid, transfer functions are used, which interpolate the material information from the particle positions to the grid nodes. In this study, we use an initial particle density of six particles per gridcell and the affine particle-in-cell method (Jiang and others, 2016, 2017) as a transfer scheme. In our scheme, we use quadratic B-splines as transfer function which have a span of 1.5 dx on both sides. Because in MPM the material is represented by particles moving in space with a non-deformable background mesh, this method allows us to simulate large material deformations, whereas in other methods large mesh distortion may lead to numerical instability.

For more in-depth information on the implementation of the numerical MPM scheme and the constitutive material laws, we encourage the reader to revisit the relevant publications, in which the solver and constitutive model were already extensively tested (e.g. Gaume and others, 2018; Li and others, 2021; Cicoira and others, 2022).

3.2 Constitutive material model for snow

To simulate a particular material with MPM, a constitutive model, which relates the deformation gradients to the stress state in the material, and the corresponding material properties are needed. In this study, we use the cohesive Cam Clay constitutive model to simulate snow (Gaume and others, 2018). This model has proven to perform well in simulating important features of the mechanical behavior of snow in avalanches, such as granulation, fracturing, hardening and softening. This capacity enables the model to capture e.g. levee formation, roll waves, erosion and deposition (e.g. Cicoira and others, 2022; Li and others, 2022b). An essential characteristic of our finite strain elastoplastic model is its capacity to encompass both the behavior of static snow cover distributed over the whole terrain for potential entrainment, and the flowing avalanche snow, which also originates from an unstable portion of the static snow cover itself. Hence, no arbitrary condition is needed to distinguish release and entrainment, but the entrainment process may occur naturally in the simulation. The cohesive Cam Clay model defines the material's yield surface as

(1)$$y( p,\; q) = ( 1 + 2\beta) q^2 + M^2( p + \beta{}p_0) ( p-p_0) $$

where p ₀ is the compressive strength, M the slope of critical state line and β the ratio of the tensile strength σ _ten and p ₀. In Eqn (1), p and q are the mean Kirchhoff stress and the von Mises equivalent Kirchhoff stress, respectively. They are defined as

(2)$$p = -{\rm tr}( {\boldsymbol \tau}) /d$$

(3)$$q = \sqrt{3/2{\, }{\rm dev}( {\boldsymbol \tau}) \, \colon \, {\rm dev}( {\boldsymbol \tau}) }$$

with the Kirchhoff stress tensor $\boldsymbol \tau$, and tr($\boldsymbol \tau$), dev($\boldsymbol \tau$) its trace and deviatoric part, respectively.

If the stress state exceeds the yield criterion in Eqn (1), the trial p–q-state outside the yield surface is projected back to the surface and a hardening law is used to adjust the yield surface. The hardening and softening of the material are calculated as follows:

(4)$$p_0 = K\, \sinh( \xi{}\, \max( -\epsilon{}^{\rm p}_{\rm v},\; 0) ) $$

In Eqn (4), K is the bulk modulus, ξ the hardening factor and $e^{\rm p}_{\rm v}$ the plastic volumetric strain. After the initial yielding of the static snow cover, the model allows us to describe the dynamic behavior of snow through a softening mechanism by changing the slope of the critical state line from the initial M to M _flow (Gaume and others, 2018).

It is important to note that in this implementation the interaction of the particles with ambient air is not captured. Hence, in our numerical model, we only reproduce the dense flow part of snow avalanche, where the physical effects of the ambient air and its interaction with the snow particles are negligible. The large powder cloud reported in the real avalanche is thus not considered here (Section 2).

3.3 Snow cover modeling

In our MPM simulations, we distribute snow cover all over the terrain along the avalanche path, mimicking an initially static snow cover as in reality. The avalanche flow is initiated by unstable sections of the snow cover, where the weight of the snow cover is not sufficiently counterbalanced by the friction forces at the ground and exceeds the yield criterion described in Section 3.2. Similar to entrainment in reality, also in our simulation, the static snow cover on the terrain can be entrained by the flowing snow, if the stress between the stationary and the flowing mass is high enough to exceed the yield limit.

In an attempt to model the snow conditions, including the snow mechanical properties and the erodible snow volume along the avalanche path, on the day of the event as close to reality as possible, we numerically simulate the layering and height of the snowpack with the SNOWPACK model (e.g. Lehning and others, 2002). We perform these simulations based on meteorological measurements for two locations near the avalanche track. The first station WFJ2 is located 1.3 km from the release (46.82945$^\circ$ N, 9.80909$^\circ$ E) at an elevation of 2540 m a.s.l., and thus, representative for the snow conditions in the release area. The second station SLF2 is located 250 m from the lake (46.81264$^\circ$ N, 9.84813$^\circ$ E) at the same elevation as the avalanche run-out at 1564 m a.s.l. Consistent with engineering guidelines (Margreth, 2007), we assume that the snow height increases linearly with altitude between these two stations. We visualize the simulated snow profiles at the two stations on the left and on the right sides in Figure 4.

Figure 4. Vertical snow profiles at stations WFJ2 (left) and SLF2 (right) with simplified snow layers (middle) interpolated linearly between the two stations.

Because we are limited to the spatial resolution of dx = 0.7 m (see Section 3.4) by our computational resources, in the MPM simulations we simplify the snow cover to only consist of two distinct layers. A lower layer (1) with older, consolidated and well solidified snow, corresponding to the red and blue colored grain types in the simulated profile. At the top of the avalanche track layer (1) has a thickness of 2 dx = 1.4 m and 1 dx = 0.7 m at the elevation of the run-out. The upper simulated snow layer (2) corresponds to the fresh snow deposited in the days just before the event, and corresponds to the grain types colored in light and dark green in Figure 4. The snow in layer (2) is fine-grained and less dense than the lower layers. At the top of the avalanche track, layer (2) has a thickness of 2 dx = 1.4 m, and 1 dx = 0.7 m at the elevation of the run-out. Hence, as a sum of layers (1) and (2), the simulated snowpack has a height of 2.8 m at the top of the avalanche track and 1.4 m in the run-out. In areas where the slope angle is larger than 50$^\circ$, we only deposit the lower layer (1) of snow on the terrain, as in reality snow cannot accumulate in considerable amounts in such steep terrain (McClung and Schaerer, 2006).

The mechanical properties of snow are notoriously difficult to assess, as the behavior depends on the complex crystalline micro-structure of snowpack which is constantly transformed by metamorphosis processes (e.g. Bader and others, 1939; Hagenmuller, 2014). Hence, the mechanical snow properties are highly sensitive to the atmospheric and load conditions, and may vary across multiple orders of magnitude as a consequence. For our simulation, we therefore use estimates of the mechanical properties of the old snow layer (1) and fresh snow layer (2), as summarized in Table 1. We estimate these mechanical parameters based on mechanical test measurement values from literature (e.g. Mellor, 1974; Jamieson and Johnston, 1990; Casassa and others, 1991; Shapiro and others, 1997; Willibald and others, 2020) and previous modeling work with MPM (e.g. Li and others, 2020, 2021, 2022b; Gaume and Puzrin, 2021; Wolper and others, 2021; Cicoira and others, 2022). While some of the parameters such as E, ν, ρ, σ _ten and M can be estimated based on a well-founded set of measurement data, others, such as the dynamic quantities M _flow and ξ, are harder to measure and therefore less measurements exist. We discuss the choice of these parameters in Section 5.

Table 1. Mechanical properties of the simplified snow layers (1) and (2)

^aβ in Eqn (1) is calculated as β = σ _ten/p ₀.

^bThe internal friction angle ϕ is calculated from M with: ϕ = asin(3M/(6 + M)).

Because layer (1) consists of old and well-consolidated snow, we implement a higher density of $\rho = 250\, {\rm kg}\, {\rm m}^{ {-}3}$, a compressive strength of p ₀ = 200 kPa and a tensile strength of σ _ten = 5 kPa, compared to the fresh snow in layer (2), for which we implement $\rho = 150\, {\rm kg}\, {\rm m}^{ {-}3}$, p ₀ = 180 kPa and σ _ten = 1 kPa (e.g. Jamieson, 1988; Jamieson and Johnston, 1990).

The rest of the parameters of the mechanical model in Eqns (1)–(4) are equal for both layers (1) and (2).

3.4 Modeling of the topography

For our case study, we simulate the avalanche flow on a terrain surface based on a DEM obtained from the Swiss Federal Office of Topography with a resolution of 2 m. For simulating the snowpack and the avalanche with MPM, we discretize the whole bounding volume of the avalanche track with a spatial resolution of dx = 0.7 m leading to a total number of 23 million particles. This is at the limit of what our current computational infrastructure (126 GiB Memory, 36× 3.00 GHz Intel^® Core^TM i9) is able to handle.

Because we explicitly simulate the erodible snow cover on the entire terrain, the avalanche front, a key determinant of avalanche dynamics, predominantly interacts with the erodible snow rather than the terrain contour, in contrast to many state-of-the-art numerical avalanche models. Consequently, in our simulations, terrain friction assumes a subordinate role in influencing flow resistance, but primarily serves to stabilize the erodible snow cover on the terrain, particularly in regions with steep slope angles.

Consequentially, as mentioned in Section 3.3, the avalanche flow is initiated where the load of the weight of the snow cover induces yielding of the material, because the weight is not sufficiently counterbalanced by friction at the ground. In order to stabilize the static snow cover on the steep terrain at an altitude above 1700 m a.s.l., we implement a ground friction coefficient f _c = 1.0. In the lower elevations, the terrain is flatter and a ground friction coefficient of f _c = 0.33 is sufficient to stabilize the snowpack.

With the spatial resolution of 0.7 m, we are not able to resolve the natural release process, which includes the collapse of a ~10 mm thick weak snow layer. Therefore, to destabilize the static snow cover in the primary release area, instead we implement a reduced ground friction coefficient f _c = 0.33 in the region of the primary and the two secondary release areas mapped in the drone survey. With this setup the avalanche flow is initiated in the region of the primary release due to the steepness of the terrain. In the secondary release areas the slope is slightly less steep and the snowpack is meta-stable. This means that the snow cover is initially stable and the snow only starts to flow due to the disturbance induced by the avalanche flowing nearby.

The two ground friction coefficients f _c = 1.0 and f _c = 0.33 used in our model are thus calibrated to capture the stability or instability of the snow cover in the real event. In this context it is also important to note that, due to the transfer functions described in Section 3.1 the boundary friction not only affects particles directly at the boundary but up to a distance of 1.5 times dx from the terrain contour away. We highlight and discuss the influence of the boundary friction on the simulation results in Sections 4 and 5.

Moreover, due to our limitations of computational power, we are not able to fully resolve the interaction of the avalanche with the lake and the ice in the run-out zone. While MPM is well suited to simulate multiple materials and their interaction in a single simulation without the need of specific coupling, the volumes of the ice and the water body in addition to the snowpack on the whole terrain make the simulations computationally too heavy to run on our current infrastructure. Hence, in agreement with the observations of the real avalanche described in Section 2.4, which show that the ice sheet does not break immediately when the avalanche crosses the lake, we assume that the avalanche head is gliding on an intact ice surface until it reaches the maximum run-out. We therefore simulate the ice surface of the lake as a solid with a reduced friction coefficient of f _c = 0.1, as we assume that the basal friction for the flow on the ice is low compared to the rest of the terrain.

4. MPM simulation results

4.1 Avalanche front approach velocity

In land-use planning, the avalanche velocity is important for practitioners to calculate the impact pressure and thus to define different hazard levels. Because avalanches are complex, 3-D and time-dependent flows, the velocities of different parts of the flowing snow within an avalanche may greatly vary even at a single instant (e.g. Sovilla and others, 2018). In the present analysis, we consider the avalanche front approach velocity v _front, representing the speed at which the avalanche front moves down-slope. Although v _front does not capture extreme local velocity peaks, this quantity is a good indicator of the dynamics at the avalanche front. In order to analyze if this crucial dynamic quantity is reproduced well in the numerical model, we compare the simulated avalanche approach velocity to the approximate approach velocity extracted from a video taken by an eyewitness (Section 2).

In order to extract the simulated front velocity shown in Figure 5a, we define the avalanche front by applying a particle velocity threshold of 1 m s⁻¹, which we use in all our analyses in the present article, to distinguish the static snowpack from the flowing avalanche mass. We define the avalanche front as the point of the flowing mass, which is furthest downstream the slope, and thus at the lowest elevation. The velocity is then calculated by dividing the distance covered by the avalanche front in a time interval Δt = 2 s by the time interval Δt. The fluctuations indicated by the error bars in Figure 5a indicate peak values of v _front if we choose Δt = 0.25 s, which is the maximum temporal resolution at which we export our simulation results. A sensitivity analysis on v _front and Δt is provided in the Supplementary material.

Figure 5. Panel a shows v _front extracted from the MPM simulation (solid blue line with fluctuations visualized by the error bars), as well as a comparison of the time-averaged simulated v _front (dashed blue line) compared to the approach velocity extracted from the eyewitness video (dashed red line) over the same time periods. The black-dashed lines and the corresponding numbers indicate the time at which the avalanche front passes the locations used to calculate the front velocity from the video. Panel b shows the same locations marked with crosses and video frames of the avalanche passing these locations in the insets. The main avalanche flow path is indicated with the red-dotted line. Map source: Swiss Federal Office of Topography.

In Figure 5a, we observe that during the first 50 s the avalanche approach velocity increases initially, with the exception of two main velocity drops after 20 and ~40 s in the simulation. These drops can be attributed to the release of secondary release areas. Indeed, our algorithm detects the accelerating particles in the release areas as the new front since they are further down the slope than the head of the avalanche itself. After the onset of the flow in the secondary releases, the front velocity in Figure 5a increases again as the mass builds up momentum.

Consistent with the steepness of the avalanche path, in Figure 5a we observe the highest avalanche approach velocities of ≈50 m s⁻¹ in the section of the gully between ~55 and 75 s. After the exit of the gully the avalanche flows on the flatter terrain between points 3 and 4, and v _front starts to decrease. The avalanche finally stops at 103 s on the other side of the lake.

In Figure 5a, the error bars indicate that the simulated v _front exhibits large fluctuations. It is important to note that the velocity peaks up to 100 m s⁻¹ are short lived and are most probably generated by a transient structure forming at the avalanche front, and therefore, do not necessarily correspond to the avalanche's approach speed. This peak velocity is the maximum of a velocity fluctuation and is representative of snow particles moving in transient flow structures, such as surges, which are faster than the avalanche approach velocity.

In order to make a direct comparison between the simulated v _front and the avalanche approach velocity extracted from the video, we average the simulated v _front (blue-dashed lines in Fig. 5a) in the same segments as in the video (red-dashed lines in Fig. 5a). We find a good agreement between the simulated and recorded average front velocities in all three segments. For the first segment (points 1 and 2), the simulated velocity is 5.1 m s⁻¹ higher than the one extracted from the video, which is the maximum absolute error and is almost within the error of 5 m s⁻¹ we estimate for the approach velocity extracted from the eyewitness video. The relative error between the velocity extracted from the eyewitness video and the simulated front velocity averaged over the same period is 12.1, 6.3 and 10.1% between points 1 and 2, 2 and 3 and 3 and 4 in Figure 5b, respectively.

4.2 Avalanche flow outline and velocity distribution

Figure 6a shows a comparison between the flow outline of the dense avalanche flow mapped from the orthophoto compared to the simulated flow outline. Identical to the analysis in the previous Section 4.1, we use a velocity threshold of 1 m s⁻¹ to distinguish the static snowpack from the flowing avalanche mass, and thus define the simulated flow outline. Figure 6b shows the distribution of the maximum avalanche flow velocity magnitude max(|v|). To be able to visualize the depth and time resolved velocity data on the 2-D map, we calculate the depth-averaged velocity magnitude of the particles within 2 m × 2 m northing and easting aligned cells for every time step and take the maximum over all simulation time frames.

Figure 6. Panel a shows the outline of the simulated flow (white area, delimited by black line) compared to the dense flow outline (purple line). The domain boundary of the simulated snow cover is marked with the gray-dashed line. Panel b shows the distribution of max(|v|). Map source: Swiss Federal Office of Topography.

Overall, there is a good match between the simulated and measured flow outline shown in Figure 6a. The simulated avalanche reproduces the correct run-out distance, with the avalanche coming to rest at the other shore of the lake, as well as minor details such as small side arms breaking away from the main flow path. The markers (1)–(4) in Figure 6a highlight a selection of points, where we find major differences between simulation and measurements or which we consider important to evaluate and discuss the capacity of MPM to capture relevant dynamical processes. The most significant difference is the lateral spreading of the avalanche in the run-out area between the gully and the lake (point 1 in Fig. 6a), where the flow is narrower in the simulation compared to the drone survey. We identify another difference close to the houses of the settlement ‘Meierhof’ (point 2 in Fig. 6a), where in the simulation a small area of snow releases, but remained stable in the real avalanche event. Further minor differences can be found at the entrance of the gully and at the starving arm of the avalanche close to the upper secondary release area (points 3 and 4 in Fig. 6a, respectively), where the simulated avalanche eroded less snow than the real one.

When reporting relatively small errors between the flow outlines from the dense flow avalanche simulation and from the drone mapping of a powder-snow avalanche, it is important to be aware that the distinction of the dense flow and the powder part is not always obvious from the orthophoto of the drone mapping as mentioned in Section 2.1.

The distribution of the simulated maximum avalanche velocity magnitude over all time frames max(|v|) in Figure 6b shows a similar trend as the avalanche front approach velocity v _front in Figure 5a. In the upper part of the path above the gully, the avalanche is building up momentum, which is however interrupted by the secondary releases. In the middle section of the flow, where avalanche flows in the gully, the velocity maximum is high, and also exhibits large fluctuations similar to the fluctuations indicated by the error bars in Figure 5a. Consistently Figures 5a and 6b also show a rapid deceleration of the avalanche on the flatter and open slope between the gully and the lake.

4.3 Snow erosion and deposition

To check how well the 3-D MPM model is able to reproduce erosion and deposition patterns of the real avalanche event, we compare the snow height distribution measured in the photogrammetric drone survey (Section 2) shown in Figure 7a to the simulated snow height distribution in Figure 7b.

Figure 7. Panels a and b show the measured and simulated snow deposition height distribution, respectively. Panel c shows a comparison of the measured (black solid line) and simulated deposition heights (scattered data points, colored according to the density), along the transect marked with the red line in panels a and b. Map source: Swiss Federal Office of Topography.

In both panels a and b in Figure 7, we can see that outside of the white avalanche flow outline, there is the general trend of increasing snow height at increasing elevation. If we compare panel a to panel b, we see that although this tendency is captured in our model setup, the snow cover is clearly idealized in the numerical model. In reality (panel a), the snow height distribution outside of the avalanche outlines is not homogeneous. Toward the summit of Salezer Horn the variability of the snow height increases and varies between 0.1 and 6.0 m over a distance of 30 m in extreme locations. In contrast, in the numerical model, we implement a homogeneous snow height of 2.8 m near the summit of Salezer Horn, corresponding approximately to the average of the snow heights reported from the drone measurements. The agreement between our simplified snow cover in the model is better toward the bottom of the slope, where the real snow is distributed more homogeneously.

Figure 7 shows that the simulated avalanche eroded nearly all of the snow cover in large parts of the avalanche track, which is in good agreement with the measurements. Moreover, the location and height of large snow deposits in the simulation mostly coincide between panels a and b. We identify the largest differences between measured and the simulated snow deposition heights in the run-out zone on the slope between the gully exit and the tunnel protecting the road (points 3 and 4 in Fig. 5). There, the simulated deposition heights reach a maximum of 8.5 m, and are therefore a factor 1.5–2 higher than in the drone measurement.

Figure 7c shows a comparison of the measured and simulated snow deposition in a 0.7 m = 1 dx wide transect along the main flow path of the avalanche, visualized by the red line in panels a and b. Especially where the deposits are high, we can clearly identify that the numerical model captures material densification, as the snow density increases from the top of the deposits toward the ground. In locations, where deposits are up to 4 m high, the simulated density in the deposits can reach up to 483 ${\rm kg}\, {\rm m}^{ {-}3}$ close to the ground on average, while the maximum implemented snow density of the initial snowpack is $250\, {\rm kg}\, {\rm m}^{ {-}3}$. In the first 500 m of the avalanche path, we can identify the primary and the upper secondary release areas, where simulated and measured snow height suddenly drops by several meters.

At the entry and the exit of the Salezer Tobel gully, located at 800–1300 m on the x-axis in Figure 7c, the measured deposition heights vary between 2 and 4 m, while in the middle section of the gully almost no deposits are present. The simulated snow deposits are in good agreement except in the middle section, where the numerical model computes depositions heights of 2–4 m. On the flatter slope between the gully exit (point 3 in Fig. 5) and the lake, the numerical model and the drone survey both show that most of the snow on the main avalanche track is eroded and almost no deposits are present, as shown in Figure 7c.

4.4 Intermittent and transient flow features

Our simulations allow us to also closely investigate complex and time-dependent dynamic flow features, which evolve naturally during the avalanche descent. Figure 8a shows the temporal evolution of the simulated avalanche flow velocity at a fixed location in the Salezer Tobel gully between points 2 3 in Figure 5b. Figure 8b shows a rendered 3-D view of the flow shown in Figure 8a at the t = 63 s in the simulation when the avalanche front is at the location corresponding to the velocity data in panel a.

Figure 8. Analysis of the simulated avalanche front flow behavior in a fixed location. Panel a shows the temporal evolution of flow velocity near point 2 in Figure 5b as a function of the flow height. The inset shows a close-up of the same data at the flow front. Panel b shows a rendering of the avalanche front at the location where we extract the velocity in panel a. Panels c and d show the vertical velocity and density profiles at t ₁, t ₂, t ₃ indicated in panel a, respectively.

Each pixel in Figure 8a is colored according to the averaged particle velocity in cells of 2 dx × 1 dx × 0.5 dx, in the main flow direction, the transverse and the vertical direction, respectively. The flow velocity is highest at the free surface of the flow near the avalanche head with a maximum velocity of 57 m s⁻¹. The close-up of the flow front in the inset shows how the static snow cover colored in gray is entrained by the avalanche. The entrainment is also visible in panel c, showing the pixel velocity as a vertical profile, where at t ₁ only particles on top of the static snow cover and at t ₂ also the particles closest to the ground are moving. For the time step t ₁, panel d shows that the snow which is only just entrained by the avalanche remains relatively loose with densities smaller than 100 ${\rm kg}\, {\rm m}^{ {-}3}$ in the flowing part. At instant t ₂, the density is almost constant with some fluctuations around a mean value of 174 ${\rm kg}\, {\rm m}^{ {-}3}$, which is in between the initial densities of 150 and 250 ${\rm kg}\, {\rm m}^{ {-}3}$ of the two snow layers. Considerable compaction only occurs later, between t ₂ and t ₃, where the snow density increases up to a maximum of 475 ${\rm kg}\, {\rm m}^{ {-}3}$ at the bottom of the dense flow. As indicated by the density in the final snow deposits in Figure 7c, later the snow is not further compacted.

Intermittent flow structures in the frontal region of the avalanche, similar to the ones shown in Figure 8a, where a part of the snow mass is detached from the ground contour and the dense flow, are also observed from real scale experimental measurements of large powder-snow avalanches (Sovilla and others, 2018) flowing in a similar configuration in a gully in the Vallée de la Sionne (VdlS) full-scale test site in Switzerland (Ammann, 1999). For better visualization of these intermittent flow structures in the frontal region of the avalanche, Figure 8b shows a 3-D spatial rendering of the flow shown in Figure 8a, at the moment when the avalanche front is at the location where the velocity is analyzed in Figure 8a.

Figure 9a shows a qualitative comparison between the temporal evolution of the simulated vertical slope-normal component of the flow velocity for the same location used in Figure 8a, and measurements performed at VdlS at the front of a powder-snow avalanche using an upward-facing Frequency-Modulated Continuous-Wave (FMCW) radar measurement (e.g. Gubler and Hiller, 1984), which are displayed in the inset of the same figure. Two striking similarities are evident when comparing the simulated and measured flow features in panel a and the inset in Figure 9. First, we observe that the surface of the dense flow exhibits an undulated shape in both plots. Second, the comparison also shows that in the simulation, as well as in the FMCW radar measurements, large snow clusters are detached at a distance above the basal dense flow. The simulated slope-normal velocities in Figure 9a are overall small, up to ${\sim }10\%$ compared to the velocity magnitude (Figs 6b and 8a). Positive and negative velocities in Figure 9a, indicate that clusters of snow are moving upward and downward, respectively.

Figure 9. Panel a shows the simulated time evolution of the flow height (y-axis) and slope-normal velocity (color map) corresponding to the same location as in Figure 8a. The inset shows the temporal evolution of flow depth measurements from an upward-looking FMCW radar, installed in the gully of VdlS. Panel b shows the slope-normal velocity at t = 63.0 s and t = 92.5 s in a 500 m long transect in the gully and the terrain slope, in the top and bottom plots, respectively. The gray-dashed lines highlight the correlation of exemplary peak values in both plots. Panels c and d show the distribution of the slope-normal velocity v _n at time t = 63.0 s and t = 92.5 s. The transect for which the slope-normal velocity and the slope angle are visualized in panel b is a straight line between the two red arrow tips. Map source: Swiss Federal Office of Topography.

To better understand the relevance of velocity component in the flow-depth direction, in Figures 9b–d, we visualize the simulated spatial distribution of the slope-normal flow velocity v _n at the t = 63 s when the avalanche front is at the location for which the velocity data in Figures 8a, b and 9a is plotted, as well as at t = 92.5 s. Similar to Figure 5b, we calculate the depth-averaged slope-normal velocity of the particles within 2 m × 2 m northing and easting aligned cells. Figure 9b shows the variation of the terrain slope and the slope-normal velocity at t = 63.0 s and t = 92.5 s in the simulation in a 500 m long transect in the gully, which is indicated between the tips of the red arrows in panels c and d. The gray-dashed lines highlight the correlation between the peaks of slope-normal velocity and sharp changes of terrain slope in the top and bottom plots, respectively. While v _n is mostly smaller than 5 m s⁻¹, both curves of slope-normal velocity exhibit peaks of slope-normal velocity in the range of 5–10 m s⁻¹. A comparison of v _n at t = 63.0 s and t = 92.5 s in Figure 9b reveals that the peaks, particularly in the distance range of 120–170 m, tend to be higher for the green curve at t = 63.0 s. The green curve corresponds to a phase when the avalanche front is traversing the terrain at a higher absolute velocity (Fig. 8a), as opposed to t = 92.5 s when the same terrain section is being traversed by the tail of the avalanche at a lower velocity, and thus with lower kinetic energy.

In Figures 9c and d, it is evident that both upward and downward particle movements occur throughout the entire avalanche flow. Notably, elevated values of v _n are predominantly observed at the avalanche front. Moreover, high absolute values of slope-normal velocity v _n are most pronounced in the steep, channeled flow section in the gully, where the avalanche reaches the maximum velocities. Meanwhile, both panels c and d show concurrently that the magnitude of the slope-normal velocity is relatively small in the upper part of the avalanche path before the gully, where the avalanche is accelerating. Finally, to help the interested reader to gain a better insight in these complex and temporally and spatially highly variable flow structures, we include a rendered video of the simulated avalanche in the Supplementary material. In this video, we visualize the slope-normal velocity component of the flow.

5. Discussion

5.1 Model novelty and physical relevance

While previous studies tackled simulating 3-D depth-resolved avalanches on a full-scale real topography (e.g. Sampl and Granig, 2009), to the best of our knowledge, in this article we present the first simulation, where we additionally explicitly simulate the snow entrainment. The simulation domain is ~2.5 km long and 800 m wide, which results in a total volume of the simulation domain of 570 M m³ and 23 M simulated snow particles. We also simulate the snow conditions on the day of the event using measurement-driven SNOWPACK simulations and use corresponding estimates of the mechanical snow properties from literature.

Despite the aim to simulate physical processes as close to reality as possible, we have to simplify the simulations to keep the calculation time within reasonable limits. Due to limited computational power (see also Section 5.3), we do not explicitly resolve the ice sheet and the lake water in the simulation, but we consider the ice sheet as a rigid boundary. Because observations from the helicopter crew presented in Section 2.4 indicate that the ice only cracked with some delay after the avalanche head already reached the maximum run-out, we think this approximation is acceptable and should not influence the simulated run-out considerably.

Furthermore, due to the coarse-grid resolution dx = 0.7 m, we simplified the snow cover stratification in only two layers and assume that the collapsed weak layer where the avalanche releases was near the ground. This solution was acceptable in our case, but the course definition of the snow cover can become a problem, if in another event the snowpack is composed of thin layers with markedly differing mechanical properties or the weak layer is further from the ground.

Another simplification in the numerical model is the elevation-dependent snow height distribution according to engineering guidelines (Margreth, 2007). The comparison shown in Figure 7 with the photogrammetric drone survey indicates that the real snow distribution is characterized by a large variability in snow height, as a result of both, wind-induced preferential deposition of snow (e.g. Dadic and others, 2010), and previous avalanche activity. This may influence the avalanche flow, because in locations with large wind drift deposits the snowpack may easily become unstable, while in locations, where the snow is blown off by the wind or transported away due to previous avalanche activity, an avalanche may starve or not release. In order to improve this, the model could for example be coupled with an algorithm calculating the snow drift based on meteorological data in the specific topography (e.g. Alpine3D; Lehning and others, 2006), which would, however, significantly increase the model's complexity.

Despite these simplifications and rough estimates of the mechanical snow properties based on literature data, the simulation results are overall in good agreement with the real avalanche event, as shown in Figures 5–7, suggesting that even with these assumptions the 3-D MPM model is able to capture the most important flow processes in our case study. The good agreement between the simulated and mapped flow outline despite the major simplifications made to the snow cover definition also indicates that the detailed layering simulated with SNOWPACK only plays a minor role in the overall dynamic behavior. Based on a sensitivity analysis with altered snowpack characteristics (see the Supplementary material), we assume that the most important factor is the presence of an erodible snow cover of sufficient height and erodibility, i.e. low compressive strength. The interaction of the flowing avalanche with the static snow cover on the terrain allows for volume gain or loss of the avalanche by eroding or depositing snow on the path, which also governs the overall dynamics of the avalanche (Schweizer and others, 2009).

5.2 Insight into avalanche flow processes

Thanks to the 3-D nature of our MPM simulations, the model explicitly resolves snow erosion and deposition processes without the need for a conceptual model or empirical relationship. For example, Figure 8a shows how the snow cover is entrained by the avalanche front. Figure 7 shows that in our case study the model is able to reproduce the most important deposition patterns of the real event qualitatively. The simulated snow deposits are mostly located on terrain with moderate slope angles below $30^\circ$ below the gully, in agreement with the findings of Sovilla and others (2010), who state that snow deposition mainly occurs on terrain with slope angles $\leq \!33^\circ$. However, in the steep middle section of the gully, the model also simulates large snow depositions, which are not observed in the real event and are not likely to occur anywhere else in such steep terrain (Sovilla and others, 2010). A probable explanation for the large simulated deposits in the gully is the boundary friction in the model, which also acts at a distance as far as 1.5 dx due to the transfer functions used in the numerical scheme. Hence, where the terrain is concave and curvature is high enough, such that the 1.5 dx distance bands from both sides of the gully overlap, the boundary friction is applied twice to the particles in the overlapping zone. Similarly, the difference in lateral spreading in the run-out zone of the simulated avalanche compared to the lateral extent mapped from the drone data can partly be attributed to the influence of the boundary condition. Because the snow height decreases with elevation, the boundary condition, acting on particles at the same distance from the terrain ≤1.5 dx everywhere, influences a larger fraction of the snowpack in the run-out zone compared to higher elevations, where the snow cover is thicker. Finally, not only the simulation, but also the drone measurement may be fraught with error due to inaccuracies including e.g. the presence of high grass or bushes in the summer DSM, from which the snow surface height registered by the drone is subtracted (e.g. Vander Jagt and others, 2015). This may lead to a small underestimation of the measured snow deposition height, which is, however, considerably smaller than the difference in deposition height we observe between Figures 7a and b.

Figures 7c and 8d show that the model captures snow densification in a realistic way. Indeed, the range of density values with the highest densities near the bottom and the densification occurring progressively during the avalanche flow is consistent with field observations (e.g. Sovilla and others, 2006; Gauer and others, 2007). In the cohesive Cam Clay constitutive model, which we used for snow, the densification mainly depends on the hardening factor ξ in Eqn (4). In the present case study, we choose ξ = 0.1 according to Cicoira and others (2022), which results in density values of the simulated snow depositions close to values measured from real avalanche deposits (e.g. Sovilla and others, 2006; Gauer and others, 2007; Steinkogler and others, 2014; Issler and others, 2020).

As shown in Figure 5a, the averaged approach velocity extracted from the simulations matches with the front velocity extracted from the video. Figures 8a and 9 suggest that short-lived velocity peaks, akin to those in Figure 5a, could be generated by transient processes. These may include material jets expelled from the basal dense layer or pulsating activity at the surface of the basal dense layer induced by waves or surges. Indeed, such intermittent activity has also been observed in the frontal region of powder-snow avalanches at the VdlS (Köhler and others, 2018; Sovilla and others, 2018) and at the Ryggfonn full-scale test site in Norway (Gauer and others, 2007). While the origin of these transient structures in the measurements is not yet fully clarified, it is often assumed that turbulence in the suspension layer may play an important role for their origin and dynamics.

Although in our simulation we do not include the interaction with the ambient air, we still observe material clusters detached from the dense flow similar to full-scale powder-snow avalanches shown in Figure 9a. Moreover, panels b–d in Figure 9 show non-negligible slope-normal velocity components up to ~10 m s⁻¹. The plots of the slope angle and v _n for two different instants in a transect in the gully in Figure 9b imply that peak values of positive and negative slope-normal velocity are attained if the variations in slope angle and the kinetic energy of the flow are high. Consequently, our simulation results imply that a significant portion of the snow clusters observed in intermittent structures within powder-snow avalanches probably originates from the ejection of particles from the basal dense flow. The ejection takes place due to the interaction of the snow mass, flowing with high kinetic energy, and the terrain, characterized by large slope variations that redirect the momentum of specific portions of the flowing mass in the slope-normal direction.

Although a more comprehensive analysis of the simulated flow features could be conducted with improved field measurement data from the event, our case study already demonstrates the potential of 3-D MPM as a valuable tool to enhance the comprehension of the processes contributing to the particle-lading of the powder cloud. These processes may influence the particle concentration and frequency of the intermittent structures at the avalanche front, where the largest part of the flow energy and destructiveness are concentrated, and are, therefore, important for engineers to identify critical pressure peaks avalanches exert on infrastructure (Eglit and others, 2007; Mast and others, 2014; Brosch and others, 2021; Gorynina and Bartelt, 2023). In addition, equally important velocity profiles including the slope-normal component (Fig. 9), which is relevant for uprooting structures, can be extracted 3-D MPM simulations. Moreover, the level of physical detail in our results of this case study highlights that physics-based 3-D MPM have the potential to be used in research to increase the understanding in avalanche flow dynamics.

5.3 Current limitations and future developments

To date, and even with a computationally efficient method such as the 3-D MPM model we use in this study, fully 3-D simulations of real-scale events with vast extents such as in our case study are still challenging. In addition to the simplifications mentioned previously, e.g. the assumption that the ice sheet on the lake is rigid, we use the maximum grid resolution of dx = 0.7 m achievable with our computational resources. As mentioned earlier, this implies simplifications in the representation of the snow cover layering and the avalanche release mechanism (Section 5.1).

Furthermore, the coarse-grid resolution also influences the dynamic behavior, because the boundary condition affects particles up to 1.5 dx = 1.05 m away from the terrain surface due to the transfer functions used in the numerical scheme. In order to stabilize the initial snowpack on the terrain and create meta-stable snow cover conditions for the secondary releases, we implement calibrated friction values of f _c = 1.0 and f _c = 0.33 at the boundary. In the future, a more prediction-oriented numerical model could involve an implementation where the boundary condition exclusively influences adjacent particles. Additionally, incorporating pre-computed boundary friction values based on a hysteretic friction model (e.g. Daerr and Douady, 1999; Pouliquen, 1999) could facilitate the representation of meta-stable snow cover conditions. Once computed, these friction values may be applied in potential release areas identified through appropriate techniques, such as those outlined by Bühler and others (2018).

To compensate for the relatively high boundary friction values f _c = 1.0 and f _c = 0.33 needed to stabilize the snowpack, we choose M _flow as low as physically reasonable. Hence, we choose M _flow = 0.37 such that it corresponds to an internal friction angle of $\phi _{\rm flow} = 10^\circ$, which is in the lower range suggested by Casassa and others (1991) for very cold and dry snow. The static values M = 0.98 and $\phi = 25^\circ$ are in a normal range compared to measurements (Casassa and others, 1991; Platzer and others, 2007; Willibald and others, 2020).

Another major shortcoming of our model is that the interaction between the snow particles and the ambient air is not captured. This implies that the 3-D MPM model inherently simulates dense flow avalanches. Hence, air turbulence or fluidization, which debatably may occur due to pore pressure increase near avalanche head, are not taken into account. However, these processes have a considerable influence on the erodibility of the snowpack (Louge and others, 2011; Issler, 2022). The reduced snow particle mobility due to lacking fluidization is probably also an important reason, why in the run-out zone between the gully and the lake (points 3 and 4 in Fig. 5), we observe less lateral spreading in the simulations compared to the mapped outline from the drone survey. Probably the non-fluidized particles in the simulation are less mobile than in reality, leading to a channeling of the flow instead of lateral spreading. The smaller lateral spreading further results in an overestimation of simulated snow deposition heights, as the avalanche mass is distributed over a smaller area than in reality, and thus, leaving higher deposits. In the future the issue of the boundary condition influencing snow particles up to a distance of 1.5 dx from the terrain could be avoided by implementing an algorithm similar to BFEMP at the boundary (Li and others, 2022a).

Furthermore, we address the challenge of the computational cost of our challenges already now and in the future, e.g. by developing an ‘activation’ based simulation strategy, where the relevant equations are only solved for particles currently involved in physical action, instead of the whole static erodible snowpack. Moreover, in the future the MPM model should support highly parallelized GPU-based simulations in addition to CPU. A recent study showed that GPU implementation of MPM could make simulations, currently lasting up ~53 h, up to 16 times faster than CPU-based codes (Gao and others, 2018).

6. Conclusions

In this article, we tested the potential and challenges to simulate large, full-scale snow avalanches with a novel depth-resolved and fully 3-D MPM model. To get an indication of how well the model performs, we compare the simulation results to the well-documented Salezer snow avalanche, which occurred in January 2019 in Davos, Switzerland. Despite these simplifications, we find that the simulation results are in good agreement with the observations from the real avalanche, particularly the avalanche approach velocity extracted from an eyewitness video and the flow outline mapped in a drone survey. Furthermore, the model reproduces the most important erosion and deposition patterns of the real event in a qualitative manner. However, quantitatively the simulated snow deposits are, locally, up to twice as high as the deposits mapped in the photogrammetric drone survey. We identify two reasons for this discrepancy, which also highlight the two main limitations of the model. First, as the model does not include the ambient air in the simulation, we are limited to simulate the basal dense flow of the avalanche. Second, due to the particle-to-grid transfer functions we use, the boundary friction affects snow particles up to 1.5 dx away from the terrain, which introduces an artificially high resistance to the flow, which could be avoided in the future by implementing another transfer algorithm at the boundary (e.g. BFEMP, Li and others, 2022a).

Furthermore, we explore the potential of the 3-D depth-resolved model by analyzing intermittent flow structures near the flow front in the gully and find that numerous snow particle clusters are ejected from the dense basal layer and remain at a distance above the basal dense flow for few seconds, even though the turbulent interaction of the snow particles with the air is not simulated. We speculate that these flow structures are generated at sudden changes in the topography in the gully, where the avalanche flow passes with high kinetic energy.

Considering the level of physical detail of the results, especially concerning the transient flow structures at the avalanche front, we conclude that the model has a high potential to be used to perform in-depth analyses, particularly to identify critical impact pressure peaks due to transient flow structures. Furthermore, the model could also be used in research to investigate on dynamic flow features, which are difficult to measure in the field. Finally, in the context of a warming climate with changing frequency and characteristic of the snow avalanche hazard, physics-based modeling approaches such as the MPM model presented here, will become increasingly important for hazard assessment, as models calibrated with historic data, while valuable, may have limitations in capturing the evolving physical processes.