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, Reference Eberhard2021).

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.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, Reference Jiang, Schroeder, Teran, Stomakhin and Selle2016, Reference Jiang, Schroeder and Teran2017) 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, Reference Gaume, Gast, Teran, Van Herwijnen and Jiang2018; Li and others, Reference Li, Sovilla, Jiang and Gaume2021; Cicoira and others, Reference Cicoira, Blatny, Li, Trottet and Gaume2022).

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, Reference Gaume, Gast, Teran, Van Herwijnen and Jiang2018). 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, Reference Cicoira, Blatny, Li, Trottet and Gaume2022; Li and others, Reference Li, Sovilla, Ligneau, Jiang and Gaume2022b). 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, Reference Gaume, Gast, Teran, Van Herwijnen and Jiang2018).

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, Reference Lehning, Bartelt, Brown and Fierz2002). 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, Reference Margreth2007), 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, Reference McClung and Schaerer2006).

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, Reference Bader1939; Hagenmuller, Reference Hagenmuller2014). 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, Reference Mellor1974; Jamieson and Johnston, Reference Jamieson and Johnston1990; Casassa and others, Reference Casassa, Narita and Maeno1991; Shapiro and others, Reference Shapiro, Johnson, Sturm and Blaisdell1997; Willibald and others, Reference Willibald, Löwe, Theile, Dual and Schneebeli2020) and previous modeling work with MPM (e.g. Li and others, Reference Li, Sovilla, Jiang and Gaume2020, Reference Li, Sovilla, Jiang and Gaume2021, Reference Li, Sovilla, Ligneau, Jiang and Gaume2022b; Gaume and Puzrin, Reference Gaume and Puzrin2021; Wolper and others, Reference Wolper2021; Cicoira and others, Reference Cicoira, Blatny, Li, Trottet and Gaume2022). 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, Reference Jamieson1988; Jamieson and Johnston, Reference Jamieson and Johnston1990).

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.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, Reference Sovilla, McElwaine and Köhler2018). 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, Reference Sovilla, McElwaine and Köhler2018) flowing in a similar configuration in a gully in the Vallée de la Sionne (VdlS) full-scale test site in Switzerland (Ammann, Reference Ammann1999). 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, Reference Gubler and Hiller1984), 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.