Simulation of cold-powder snow avalanches considering daily snowpack and weather situations

Snow avalanches are rapid gravitational mass movements that pose a significant hazard to both humans and infrastructure, including traffic lines. Risk management in mountainous regions usually relies on the experience of avalanche experts, observations in the field, weather and snowpack measurements, and numerical simulations.

Ensuring road safety requires considering daily weather conditions, snowpack characteristics and terrain features. To include a numerical model in the decision-making process for road safety, it is essential to incorporate all these factors and utilize in situ measurements as input parameters for the simulations.

This study investigates the predictive capabilities of the numerical simulation model RAMMS::EXTENDED, an extended version of the well-established RAMMS (Rapid Mass Movement Simulation) software, to estimate avalanche runout distances along an important infrastructure corridor in the Dischma valley near Davos, Switzerland. Specifically tailored to cold-powder avalanche dynamics and taking into account the temperature of the snowpack and entrainment, our inquiry utilizes meteorological station measurements as an input to evaluate the model's performance.

In this paper, we begin by providing an overview of the model, examining its physical and practical aspects. We then conduct a sensitivity analysis on input and system parameters, focusing on avalanche dynamics representation. Leveraging drone-based observational data, we perform a comparative analysis to validate the simulation results.

In addition to recalculating avalanches due to the sensitivity analysis, we show that we achieve meaningful predictions of the avalanche runout distance for cold-powder avalanches by incorporating snow height and snow temperature measured by weather stations at two different altitudes near the avalanche release zone. In the future, a refined version of this approach could allow for near real-time hazard assessments, which has the potential to significantly improve the decision-making protocol for road closures and reopenings. Furthermore, we plan to calibrate the model for wet-snow avalanches to cover a larger range of weather and snowpack scenarios.

1 Introduction

In some avalanche-endangered regions, a robust network of measurement stations provides point information on snow depth, temperature, and wind speed. In this publication, we aim to explore how these data sources can be leveraged to simulate avalanche runout distances to represent daily conditions. Using the example of road safety, we will demonstrate how this information can be applied to protect infrastructure. The same approach is adaptable for other safety frameworks, such as ski slopes or buildings in avalanche-prone areas.

Due to economic and environmental constraints, many mountain roads cannot be effectively protected using long-term technical measures to prevent avalanche release (avalanche defence structures) or inundation (road alignment, snow sheds and tunnels). Therefore, local hazard experts must make decisions to close roads and stop all traffic during avalanching periods. These decisions rely on information from the warning services, the interpretation of measurement data and experience (Stoffel and Schweizer, 2008). Increasingly computer-based expert systems such as the nearest-neighbour model for regional avalanche forecasting called NXD (Brabec and Meister, 2001) or AI systems are applied to help the hazard experts.

In this work, we focus on adding information to the decision-making process by combining data from weather stations, numerical modelling and drone measurements. The primary goal is to answer the question of whether an avalanche could reach a road under specific snowpack and weather conditions. Having this information at hand could enhance road safety mitigation and reduce the road closure times to a minimum (Sheng et al., 2022). For this approach to be successful and to include all possible avalanche paths along a road, accurate reports of snowpack and weather conditions are needed. These data must be collected as close as possible to the specific avalanche path.

The utilization of numerical avalanche dynamics modelling to enhance road safety signifies an important paradigm shift in avalanche engineering. While numerical models have traditionally and extensively been adopted for generating hazard maps and designing avalanche defence structures along specific avalanche paths, they often do not include crucial snow properties such as snow-cover layering, density, temperature or moisture content to represent daily conditions. Avalanche fracture heights are typically determined through statistical analysis of long-term snow accumulation data from measurement stations (Salm et al., 1990). Following an approach pioneered by Voellmy (Voellmy, 1955), extreme avalanche events are typically addressed using calibrated parameters derived from historical avalanche occurrences (Gruber and Bartelt, 2007). While this approach is suitable for hazard mapping, it fails to leverage recent advancements in automatic weather stations or drone measurements (Bühler et al., 2017). Consequently, the output from numerical models available to local hazard engineers for deciding whether to close a road is limited, as it becomes challenging to correlate specific measured data with potential avalanche runout distances.

A first approach to use numerical simulations for road safety is presented in Keylock et al. (1999), where typical avalanche scenarios are precalculated. For operational daily predictions, a system is implemented in Chile as described in Vera Valero et al. (2018). It uses avalanche dynamics modelling based on RAMMS to predict whether an avalanche reaches a road. For the input data, it relies additionally on the simulation tool SNOWPACK (Cerda et al., 2016; Lehning et al., 1999).

As we experiment with new applications and develop more complicated modelling chains, we also place new demands on these existing numerical models. To accurately represent snow and weather conditions, a model must have the capability to encompass avalanches with different flow regimes (including wet-, mixed- and dry-snow avalanches) and consider snow-cover entrainment and mass growth, the braking effects of different forest compositions, and (most importantly) the influence of snow temperature. Existing avalanche dynamics models which focus on the flowing regime, ELBA (Keiler et al., 2006), OpenFOAM (Rauter et al., 2018), SAMOS-AT (Sampl and Granig, 2009), AvaFrame (ava, 2023) and RAMMS::Avalanche (Christen et al., 2010), only meet parts of these requirements.

In this publication, we focus on the investigation of a system for cold avalanches. We utilize three well-documented avalanches that overflowed a mountain road near Davos (Switzerland) to investigate how avalanche dynamics models can effectively be used with weather station data. The avalanches were artificially released and developed into a mixed flowing-powder type. Post-event drone scans provided detailed information on runout and snow-cover distribution. We apply an extended RAMMS model that includes snow temperature (Vera Valero et al., 2016), entrainment (Bartelt et al., 2018), and formation and propagation of the powder cloud (Zhuang et al., 2023a). The model was calibrated using avalanches observed at Vallée de la Sionne (VdlS) (Ammann, 1999), considering only those that did not reach the counterslope according to P. Bartelt (Perry Bartelt, personal communication, 26 January 2024). Additionally, avalanches from winter 1999 in Switzerland were used, which are presented in Vallet et al. (2001). In the first part of this paper, we provide an overview of the current version of the model, summarize the recent literature on RAMMS::EXTENDED and present the current calibrations. In the second part, we back-calculate the observed avalanches using temperature data from nearby snow monitoring stations and the parameter set from VdlS. We show how the model reacts to changing boundary conditions and the sensitivity of model performance to variation of parameters. Our results highlight the challenges of using avalanche dynamics models for road safety applications.

2 Observations and methods

2.1 Avalanche events, Davos, 15 January 2019

We examine three separate avalanche events that took place in the vicinity of Davos, Switzerland. These incidents occurred in mid-January 2019 in the Dischma valley during a cold-weather period, leading to the formation of mixed flowing-powder avalanches (see Fig. 1). The three avalanche tracks are located on the northeastern slope of Brämabühl (Davos) and have the names Wildi, Rüchi and Chaiseren. The release zones are all located at roughly 2300 m a.s.l. on a northeastern aspect. The tracks drop between 650 and 750 m in elevation, running over a well-used crosscountry skiing track and a local road. The tracks are somewhat channelized below the release zone, but they open to wide, laterally unconstrained runout zones at the valley bottom. In the past, avalanches from these tracks have blocked the road connecting the inhabitants of the valley and Davos, wooden buildings have been destroyed and trees in the surrounding forests have been blown over by avalanche air blasts. Hazard maps exist for all three tracks.

Figure 1 On the left side is the depiction of the measured release zones of the three avalanche paths named (1) Wildi, (2) Rüchi and (3) Chaiseren. The measured outline of the three avalanches at Brämabühl on 15 January 2019 from the drone data are marked in red (map source: Federal Office of Topography). On the right side is the deposition of the three powder avalanches in the Dischma valley (Davos, Switzerland) originated from northeast-facing slopes. Powder avalanches often reach the valley road after traversing a flat runout zone (photos made by Vali Meier, SOS Davos Klosters). The grid lines show a distance of 1 km.

The 2019 avalanche events are unique since substantial snowfall preceded the cold spell, with snow depths as high as 2.5 m measured at nearby snow-cover monitoring stations. A strong winter storm passed through the valley on the 14 January 2019 with strong winds that redistributed snow on the slope (Glaus et al., 2024). The avalanches eroded a deep, cold snowpack, which contributed to the formation of powder avalanches. The avalanches were artificially released, and additional photos from the helicopter during the avalanche control provided estimates of the powder cloud speed (around 30 m s⁻¹ in the runout zone) and height (approximately 40 m). For the Rüchi and Chaiseren paths, the approximate cloud speed can be estimated by analysing the position of the powder cloud front over time as illustrated in Fig. 2. The avalanches considerably increased in mass after release due to snow entrainment. In the following days, a field campaign was carried out to gather data concerning the location of the release zones, entrainment heights and avalanche runout lengths. The data were coupled with temperature data from nearby snow monitoring stations (Weissfluhjoch and IMIS SLF2), which are situated 6 and 3 km away. The Weissfluhjoch station (2536 m a.s.l.) is located roughly at the same elevation as the avalanche release zones, and the IMIS SLF2 station (1570 m a.s.l.) is at the same elevation as the runout zones (1600 m a.s.l.).

Figure 2 Estimation of the powder cloud velocity from the helicopter-captured images for the Rüchi path. The velocity was calculated by comparing the distance travelled by the cloud to the time interval between the images (photos made by Vali Meier, SOS Davos Klosters).

On-site data collection was conducted with structure-from-motion photogrammetry (Bühler et al., 2011) using drones, allowing for the measurement of snow heights by comparing post-avalanche elevation surfaces with bare ground surfaces during the summer obtained from the Swiss Federal Office of Topography (SwissTopo) (Swisstopo, 2024). In Fig. 3, the measured post-avalanche snow heights for all three avalanche tracks are shown. These drone data allowed for the delineation of the release zones, thanks to a clearly visible stauchwall. Avalanche fracture height (d₀) perpendicularly to the terrain could be estimated by comparing measured snow heights at similar altitudes to the average snow height in the release zone after avalanche release. Additionally, from the drone measurements, the snow distribution gradient ∇D could be determined. These values are reported as the average decrease in snow-cover height per 100 m drop in elevation (Fig. 6).

Figure 3 Snow-depth distribution mapped photogrammetrically with the eBee RTK drone on 16 January 2019 for the entire area (left) and zoomed in to the avalanche release zone (right) after the avalanches were triggered. Significant wind redistribution effects are visible in the image (map source: Federal Office of Topography).

Figure 4 Visualization of the idealized plane inspired by the steepness of the Rüchi avalanche path. (a) Bird's-eye view of the plane with contour lines and altitude marked. The colour gradient indicates slope angles, ranging from low (flat areas) to high (steep areas). The scale for the runout distance shows the avalanche length for a release zone which is placed at 2200 m a.s.l. (b) Cross-section of the plane, illustrating the elevation profile along a selected transect.

To better understand temperature gradients ∇T along the path (Fig. 6), we analysed snow pits concurrently with temperature readings at Weissfluhjoch and SLF stations (Appendix Fig. A1). The average temperature of the released snowpack was used to infer temperature gradients and snow density in the release zone, further interpolated across stations at varying altitudes to understand snow-cover gradients.

Our approach was to simulate the avalanche events from 2019 based on meteorological data and snow pits and to validate the results for the avalanche outlines and dimensions that we measured based on the post-avalanche drone data. We started by simulating one avalanche and then applying the parameter set of this avalanche to the other two tracks as they were triggered almost at the same time and hence should have the same input data. A summary of the snow cover and temperature input data of the Brämabühl events is presented in Table 1. The same avalanching period in January 2019 produced a well-documented event on the nearby Salezer avalanche track in Davos (on the same day) as well as a powder avalanche at the experimental VdlS test site, which we could additionally use for validation. In this publication, we keep the focus on the Brämabühl event to describe our methodology.

2.2 Evaluation method

Our evaluation approach focused on determining if an avalanche can reach the road. We also estimated the extent of both the dense avalanche core and powder cloud impact pressures. To do so, we have developed a post-processing tool to assess model outputs based on the maximum values per pixel of velocity, flow height, and pressure per calculation cell as described in Glaus et al. (2023). We extract the outlines of both the core and cloud, determining the longest distance by identifying the two most distant points using the convex hull algorithm combined with the rotating calipers method (in Python via scipy.spatial.ConvexHull). The resulting avalanche runout distances for the Brämabühl event are marked in Fig. 13, and the values are given in Table 1. For core outline, we identified pixels with a flow height greater than 0.1 m and a velocity less than 1 m s⁻¹. For the cloud, the outline is based on the stagnation pressure with a lower threshold set at 0.5 kPa computed by $\mathrm{0.5}{\mathit{\rho }}_{\mathrm{\Phi }}{u}_{\mathrm{\Phi }}^{\mathrm{2}}$ with the powder cloud profile described in Zhuang et al. (2023a). While this method works well for simple avalanches, it requires careful consideration in cases where avalanches exhibit finger formation or the avalanche strongly deviates in the lateral direction. Additionally, the code accounts for the case of perfectly symmetric avalanches (as observed on idealized slopes) by measuring the distance between the two outermost points in the flow direction.

We conducted a sensitivity analysis to determine how the avalanche responds to variations in the initial conditions and model parameters. The goal of that study was to quantify the impact of uncertainties in input parameters on simulation outcomes. Given the complexity of avalanches, a large number of model parameters required calibration, but not all data can be presented in detail in this publication due to space limitations.

First we back-calculate the avalanche data from the 2019 Brämabühl event with the given parameter set from the VdlS calibration for all the hard-coded parameters. We vary one input parameter at a time and compare it to the observed avalanche outlines that are presented in Fig. 1. To set the input parameter set, we use the weather station data, with an emphasis on parameters that practitioners can measure. Subsequently, we expand to varying pairs of parameters, such as snow-cover temperature and temperature gradient. Details of these initial findings are discussed in Glaus et al. (2023).

Table 1 Overview of the input data for snow height and temperature data for the Brämabühl events.

Figure 5 Simulated erodible snow cover (a) and temperature distribution (b) applied for the modelling. The observed outlines from the Brämabühl events are marked in red.

2.3 A method for modelling snow-cover distribution d_Σ(Z) and temperature T_Σ(Z)

The underlying idea behind the SLF procedure on avalanche dynamics calculations is to exploit long-term, measured frequency–magnitude snowfall data to determine avalanche fracture heights d₀. Avalanche fracture heights are explicitly related to measured extreme 3 d snow-depth increase (Salm et al., 1990). This procedure underscores two salient assumptions of the Swiss guidelines. Firstly, regional variations in snowfall climatology are included via the measurement data (snow height frequency); secondly, extreme avalanche activity is directly related to intense new snowfall. In the following, we develop a methodology to determine avalanche entrainment heights for road safety calculations within the framework of these existing Swiss guideline procedures that uphold these two basic assumptions.

Because the snow monitoring stations are located at different elevations and different slope angles, snow accumulation data must be adjusted to account for the specific elevation and slope of the avalanche release zone. Within the guidelines, this is performed by applying a snow height gradient ∇D. If Z_m is the elevation of the measurement station and if Z₀ is the elevation of the avalanche release zone, then the first iteration of the fracture height $d_{\mathrm{0}}^{\left(\mathrm{1}\right)}$ is found by adjusting the value obtained, d_m, from the statistical frequency–magnitude analysis of the measurement station (Salm et al., 1990):

$$\begin{array}{}\text{(1)}& d_{\mathrm{0}}^{\left(\mathrm{1}\right)}=d_{\mathrm{m}}+\mathrm{\nabla }D\left(Z_{\mathrm{0}}-Z_{\mathrm{m}}\right).\end{array}$$

The gradient ∇D has units of metres per metre (m m⁻¹). To align with the input values for RAMMS::Extended, we express it as metres per 100 m change in elevation (see Fig. 6). Higher snow accumulation heights are found at higher elevations. Typical gradient values for Switzerland (European Alps) are 0.03 m/100 m $\le \mathrm{\nabla }{d}_{\mathrm{0}}\le$ 0.05/100 m (Salm et al., 1990).

Figure 6 Graphical depiction of input values to initialize a simulation. An example of the values used for reproducing the Rüchi path is shown in (b). The values are chosen based on field measurements and weather stations.

The next iteration $d_{\mathrm{0}}^{\left(\mathrm{2}\right)}$ accounts for the slope angle of the release zone for slopes steeper than 28°. The height $d_{\mathrm{0}}^{\left(\mathrm{1}\right)}$ is adjusted with the slope reduction factor f(ψ) (Salm et al., 1990),

$$\begin{array}{}\text{(2)}& f\left(\mathit{\psi }\right)={\displaystyle \frac{\sin \left(\mathrm{28}\mathit{°}\right)-\mathrm{0.202}\cos \left(\mathrm{28}\mathit{°}\right)}{\sin \left(\mathit{\psi }\right)-\mathrm{0.202}\cos \left(\mathit{\psi }\right)}},\end{array}$$

to calculate avalanche fracture heights:

$$\begin{array}{}\text{(3)}& d_{\mathrm{0}}=(d_{\mathrm{0}}{)}^{\mathrm{2}}=f(\mathit{\psi })d_{\mathrm{0}}^{\left(\mathrm{1}\right)}.\end{array}$$

The slope reduction formula is derived by treating the new snow layer as a Mohr–Coulomb continuum governed by cohesion (c≈600 Pa) and internal friction angle (tan (ϕ) = 0.202) with density ρ₀=200 kg m⁻³ (Salm et al., 1990).

We adopt the same two-step procedure to derive the erodible snow depth for road safety calculations. Moreover, we take

$$\begin{array}{}\text{(4)}& d_{\mathrm{\Sigma }}\left(Z\right)=f\left(\mathit{\psi }\right)\left[d_{\mathrm{m}}+\mathrm{\nabla }D\left(Z-Z_{\mathrm{m}}\right)\right],\end{array}$$

where we now replace the fracture zone elevation Z₀ with the slope elevation Z. This procedure thus places less snow on very steep track segments, e.g., on cliff faces. In the following, we do not take the guideline values for accumulation gradients, but we do take the values derived directly from the measurements ∇D=0.1 m/100 m in the present case; see Table 1 as described in Sect. 2. The resulting snow distribution is shown in Fig. 5.

We see strong temperature gradients present in the upper layer of the snowpack (i.e. recently accumulated snow) in both the release and runout zones (see Fig. A1 in Appendix A for the snow-pit measurements at nearby stations at the altitude of the release zone and deposition). The release temperature is defined as the average temperature of the snowpack from the fracture depth to the surface, excluding the snow surface temperature from the calculation. The measurements indicate that low snow temperatures ($T_{\mathrm{m}}\approx -\mathrm{8.5}$ °C) exist in the avalanche release zone of Rüchi avalanche (Z=2350 m a.s.l.), and higher temperatures exist at the IMIS SLF2 station ($T_{\mathrm{m}}\approx -\mathrm{1.9}$ °C) at Z=1570 m a.s.l. The temperature gradient ∇T=1 °C/100 m (Table 1 and Fig. 5) can be determined from

$$\begin{array}{}\text{(5)}& T_{\mathrm{\Sigma }}\left(Z\right)=T_{\mathrm{m}}-\mathrm{\nabla }T\left(Z-Z_{\mathrm{m}}\right).\end{array}$$

3 Avalanche model

To back-calculate the observed avalanches, we utilize the enhanced version of the depth-averaged RAMMS model (Christen et al., 2010). The extended model encompasses the avalanche core (denoted by the Greek letter Φ), the powder cloud (designated as Π) and the underlying snow cover (Σ) (see Fig. 7). The basics of the model are presented in Zhuang et al. (2023a). In this section, we will revisit some of the key equations from that publication to provide a complete overview. Along with these equations, we will include more detailed explanations and introduce the closure relations.

To accurately model the observed avalanches and snow-cover conditions, the following model features are necessary and contained in the extended RAMMS model: (1) computation of the mean internal energy (thermal temperature) of the avalanche core given the initial temperature of the snow cover; (2) the ability to define snow-cover properties as a function of elevation, exposition, and slope-dependent terrain features; and (3) the ability to track the generation and independent propagation of the powder cloud. Development of the model has been conducted incrementally by Bartelt and Buser, as well as their collaborators (Bartelt et al., 2006; Buser and Bartelt, 2009; Bartelt et al., 2012, 2015a; Zhuang et al., 2023b). A first 1D two-layer model for powder snow avalanches including entrainment was developed at the beginning of the 1980s by Eglit (1983). Further models were developed by Russian researchers (Bozhinskiy and Losev, 1998); however, these models did not include grain flow process physics (Haff, 1983; Hutter et al., 1987; Jenkins and Mancini, 1987) or thermal effects (Valero et al., 2015; Vera Valero et al., 2018). In the following sections, we present the model equations for the core and the cloud, and we discuss the role of the snow cover. The equations are based on the change of height of the snow cover, core, cloud and air to represent the mass flow. As we calculate with the mass divided by the deposit density and unit footprint area, we will show the equations as a function of the height H. The mass flux which is also normalized over density and unit area across the boundary between two layers will be denoted by the velocity $\dot{H}$.

Figure 7 The three primary components of the extended RAMMS model are the avalanche core Φ, the powder cloud Π and the erodible snow cover Σ. The surrounding air is denoted by Λ.

3.1 Avalanche core Φ

The avalanche core is a shear flow containing mass in the form of snow clods (grain flow). The core dynamics are characterized by three state variables, namely the co-volume height ${\widehat{h}}_{\mathrm{\Phi }}$, representing the snow packing found in the deposition zone; the dispersed or flowing height h_Φ; and the slope-parallel velocity vector u_Φ. The co-volume height has an associated density ${\widehat{\mathit{\rho }}}_{\mathrm{\Phi }}$, whereas the dispersed flow height has a density ρ_Φ. The model assumes uniform density and velocity profiles where the mean values are taken. Hence, to estimate damage areas, the profile function must be applied by the user during post-processing. The mass and momentum equations for the avalanche core Φ are

$$\begin{array}{}\text{(6)}& {\partial }_t{\widehat{h}}_{\mathrm{\Phi }}+\mathbf{\nabla }\cdot \left({\widehat{h}}_{\mathrm{\Phi }}{\mathit{u}}_{\mathrm{\Phi }}\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{\Sigma }}}{{\widehat{\mathit{\rho }}}_{\mathrm{\Phi }}}}{\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}-{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Psi }}-{\displaystyle \frac{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}}{{\widehat{\mathit{\rho }}}_{\mathrm{\Phi }}}}{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }},\end{array}$$

$$\begin{array}{}\text{(7)}& {\partial }_t{h}_{\mathrm{\Phi }}+\mathbf{\nabla }\cdot \left(h_{\mathrm{\Phi }}{\mathit{u}}_{\mathrm{\Phi }}\right)=\mathcal{D}(x,y,z),\end{array}$$

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\partial }_t\left({\widehat{h}}_{\mathrm{\Phi }}{\mathit{u}}_{\mathrm{\Phi }}\right)+\mathbf{\nabla }\cdot \left({\widehat{h}}_{\mathrm{\Phi }}{\mathit{u}}_{\mathrm{\Phi }}\otimes {\mathit{u}}_{\mathrm{\Phi }}+p_{\mathrm{\Phi }}I\right)\\ & =\mathit{G}{\widehat{h}}_{\mathrm{\Phi }}-{\displaystyle \frac{{\mathit{u}}_{\mathrm{\Phi }}}{\left|\right|{\mathit{u}}_{\mathrm{\Phi }}\left|\right|}}{\mathit{S}}_{\mathrm{\Phi }}-\left[{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Psi }}+{\displaystyle \frac{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}}{{\widehat{\mathit{\rho }}}_{\mathrm{\Phi }}}}{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}\right]{\mathit{u}}_{\mathrm{\Phi }}.\end{array}\end{array}$$

The avalanche core is driven by gravity G and resisted by shear stress per unit density S_Φ. The mass and momentum balances involve the snow-cover entrainment ${\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}$, snow detrainment by trees ${\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Psi }}$ and mass/momentum transfer to the cloud ${\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}$. As a simplification, the mass exchange between the dense core and the suspension layer is considered unidirectional, with no mass from the cloud being drawn back into the core. The implications of this assumption are discussed in Vicari and Issler (2025). Parametrization of forest detrainment is discussed in detail in Feistl et al. (2015). Equation (7) describes the dilution and compression of the core; $\mathcal{D}(x,y,t)$ represents the change in core height due to dispersive pressure effects (Buser and Bartelt, 2015). We calculate volumes of air that enter the core as the core expands and contracts. Hence, the equation tracks the centre of mass of the granular ensemble. The equation is both a conservation (of air) and an evolution (centre-of-mass) equation. A critical discussion on the approach of modelling the acceleration of the bed-normal expansion of the core's centre of mass directly on the bed-normal component of gravity can be found in Issler et al. (2018). $\mathcal{D}(x,y,t)$ is found by linking the mass and momentum equations to a balance equation for the fluctuation energy (granular temperature) R_Φ; see Haff (1983), Jenkins and Savage (1983), Hutter et al. (1987), Jenkins and Mancini (1987), and Buser and Bartelt (2009) (for more details see Sect. B).

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\partial }_t\left({\widehat{h}}_{\mathrm{\Phi }}{R}_{\mathrm{\Phi }}\right)+\mathbf{\nabla }\cdot \left({\widehat{h}}_{\mathrm{\Phi }}{R}_{\mathrm{\Phi }}{\mathit{u}}_{\mathrm{\Phi }}\right)={\mathit{\alpha }}_{\mathrm{\Phi }}{\dot{W}}_{\mathrm{\Phi }}\\ & -{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}{R}_{\mathrm{\Phi }}-{\mathit{\beta }}_{\mathrm{\Phi }}{\widehat{h}}_{\mathrm{\Phi }}{R}_{\mathrm{\Phi }}+{\mathit{ϵ}}_{\mathrm{\Phi }}{\mathit{\rho }}_{\mathrm{\Sigma }}{\dot{L}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}.\end{array}\end{array}$$

The fluctuation energy R_Φ is associated with random and dispersive particle movements in the flowing granular ensemble. It is produced by shearing ${\dot{W}}_{\mathrm{\Phi }}$ (parameter α_Φ) and decaying by collisions/rubbing (parameter β_Φ) (see Haff, 1983; Jenkins and Savage, 1983; Bartelt et al., 2006). These parameters are temperature dependent. Due to limited data between cold (entrained snow colder than −1 °C) and warm (entrained snow warmer than −1 °C) avalanches (Köhler et al., 2018), a constant α of 0.07 is assumed for cold avalanches and 0.05 for warm avalanches. The calibrated curve for β is shown in Fig. A1 in Appendix A. Random particle movements are likewise produced during the entrainment process at the rate $L_{\mathrm{\Sigma }\to \mathrm{\Phi }}=\mathrm{1}/\mathrm{2}{\dot{H}}_{\mathrm{\Sigma }}{u}_{\mathrm{\Phi }}^{\mathrm{2}}$. The parameter ϵ_Φ defines the linear partitioning of the energy dissipated during the entrainment process that is converted into heat (see Bartelt et al., 2018). The counterpart to the macroscopic random fluctuations is identified as another form of stochastic energy, denoted as the internal energy E_Φ, which is the complementary part of the macroscopic random fluctuations:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\partial }_t\left({\widehat{h}}_{\mathrm{\Phi }}{E}_{\mathrm{\Phi }}\right)+\mathbf{\nabla }\cdot \left({\widehat{h}}_{\mathrm{\Phi }}{E}_{\mathrm{\Phi }}{\mathit{u}}_{\mathrm{\Phi }}\right)=\left[\mathrm{1}-{\mathit{\alpha }}_{\mathrm{\Phi }}\right]{\dot{W}}_{\mathrm{\Phi }}-{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}\\ & E_{\mathrm{\Phi }}+{\mathit{\beta }}_{\mathrm{\Phi }}{\widehat{h}}_{\mathrm{\Phi }}{R}_{\mathrm{\Phi }}+\left[\mathrm{1}-{\mathit{ϵ}}_{\mathrm{\Phi }}\right]{\mathit{\rho }}_{\mathrm{\Sigma }}{\dot{L}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}+{\mathit{\rho }}_{\mathrm{\Sigma }}{c}_{\mathrm{\Sigma }}{T}_{\mathrm{\Sigma }}\\ & {\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}-{\dot{Q}}_{\mathrm{m}}-q_{\mathrm{\Phi }\to \mathrm{\Lambda }}.\end{array}\end{array}$$

The model therefore predicts the mean avalanche temperature T_Φ which is related to the internal energy $E_{\mathrm{\Phi }}={\widehat{\mathit{\rho }}}_{\mathrm{\Phi }}{c}_{\mathrm{\Phi }}{T}_{\mathrm{\Phi }}$, where c_Φ is the specific heat capacity of snow at the density ρ_Φ, and ${\dot{Q}}_{\mathrm{m}}$ is a heat flux. We assume frictional heating processes due to shearing and rubbing between the particles. More details on the implementation of the melting process can be found in Vera Valero et al. (2018) and Zhuang et al. (2023a).

3.2 Flow friction

The inclusion of the state variables (R_Φ, T_Φ) allows us to define a process-based frictional resistance for the avalanche core Φ which is governed by material constants (Table 2). We apply a modified Voellmy-type friction law for flowing snow,

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& S_{\mathrm{\Phi }}\left(R_{\mathrm{\Phi }}\right)={\mathit{\mu }}_{\mathrm{\Phi }}{R}_{\mathrm{\Phi }}{N}_{\mathrm{\Phi }}+(\mathrm{1}-{\mathit{\mu }}_{\mathrm{\Phi }}{R}_{\mathrm{\Phi }})N_{\mathrm{0}}\\ & \left[\mathrm{1}-\exp \left(-{\displaystyle \frac{N_{\mathrm{\Phi }}}{N_{\mathrm{0}}}}\right)\right]+{\mathit{\rho }}_{\mathrm{\Phi }}g{\displaystyle \frac{\left|\right|{\mathit{u}}_{\mathrm{\Phi }}|{|}^{\mathrm{2}}}{{\mathit{\xi }}_{\mathrm{\Phi }}{R}_{\mathrm{\Phi }}}},\end{array}\end{array}$$

where N_Φ is the basal normal stress, μ_Φ(R_Φ) the Coulomb friction coefficient, ξ_Φ(R_Φ) the velocity-dependent Coulomb friction coefficient and N₀ the so-called cohesion (Bartelt et al., 2015b). This empirical formulation was calibrated based on chute experiments with flowing snow (Platzer et al., 2007b, a; Bartelt et al., 2015b). When N₀ = 0, the formula reduces to the traditional Voellmy friction law (Salm, 1993). The formula therefore allows us to exploit, if necessary, the long historical knowledge and well-calibrated sets of Voellmy parameters used by practitioners; see Salm et al. (1990).

The same experiments with flowing snow reveal a strong frictional hysteresis between the front and tail of the flow, indicating a process or flow-dependent relation (Platzer et al., 2007b). Avalanche flow structure, now readily observed in field experiments (Sovilla et al., 2008), is likewise controlled by the frictional hysteresis between front and tail. Moreover, flow resistance at the front of the avalanche differs from the friction at the avalanche tail (Bartelt et al., 2007, 2012). This has significance for the determination of the frictional constants. We note that when R_Φ=0 we have the co-volume or non-dispersive (dense, plug, tail) friction values ${\mathit{\mu }}_{\mathrm{0}}={\mathit{\mu }}_{\mathrm{\Phi }}(R_{\mathrm{\Phi }}=\mathrm{0})$ and ${\mathit{\xi }}_{\mathrm{0}}={\mathit{\xi }}_{\mathrm{\Phi }}(R_{\mathrm{\Phi }}=\mathrm{0})$. The friction is lower at the avalanche front because of the larger random kinetic energy:

$$\begin{array}{}\text{(12)}& {\mathit{\mu }}_{\mathrm{\Phi }}\left(R_{\mathrm{\Phi }}\right)={\mathit{\mu }}_{\mathrm{0}}\exp \left[-{\displaystyle \frac{R_{\mathrm{\Phi }}}{A_{\mathrm{\Phi }}}}\right]{\mathit{\xi }}_{\mathrm{\Phi }}\left(R_{\mathrm{\Phi }}\right)={\mathit{\xi }}_{\mathrm{0}}\exp \left[{\displaystyle \frac{R_{\mathrm{\Phi }}}{A_{\mathrm{\Phi }}}}\right],\end{array}$$

where A_Φ is the so-called activation energy (Bartelt et al., 2012). In this model, the avalanche front dynamics, responsible for the formation of the powder cloud, are mathematically represented as the region of the avalanche with higher fluctuation energies R_Φ. The parameters (μ₀, ξ₀, A_Φ) are found via experiments but, more importantly, by back-calculation of measured avalanche deposits (Bartelt et al., 2012). The spatial distribution of avalanche deposits in the field provides the additional needed information to determine friction parameters. For example, they can be immediately estimated in the field by noting the steepest slope ψ with avalanche snow tan (ψ)≈μ₀, as deposition begins when R_Φ→ 0. The location of the frontal deposits (runout) and the terminal velocity of the avalanche are necessary to calibrate ξ₀ and the activation energy A_Φ. Table 2 lists the recommended frictional values we take for avalanching after 3 d of new snowfall accumulations.

Platzer et al. (2007b)Bartelt et al. (2012)Zhuang et al. (2023b)Bartelt et al. (2015b)Bartelt et al. (2012)Dreier et al. (2016)Zhuang et al. (2023b)Bartelt et al. (2012)

Table 2 Overview of the fixed parameters in RAMMS::EXTENDED with the corresponding publications where the calibrations are described in more detail.

For our present purposes to investigate cold, mixed flowing avalanches appearing after new snowfall periods, we will take the model parameters (μ₀, ξ₀, N₀, A_Φ, α_Φ) to be temperature-independent constants (Table 2). The only temperature-dependent parameter will be the decay of fluctuation energy β_Φ(T_Φ) with T_Φ (given in K). With the arctangent relationship, we model it as

$$\begin{array}{}\text{(13)}& {\mathit{\beta }}_{\mathrm{\Phi }}\left(T_{\mathrm{\Phi }}\right)=\mathrm{1.40}+{\displaystyle \frac{\mathrm{1.6}}{\mathit{\pi }}}\arctan \left(\mathrm{1.6}\right(T_{\mathrm{\Phi }}-\mathrm{271.5}\mathit{°}K\left)\right)\end{array}$$

to ensure that the decay coefficient is within the range 0.6/ s $\le {\mathit{\beta }}_{\mathrm{\Phi }}\le \mathrm{2.0}$/s. As the inverse of β_Φ physically represents the lifetime of the fluctuation energy R_Φ, it is linked to the onset of deposition and the flow structure of the avalanche (formation of the avalanche tail). It can therefore be determined by measuring the distribution of deposits in the avalanche runout zone (Bartelt et al., 2012). The lifetime of the fluctuation energy decreases as the avalanche temperature increases; it is approximately 4 times longer in a cold avalanche than a warm avalanche. This ensures that warm, moist avalanches have plug-like flow regimes (Köhler et al., 2018).

Finally, we presently do not consider the influence of generated meltwater on the frictional constants, assuming that the snow temperature remains below the melting temperature of ice.

3.3 Powder cloud Π

A comparable set of partial differential equations is proposed to model the powder cloud Π. The powder cloud is simulated by equations governing mass (Eqs. 14 and 15) and momentum balance (Eq. 16), along with supplementary equations related to the generation and dissipation of turbulent fluctuations (Eq. 17):

$$\begin{array}{}\text{(14)}& {\displaystyle }{\partial }_t{\widehat{h}}_{\mathrm{\Pi }}+\mathbf{\nabla }\cdot \left({\widehat{h}}_{\mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Pi }}\right)={\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }},\text{(15)}& {\displaystyle }{\partial }_t{h}_{\mathrm{\Pi }}+\mathbf{\nabla }\cdot \left(h_{\mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Pi }}\right)={\dot{H}}_{\mathrm{\Lambda }\to \mathrm{\Pi }}+{\displaystyle \frac{{\mathit{\rho }}_i-{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}}{{\mathit{\rho }}_i-{\mathit{\rho }}_{\mathrm{\Lambda }}}}{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }},\text{(16)}& {\displaystyle }\begin{array}{rl}& {\partial }_t\left({\widehat{h}}_{\mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Pi }}\right)+\mathbf{\nabla }\cdot \left({\widehat{h}}_{\mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Pi }}\otimes {\mathit{u}}_{\mathrm{\Pi }}+p_{\mathrm{\Pi }}\mathbb{I}\right)={\displaystyle \frac{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}-{\mathit{\rho }}_{\mathrm{\Lambda }}}{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}}}\\ & +\mathit{G}{\widehat{h}}_{\mathrm{\Pi }}{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Phi }}-{\displaystyle \frac{{\mathit{u}}_{\mathrm{\Pi }}}{\left|\right|{\mathit{u}}_{\mathrm{\Pi }}\left|\right|}}{\mathit{S}}_{\mathrm{\Pi }}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{\Lambda }}}{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}}}{\dot{H}}_{\mathrm{\Lambda }\to \mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Lambda }},\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}& {\partial }_t\left({\widehat{h}}_{\mathrm{\Pi }}{R}_{\mathrm{\Pi }}\right)+\mathbf{\nabla }\cdot \left({\widehat{h}}_{\mathrm{\Pi }}{R}_{\mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Pi }}\right)={\dot{W}}_{\mathrm{\Pi }}+{\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}{R}_{\mathrm{\Phi }}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{\Lambda }}{\dot{H}}_{\mathrm{\Lambda }\to \mathrm{\Pi }}\left|\right|{\mathit{u}}_{\mathrm{\Pi }}|{|}^{\mathrm{2}}-{\mathit{\beta }}_{\mathrm{\Pi }}{\widehat{h}}_{\mathrm{\Pi }}{R}_{\mathrm{\Pi }}.\end{array}\end{array}$$

The initial cloud height is denoted as ${\widehat{h}}_{\mathrm{\Pi }}$ with the corresponding initial cloud density, ${\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}$, representing the cloud state before the expulsion from the core. The variable h_Π denotes the actual cloud height which results from dust–air mixture expelled from the core ${\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}$ and air entrainment ${\dot{H}}_{\mathrm{\Lambda }\to \mathrm{\Pi }}$. The cloud density decreases due to air entrainment denoted by ρ_Π. The actual cloud density must fulfil the following relation: ${\mathit{\rho }}_{\mathrm{\Pi }}={\mathit{\rho }}_i\frac{{\mathit{\varphi }}_i{\widehat{h}}_{\mathrm{\Pi }}}{h_{\mathrm{\Pi }}+{\mathit{\varphi }}_i{\widehat{h}}_{\mathrm{\Pi }}}+{\mathit{\rho }}_{\mathrm{\Lambda }}\frac{h_{\mathrm{\Pi }}}{h_{\mathrm{\Pi }}+{\mathit{\varphi }}_i{\widehat{h}}_{\mathrm{\Pi }}}$, where ρ_i=917 kg m⁻³ is the ice density, ρ_Λ=1.225 kg m⁻³ is the air density and ${\mathit{\varphi }}_i=\frac{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}-{\mathit{\rho }}_{\mathrm{\Lambda }}}{{\mathit{\rho }}_i-{\mathit{\rho }}_{\mathrm{\Lambda }}}$ represents the ice fraction in the initial cloud. The density of air varies with changes in temperature, pressure and humidity, but this variation has a negligible effect on the simulation. Therefore, we have chosen to keep the value fixed. The cloud is propelled by the momentum imparted from the core ${\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Phi }}$ and gravity $\frac{{\mathit{\rho }}_{\mathrm{\Pi }}-{\mathit{\rho }}_{\mathrm{\Lambda }}}{{\mathit{\rho }}_{\mathrm{\Pi }}}{\mathbf{G}}_{\mathrm{\Pi }}$. Generally, we observe ${\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}{\mathit{u}}_{\mathrm{\Phi }}\gg \frac{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}-{\mathit{\rho }}_{\mathrm{\Lambda }}}{{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}}\mathit{G}{\widehat{h}}_{\mathrm{\Pi }}$. We neglect the gravitational component even if we agree that $\mathit{G}{\widehat{h}}_{\mathrm{\Pi }}$ can become significant for a large powder cloud. This simplification remains valid as long as the dense core and the suspension layer move together. However, once the core comes to a stop and the cloud ascends a counterslope, the gravitational term will once again become more significant, as discussed in Issler et al. (2018). Additionally, it is important to note that the mass flow between the cloud and core is assumed to be unidirectional. Consequently, the vertical expansion of the core, driven by increased granular temperature, prevents the formation of a vacuum within the core. As a result, a portion of the air–snow mixture from the cloud is drawn back into the core (Vicari and Issler, 2025).

The turbulence in the cloud due to the fluctuation energy is described by Eq. (17). The structure of the equation is based on earlier one-layer 3D models (Hermann et al., 1994; Gauer, 1995; Sampl and Zwinger, 2004). The fluctuation energy is produced by three sources (Eq. 17): internal shearing $W_{\mathrm{\Pi }}=\left[{\widehat{\mathit{\rho }}}_{\mathrm{\Pi }}{S}_{\mathrm{\Pi }}\right]\Vert {\mathit{u}}_{\mathrm{\Pi }}\Vert$, fluctuation energy transferred from the core ${\dot{H}}_{\mathrm{\Phi }\to \mathrm{\Pi }}{R}_{\mathrm{\Phi }}$ and air entrainment $\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\rho }}_{\mathrm{\Lambda }}{\dot{H}}_{\mathrm{\Lambda }\to \mathrm{\Pi }}{u}_{\mathrm{\Pi }}^{\mathrm{2}}$. β_Π is the parameter that controls the decay of turbulence and, therefore, the lifetime of the fluctuation energy $-{\mathit{\beta }}_{\mathrm{\Pi }}{\widehat{h}}_{\mathrm{\Pi }}{R}_{\mathrm{\Pi }}$. The pressure p_Π includes both the hydrostatic and turbulent parts. More details of the powder model equations, including the entrainment function ${\dot{H}}_{\mathrm{\Lambda }\to \mathrm{\Pi }}$, friction S_Π and turbulence parameters is contained in the publication (Zhuang et al., 2023b) and in Appendix B. The air entrainment function and friction are calibrated based on observed avalanches in VdlS and shown in Appendix B.

In the cloud, we assume inelastic collisions between particles, such that all energy is converted into random motion, with no energy dissipated as heat, as a simplification. Hence, the constant ϵ_Π is set at 1.

3.4 Entrainment

Entrainment in the extended model equations is treated as a plastic collision between the avalanche core and the snow cover (Bartelt et al., 2018). We initially define the snow volume per unit footprint area and unit time, which results in a velocity ${\dot{H}}_{\mathrm{\Sigma }}$ that is affected by the passage of the avalanche core:

$$\begin{array}{}\text{(18)}& {\dot{H}}_{\mathrm{\Sigma }}={\mathit{\kappa }}_{\mathrm{\Sigma }}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{\Sigma }}}{{\widehat{\mathit{\rho }}}_{\mathrm{\Phi }}}}\left|\right|{\mathit{u}}_{\mathrm{\Phi }}\left|\right|.\end{array}$$

We represent this interaction rate as proportional to the avalanche speed $\left|\right|{\mathit{u}}_{\mathrm{\Phi }}\left|\right|$, as it determines the distance the avalanche travels during the interaction time. Avalanches moving at higher speeds cover more ground, leading to an increase in the amount of snow-cover mass affected by the avalanche. In this model, even very thin avalanches could potentially erode the same amount of snow as a thick one (Issler et al., 2024). To prevent this, a shear stress cutoff is implemented, ensuring that an avalanche must exceed a certain energy threshold before it can erode snow.

We define the erodibility coefficient as the dimensionless parameter κ. Essentially, κ determines the ratio of erosion speed of the erosion front in the snow cover in direction normal to the terrain relative to the flow velocity. Low κ values indicate that only the surface of the snow cover is affected by the avalanche passage (basal erosion), while, conversely, high values of κ indicate that the core affects the entire depth of the snow cover (frontal erosion). The value of κ can be adjusted to incorporate an internal friction coefficient of the snow μ_b. Defining g_s to be the slope parallel acceleration and g_z as the slope normal acceleration, we modify κ to be

$$\begin{array}{}\text{(19)}& \mathit{\kappa }={\displaystyle \frac{{\mathit{\kappa }}^{\prime }}{g}}\left[g_{\mathrm{s}}-{\mathit{\mu }}_{\mathrm{b}}{g}_z\right]\text{, where}\mathit{\kappa }\ge \mathrm{0}\text{always}.\end{array}$$

This model provided good results, even though we recognize the limitation of no entrainment below a limiting slope angle for cases as in flat terrain, in level terrain or on counterslopes. The avalanche observed in a VdlS avalanche on 10 February 1999 (Gruber and Margreth, 2001) shows, for example, a case where an avalanche eroded snow in a flatter part. For 3 d accumulation periods with new snow, we take ${\mathit{\kappa }}^{\prime }=\mathrm{0.015}$ and μ_b≈0. These values ensure that on steep slopes with high avalanche velocities we model frontal entrainment, whereas on more gentle slopes and at lower speeds the avalanche enters a mode of basal erosion. The bonding strength model is motivated by observations of eroded segments in avalanche tracks. On track segments where there are no depositions, it can be ascertained that erosion has occurred. In this case, the parameter μ_b must be smaller than the tangent of the slope angle.

Presently, ${\dot{H}}_{\mathrm{\Sigma }}$ represents the change in snow-cover mass due to the avalanche core – not the total amount of snow taken in by the avalanche. We now partition the affected mass into two parts: a part of mass which is entrained by the avalanche and a part of mass which is not entrained, i.e. possibly splashed in front of the core to build a pre-front or frontal saltation layer Γ:

$$\begin{array}{}\text{(20)}& {\dot{H}}_{\mathrm{\Sigma }}={\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}+{\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Gamma }}.\end{array}$$

The mass flux ${\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}$ represents the snow-cover mass that is accelerated to the avalanche velocity and can be found on the right-hand side of model equations (Eqs. 6 and 8). We apply a partitioning parameter γ to separate the entrained/non-entrained fractions of the snow cover:

$$\begin{array}{}\text{(21)}& {\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Gamma }}=\mathit{\gamma }{\dot{H}}_{\mathrm{\Sigma }}{\dot{H}}_{\mathrm{\Sigma }\to \mathrm{\Phi }}=(\mathrm{1}-\mathit{\gamma }){\dot{H}}_{\mathrm{\Sigma }}.\end{array}$$

The parameter γ, which we term the splashing parameter, could also represent the non-entrained mass in the disrupted snow cover that is simply accelerated by the passage of the avalanche front. Different snow covers will be governed by different entrainment parameters (κ^′, μ_b, γ). For 3 d of new snow accumulation periods, we take γ=0.2 as a simplification. For moist avalanches, this value goes to zero.

4 Results and discussion

In the following, we present the results of the sensitivity analysis on model input parameters. We delve deeper into the mathematical model’s representation of the effects of snow temperature and erosion (Sect. 4.1), release-zone properties (Sect. 4.2) and friction parameters (Sect. 4.3). We conclude with a comparison to the measured avalanches (Sect. 4.4). All other system parameters were calibrated based on observations of avalanches in VdlS as presented in Sect. 3. As a base, we used the model parameters shown in Table 2.

4.1 Snow temperature and erodible snow depth

In most avalanching situations, alpine snow is within a few degrees of its phase transition point (T=0 °C). The physical properties of flowing snow undergo rapid transformations as temperatures edge towards the melting threshold. In the model equations, the temperature dependence is contained in the decay parameter of granular temperature β (Eq. 13), which describes the decay of the random kinetic energy (granular temperature) and therefore the dispersion of the snow granules. The parameter is set such that colder avalanches exhibit the tendency to form mixed flowing avalanches (Bartelt et al., 2012), while warmer avalanches will exhibit more plug-type flows (Li et al., 2021). By defining the avalanche release temperature, not only do we set the initial thermal energy, but we also dictate the predisposition of the avalanching snow towards dry, mixed flowing avalanche or moist flow regimes.

Figure 8 Simulation of runout distance of an avalanche on an idealized plane. (a) Variation of maximum erosion depth for different erosion gradients. (b) Variation of release temperatures. Simulation points for release temperatures near the melting point should be interpreted with caution, as the simulation tool is not yet well calibrated for this range.

In a first series of numerical experiments, we applied the model to an idealized slope (Fig. 4). The slope inclination was set to approximate the Brämabühl slopes under investigation. In this way, secondary terrain features inducing flow channels and secondary flow fingers could be removed from the analysis and model performance gauged in idealized conditions. We varied release temperatures from extremely cold temperatures to the melting point (−20 °C ≤ T₀≤ 0 °C).

In the first simulations, we included no entrainment, $d_{\mathrm{0}}^{*}=\mathrm{0}$, meaning no additional snow was eroded by the avalanche. We calculated the runout distance according to our post-processing procedure of minimum heights and velocities (green dots, Fig. 8a). In this case, the runout distances remained constant over a wide temperature range (−20 °C ≤ T₀≤ −8 °C). At higher temperatures the decay parameter of granular temperature increases, with the resulting effect of reduced fluidization and shorter runout distances.

By incorporating erosion in our simulations, we add entrained mass to the flowing avalanches, and their responses become more intricate. In a subsequent series of simulations, we assume a deep snow cover, $d_{\mathrm{0}}^{*}=\mathrm{1.9}$ m, with a small snow-depth gradient ∇D=1 m/100 m and temperature gradients ∇T=0.01 °C/100 m (blue dots, Fig, 8 b). At extremely cold temperatures of −20 °C ≤ T₀≤ −13 °C, the calculated runout distances are independent of T. We assume the reason for this behaviour is the decay parameter of granular temperature β. For very cold avalanches, the formation of the powder cloud dominates, extracting mass and energy from the core and ultimately leading to dispersion and dissipation of the avalanche. As the temperatures increase, an optimal balance between the core and cloud emerges, yielding far-reaching flows. With the specification of more realistic release temperatures, $T>-\mathrm{13}$ °C, the runout distances increase. This phenomenon underscores the counteracting effects of frictional heating (rise in temperature) and the entrainment of cold snow (fall in temperature). The avalanche temperature remains lower for longer, fostering long-lasting fluidized regimes and more potent powder avalanches. The model predicts that the entrainment of cold snow at lower elevations sustains the fluidized regime and the formation of powder avalanches. However, the runout distances decrease again with higher release and entrainment temperatures, $T_{\mathrm{0}}>-\mathrm{8}$ °C (Fig. 8b). At higher temperatures, the decay of fluctuation energy increases, leading to dense, less fluidized flows and therefore an increase in friction which curbs runout distances. It becomes evident that the temperature-dependent decay parameter of granular temperature, β(T_Φ), controls the flow regime of the simulated avalanches with and without entrainment.

In the next series of simulations on the idealized slope, we vary the snow height gradient ∇D=0.0, 0.05 and 0.10 m/100 m; the release temperature remains set at $T_{\mathrm{0}}=-\mathrm{6}$ °C, and the temperature gradient is zero, i.e. ∇T=0. Less snow is subsequently encountered by the avalanche at lower elevations and into the runout zone. We perform simulations with different values of snow heights, 0.0 m $\le d_{\mathrm{0}}^{*}\le$ 2.0 m. This situation mirrors actual snow conditions in road safety applications, and there is a large uncertainty in the snow distribution. We simulate the effect of entrained snow of different height gradients for shallow and deep snow covers. The results indicate that for shallow snow covers ($d_{\mathrm{0}}^{*}\le \mathrm{1.0}$ m) the different gradients produced large differences in avalanche runout (Fig. 8a). Avalanches that encounter snow along the entire track (∇D=0.0 m/100 m) run longer than those that encounter regions of no snow in the runout zone. This result corresponds well with the experience that a deep snow cover from initiation to runout is needed for extreme avalanche events, especially powder snow avalanches. For the case of $d_{\mathrm{0}}^{*}>\mathrm{1.0}$ m, the gradients appear to have no influence on the avalanche runout; the snow is so deep that even a strong gradient still does not entirely remove the snow cover for this parameter set.

Figure 9 Simulation of avalanche runout on the idealized slope with and without erosion. We consider three initial temperatures of (a) $T_{\mathrm{0}}=-\mathrm{6}$ °C, (b) $T_{\mathrm{0}}=-\mathrm{9}$ °C and (c) $T_{\mathrm{0}}=-\mathrm{12}$ °C for 21 different temperature gradients, ∇T. The colours on the bars depict the snow temperature varying from cold (blue) to warm (green).

Next, we investigated the effects of elevational snow temperature gradients in Fig. 9. We show the decrease in avalanche runout for variable temperature gradients ∇T as we simulated it on the idealized slope. The colours in each panel depict the core temperature at various time steps throughout the avalanche flow. The analysis is done for different release temperatures. The results indicate the reduction in avalanche runout as a cold avalanche runs into a warm snow cover. The colder the initial temperature, the smaller the reduction. The results underscore the complex interplay between initial release conditions, the process of entrainment and the responding avalanche temperature.

Figure 10 Temperature evolution in the Rüchi avalanche over time from t=0 s (a) to t=60 s (e). The initial temperature distribution is presented in Fig. 5 (map source: Federal Office of Topography).

The simulations on the idealized slope are useful because they highlight how the various snow parameters play an integral role for the avalanche outcome without localized terrain effects. On an actual slope the snow-cover height is given not only by the elevation gradient ∇D but also by the local slope inclination. The simulated evolution of avalanche temperature over time for the Rüchi path is depicted in Fig. 10. The initial temperatures are specified using the weather station and snow-pit data. This figure displays the calculated mean avalanche temperature T_Φ of the avalanche core, which increases by 6 °C from $T_{\mathrm{\Phi }}=-\mathrm{8}$ °C to $T_{\mathrm{\Phi }}=-\mathrm{2}$ °C in the runout zone. The calculated avalanche temperature results from the competition between frictional heating (avalanche velocity and therefore terrain) and intake of cold snow (Valero et al., 2015). The model equations assume that the entrained snow mixes with the avalanche snow instantaneously, producing a new mean temperature. In reality, the energy exchange between the avalanche snow and the entrained snow will happen fast but takes some time, as shown in Köhler et al. (2018). Temperature variations will exist in the avalanche; heat concentrations will most likely exist on the surface of the granules, while the interior of the granules remains cold (Jomelli and Bertran, 2004; Steinkogler et al., 2015).

4.2 Release zone

In this section, we investigate how the shape and location of the defined release zone influence avalanche simulations by varying the length and steepness of the release zone. We conduct our simulations on the idealized slope, as shown in Fig. 4, but with an infinite runout of 28° slope to remove the influence of the flat runout zone. For all simulations, we used the parameter set for the Rüchi avalanche as described in Table 1.

Figure 11 Representation of the avalanche runout length if the release length and the release angle are varied for the core and the cloud. The release volume is represented in grey.

We used an idealized version of the Rüchi release zone which is a square with a width of 120 m and a length in the flow direction of 75 m. We varied the length up to 2.5 times the width to stay in the ratio of avalanche width compared to length proposed by the Swiss guidelines and the release shape proposed by McClung (2009). The release zones have an average release angle of 35°. As the release volume increases, RAMMS proposes a different volume-dependent friction parameter set for μ and ξ. As we want to analyse RAMMS with its proposed calibration, we plot the values with the proposed friction values according to the release volume and do not hold them constant. By increasing the release length, we automatically increase the release volume, and the avalanche starts with a higher potential energy (see Fig. 11). The simulations show that we get longer avalanches with an increased release zone length.

For the variation in release angle, we used the same release zone throughout. In RAMMS::EXTENDED, the release zone is defined by its projected area, meaning that the actual release area becomes slightly larger as the release angle increases. Figure 11 shows that the runout distance increases almost linearly with an increase in release angle. Avalanches initiating at lower release angles show a powder cloud that is shorter in length compared to the core. Conversely, at release angles exceeding 48°, the powder cloud can achieve greater acceleration, flowing further than the core. On real terrain, the terrain roughness (caused by, for example, vegetation and rocks) would strongly decelerate the avalanche core and reduce its runout distance.

4.3 Friction parameters μ₀, ξ₀ and N₀

In the preceding section, all calculations were performed with constant friction and process parameters (see Table 2) representing avalanche situations governed by periods of new snowfall. The simulation results indicate that given an initial release location, mass and temperature, the calculated terminal velocity and avalanche runout are governed by snow-cover disposition and temperature. Traditionally, runout and velocity are reproduced in avalanche dynamics calculations by changing the values of the friction parameters from avalanche to avalanche, creating an envelope of extreme values (Gruber and Bartelt, 2007). Here, we do not adopt this approach. Friction parameters change dynamically as a function of temperature according to the following process chain:

$$\begin{array}{}\text{(22)}& T_{\mathrm{\Phi }}\left(t\right)\to R_{\mathrm{\Phi }}\left(t\right)\to \text{Voellmy parameters}\mathit{\mu }\left(t\right),\mathit{\xi }\left(t\right).\end{array}$$

This chain of relations indicates that the temperature of snow influences the mean fluctuation energy (via the decay parameter β(T_Φ)), which controls the dispersion of snow granules and therefore the avalanche flow regime. This fluctuation energy, being a stochastic variable, signifies the inherent randomness in the movement of all granules within the flowing snow ensemble. The momentary state of friction is influenced by this fluctuation energy. Thus, the basic model assumption is that the temperature of snow governs the stochastic dynamics of its granular ensemble, ultimately impacting frictional behaviour. The grain flow process parameters controlling the relationship between R_Φ(t) and flow friction (avalanche deposition) have been identified by Bartelt et al. (2012) in the analysis of experimental data from the test site Vallée de la Sionne.

It is now necessary to validate this approach using the Brämabühl avalanches using in situ information of snow distribution and weather station records. For this, the initial values of the friction parameters μ₀ and ξ₀ were varied to identify the combination which results in the measured avalanche runout distance. This was also done for different cohesion values. The results are depicted in Fig. 12. For some cohesion values, the measured runout distance could not be reproduced. Therefore, we restricted the investigated cohesion values to the range obtained from measurements conducted in snow chutes (Bartelt et al., 2012).

Figure 12 Simulation of the Rüchi avalanche. (a) Different friction values and a cohesion of N₀=150 Pa. The measured runout length of 1545 m ±5 % error (over and underprediction of the avalanche runout) is marked by the red rectangle. The best-fit parameters are near the recommended values in Table 2. (b) Evolution of avalanche volume along the Rüchi path from the simulation.

We used the observed runout distance to determine the optimal friction parameters. We find that the best fit of all three Brämabühl avalanche runout distances is provided by friction values of $\mathrm{0.50}\le {\mathit{\mu }}_{\mathrm{0}}\le \mathrm{0.55}$ and 1750 m s⁻² $\le {\mathit{\xi }}_{\mathrm{0}}\le \mathrm{2200}$ m s⁻². This is in good agreement with values found in VdlS calibrations (Bartelt et al., 2012). The range of ξ₀ could be reduced by knowing the avalanche core velocity. In Fig. 12, the Rüchi path exhibits the most symmetric pattern in terms of friction values. This is attributed to the fact that the Rüchi path features a straightforward avalanche outline without any flow fingers.

4.4 Comparison to measured avalanches

Photogrammetric data collected with the drone provided important information on the lateral flow width of the avalanche, the avalanche volume, and the height and travel distance of the powder cloud.

Figure 13 depicts the calculated extent of the model avalanches in comparison to the measured outlines. The relatively good agreement between the calculated and measured avalanche flow widths is particularly significant.

Figure 13 Comparison of the measured outlines of the Brämabühl avalanches with the model results using the measured snow-cover height and temperature data (map source: Federal Office of Topography).

Each avalanche was accompanied by a powder cloud that ran up the counterslope (Fig. 14). The calculated powder cloud widths are in good agreement with the observations. By comparing the photographs of the fully developed powder cloud with the tree heights known from a lidar-based vegetation height model and by knowing that maximum powder cloud pressures on the road never exceeded 5 kPa (as houses and windows remained undamaged), we estimate that the maximum powder cloud heights reached up to 40 m. In the simulation, the powder cloud reaches a maximum height of 50 m already in the very first 100 m of elevation drop, which surpasses expected real-world values. Regions of isolated tree damage by the powder cloud are reproduced by the model. Additionally, comparing the estimated powder cloud velocities (Fig. 2) for the Rüchi and Chaiseren paths in the runout zones shows that they are in the same range as the average velocity of the front of the powder cloud between the pictures resulting from the simulation as shown in Table 3. The simulated velocity consistently overestimates the measured powder cloud front velocity, and the deceleration of the powder cloud is not accurately captured.

Table 3 Comparison of the computed powder cloud front velocity with the simulated velocity for the example of the Rüchi avalanche. The purpose of this comparison is to verify that the results fall within a reasonable range rather than to suggest that the simulation accurately resolves velocity at such a detailed level.

Figure 14 Photo of the cloud taken from the helicopter of the Rüchi avalanche and the simulation results from the back-calculation, showing the maximum cloud height, pressure, velocity and density. All values shown are maximum values per pixel (map source: Federal Office of Topography).

In Fig. 15, we compare the snow height in the deposition zone of the Rüchi path to the simulated deposition height. In the deposition area, on the right side in the flow direction of the avalanche, a higher deposition pile has formed, which is also evident in the simulation. Photogrammetric snow-depth mapping with the drone revealed that the total snow volume present in the area overflowed by the avalanche in Fig. 15 was approximately 155 000 m³. The simulated deposition in the same area was only approximately 30 000 m³. The drone measurements assess the total volume of snow present at acquisition time. As the avalanche ran on an already present snow cover of approximately 1 m (measured next to the avalanche deposit), we need to subtract an estimation of the compressed old snow in the deposition zone. We estimate an additional volume of approximately 100 000 m³. Furthermore, a large portion of the deposition in the simulated avalanche is already deposited further up in the avalanche track, and the front stopped just outside the mapped avalanche outline and is therefore not taken into account for the volume calculation. Additionally, the trees along the avalanche path cause shading effects that lead to gaps in the snow height measurements. Qualitatively, the deposits simulated in the upper part of the track are also present in the drone measurements and the simulations at similar locations (Fig. 13) and show similar deposition heights.

Figure 15 Comparison of the deposition area of the core measured from the drone data (a) and the modelled deposition (b) (map source: Federal Office of Topography).

The drone orthophoto of the Chaiseren avalanche track showed that its snowpack had been scoured by wind, resulting in less accumulated snow than was calculated by the snow gradients used in our simulations. Therefore, the initial snow mass is overestimated in this simulation.

5 Conclusions

For road safety managers, it is essential to know, based on near real-time weather and snowpack information, whether an avalanche could reach the road and a road closure is warranted. With the goal to better support their decision-making process, we tested an approach where we adopted a numerical model (RAMMS::EXTENDED) to investigate the sensitivity of the model predictions to various snow parameters that can be measured in the field, such as, release zone location and length, release temperature, and gradients for snow temperature and snow height on the path, as well as the avalanche response to changes in these parameters. The study presented herein was based on an avalanche cycle from the Dischma valley, Switzerland, in 2019 that caused several avalanches to cover the road during a post-storm cold spell. These avalanches were well documented and could therefore be used for avalanche back-calculations. For the purpose of this study, we only simulated cold-powder avalanches; both the core and the powder cloud of these mixed avalanches could pose a hazard to cars and people. The sensitivity analysis was mainly based on simulations on an idealized plane, resembling the topographic profile of the Rüchi avalanche path. The influence of terrain effects, e.g. friction parameters in the model, was investigated in a second phase of the analysis where we reconstructed the avalanche events on Wildi, Rüchi and Chaiseren paths that took place in 2019 and compared the simulated results to drone observations. In our simulations, we incorporated snow entrainment by varying the snow characteristics on the path, mainly snow height and temperature; we demonstrated the complex relationships between initial snow conditions (volume and temperature); the various processes that take place during avalanche flow (frictional heating, entrainment, etc.); and the resulting flow behaviour, for which runout length and impact pressure are of highest concern for road safety managers. We also demonstrate how various flow regimes may develop and how some interesting flow behaviours may be expressed, based on changes in the snowpack on the path.

The applied model continues to utilize a Voellmy-based frictional approach as in the well-established models applied for hazard mapping. However, the friction coefficients are now dynamically calculated and are affected by terrain and snowpack variables, which differ considerably between different avalanche tracks and avalanche periods. To set up the simulations, we applied snowpack parameters measured at nearby automated weather stations and snow profiles. Further, we utilized photogrammetrically measured snow-depth distributions acquired by drones that were able to capture the extreme spatial variability of snow-depth distribution in our avalanche paths of interest.

The results indicate that we can use measurements from weather stations at different altitudes and locations to calculate the snow-cover distribution and snow temperature gradient and let those inform our model to simulate realistic avalanches. Choosing weather stations from a nearby valley (approximately 3 km distance) showed acceptable results with a coefficient of variation below 5 %. Comparing the modelled avalanche outlines to the observed ones, it is visible how the model represents important features such as the evolution of fingers, indicating zones with higher impact pressures or the development of the powder cloud.

The analysis presented in this publication used version 2.8.28 of RAMMS::EXTENDED. As the model will be developed and calibrated further, the presented parameter set must be treated as provisional and can only be applied to this version. For more current versions, the publication by Stoffel et al. (2024) will give more insights.

The presented approach will now be applied to calculate avalanche runout for different representative weather and snowpack scenarios for the Dischma road. These results will then be evaluated by local experts, and the applicability of this approach for future decision-making will be assessed. A probabilistic approach is currently being tested to calculate reach probabilities to the road for specific avalanche tracks. These are important steps towards a more data-based decision-making for road management in mountain regions. In the future, this concept could also be applied to assess additional areas which are endangered by avalanches, such as houses built in hazard zones and ski slopes.

6 Outlook

Building on the foundations of this project, we aim to develop a tool that integrates real-time data with numerical simulations to predict the daily avalanche runout distance. To achieve this objective, the following directions can be pursued in future work.

As nearby and representative weather station measurements are not accessible for many roads, in a next step snow-cover models such as SNOWPACK (Lehning et al., 1999) and CROCUS (Vionnet et al., 2012) could be used to calculate the relevant input parameters. For rough estimations of the snow temperature, it would be necessary to compare the cloud coverage during the hours before the avalanche occurred to the average temperature of the release snow mass. As the runout distance is not too sensitive to the snow temperature, a first approach could be to define simulation scenarios of cold and warm temperatures depending on the cloud coverage.

To enable the model to simulate additional wet-snow avalanches, we need to collect more avalanche data that are directly related to measured snowpack temperature and moisture content.

In the current project, we worked with release zones that we inferred from the drone data. In practice, we would additionally need to make an expert guess on where the potential release areas could be. This is an important source of uncertainty. A first approach could be to use automatically delineated release areas as proposed by Bühler et al. (2022) and Issler et al. (2023).

To advance the development of a practical tool, additional factors beyond avalanche runout distance must be considered. First, the system should incorporate the probability of avalanche release under daily conditions, with an initial approach being to integrate the avalanche danger level into the calculations. Second, our analysis highlights the significant influence of snow entrainment on avalanche dynamics. Observations from past avalanche events (Issler, 2020) indicate that accurately estimating the potential erodible snow is essential. In a simplified model, erodible snow height could be approximated by identifying weak layers in the SNOWPACK model.

Furthermore, a method for displaying relevant data to experts must be designed. Key questions will arise, such as how to establish thresholds for acceptable powder cloud pressure or snow depth on roads, as well as the acceptable probability levels for road openings.

Appendix A: Snow-pit data

Figure A1 Snow-pit data for the stations at Weissfluhjoch (a) and at SLF Davos (b) (Source: https://whiterisk.ch produced by SLF, last access: 26 March 2025).

Appendix B: Closure relations

In this section, we show the closure relations which result from calibrations based on data from Vallée de la Sionne. We want to point out that these relations represent the best fit with the data we have measured so far.

Decay constant

$$\begin{array}{}\text{(B1)}& \mathit{\beta }=\mathrm{1.40}+{\displaystyle \frac{\mathrm{1.6}}{\mathit{\pi }}}\cdot \arctan \left[\mathrm{1.6}\left(T_{\mathrm{\Phi }}-\mathrm{271.5}\mathrm{K}\right)\right]\end{array}$$

Entrainment

The snow erosion rate is parameterized as follows:

$${\dot{H}}_{\mathrm{\Sigma }}=k_{\mathrm{0}}{k}_T{k}_{\mathit{\psi }}\Vert U_{\mathrm{\Phi }}\Vert ,$$

with the parameters defined as,

$$\begin{array}{}\text{(B2)}& {\displaystyle }{k}_{\mathrm{0}}=\mathrm{0.005},\text{(B3)}& {\displaystyle }{k}_T=\mathrm{3.00}-\mathrm{0.6366}\cdot \text{arctan}\left[\mathrm{0.8}{T}_{\mathrm{\Sigma }}-\mathrm{213.6}\mathrm{K}\right],\text{(B4)}& {\displaystyle }{k}_{\mathit{\psi }}=\sin \left(\mathit{\psi }\right)-\tan \left({\mathit{\theta }}_b\right)\cdot \cos \left(\mathit{\psi }\right),\text{(B5)}& {\displaystyle }\mathit{ϵ}=f\left(T_{\mathrm{\Sigma }}\right)g\left(r\right)=\left[e^{\frac{-(T_{\mathrm{\Sigma }}-\mathrm{266}\mathrm{K}{)}^{\mathrm{2}}}{\mathrm{10}}}\right]\left[{\displaystyle \frac{r^{\mathrm{2}}{e}^{\frac{-r^{\mathrm{2}}}{\mathrm{2}{a}^{\mathrm{2}}}}}{a^{\mathrm{3}}}}\right],\end{array}$$

with $r=\frac{h_{\mathrm{\Phi }}}{h_{\mathrm{\Sigma }}}$ and a=1.1.

Air into cloud

$$\begin{array}{}\text{(B6)}& {\dot{H}}_{\mathrm{\Lambda }\to \mathrm{\Pi }}=\left(\mathrm{1.16}\mathit{\psi }+\mathrm{0.013}\sqrt{R_{\mathrm{\Pi }}{\widehat{h}}_{\mathrm{\Pi }}}\right)\left({\mathit{\rho }}_{\mathrm{\Pi }}-{\mathit{\rho }}_{\mathrm{\Lambda }}\right)\end{array}$$

Powder cloud drag

$$\begin{array}{}\text{(B7)}& {\displaystyle }{S}_{\mathrm{\Pi }}=\mathrm{0.1}\left[{\mathit{\mu }}_{\mathrm{lam}}{u}_{\mathrm{\Pi }}+\mathrm{0.08}{\mathit{\mu }}_{\mathrm{turb}}{R}_{\mathrm{\Pi }}{\widehat{h}}_{\mathrm{\Pi }}\right],\text{(B8)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{lam}}=-\mathrm{0.22}+\mathrm{0.28}{a}_{\mathrm{drag}},\text{(B9)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{turb}}=\mathrm{0.02}+\mathrm{0.035}{a}_{\mathrm{drag}}.\end{array}$$

The range of powder cloud drag depends on mean speed and total turbulent energy. A lower decay number, β_Π, results in higher drag and a narrower cloud flow, while a higher β_Π leads to lower drag and wider flow.

Dispersive pressure

Dispersive pressure (Bagnold, 1954) arises in dense snow avalanches and other dry, granular-type debris flows due to shear deformation in granular material (Buser and Bartelt, 2011).

It causes the core height h_z to increase and shifts the centre of mass k_z:

$$\begin{array}{}\text{(B10)}& k_z={\displaystyle \frac{h_z}{\mathrm{2}}}.\end{array}$$

RAMMS::EXTENDED models this with the following equations:

$$\begin{array}{}\text{(B11)}& {\displaystyle }\mathcal{D}(t,h_z,{\dot{h}}_z,{\ddot{h}}_z)={\dot{h}}_z(h_z,{\dot{h}}_z,{\ddot{h}}_z),\text{(B12)}& {\displaystyle }{\left(\dot{{\displaystyle \frac{h_z}{\mathrm{2}}}}\right)}_t+\text{div}\left(\dot{{\displaystyle \frac{h_z}{\mathrm{2}}}}{\mathit{u}}_{\mathrm{\Phi }}\right)=\ddot{{\displaystyle \frac{h_z}{\mathrm{2}}}},\text{(B13)}& {\displaystyle }{\left(\ddot{{\displaystyle \frac{h_z}{\mathrm{2}}}}\right)}_t+\text{div}\left(\ddot{{\displaystyle \frac{h_z}{\mathrm{2}}}}{\mathit{u}}_{\mathrm{\Phi }}\right)={\displaystyle \frac{{\dot{W}}_z}{H_{\mathrm{\Phi }}}}-\left[g_z+\ddot{{\displaystyle \frac{h_z}{\mathrm{2}}}}\right]{\displaystyle \frac{{\dot{h}}_z/\mathrm{2}}{h_z/\mathrm{2}}}.\end{array}$$

Here, u_Φ is the average core velocity, ${\dot{W}}_z$ is the shearing work rate in the z direction and M_Φ is the core mass. The equations are solved using a Harten–Lax–van Leer (HLL) finite volume scheme (Christen et al., 2010).