Numerical modelling of dense snow avalanches with a well-balanced scheme based on the 2D shallow water equations

Abstract

A common technique for simulating non–Newtonian fluid dynamics, such as snow avalanches, is to solve the Shallow Water Equations (SWE), together with a rheological model describing the momentum dissipation by shear stresses. Friction and cohesion terms are commonly modelled using the Voellmy friction model and, recently, the Bartelt cohesion model. Here, an adaptation of the Roe scheme that ensures the balance between the flux and pressure gradients and the friction source term is presented. An upwind scheme was used for the discretisation of the SWE numerical fluxes and the non–velocity-dependent terms of the friction–cohesion model, whereas a centred scheme was used for the velocity-dependent source terms. The model was tested in analytically solvable settings, laboratory experiments and real cases. In all cases, the model performed well, avoiding numerical instabilities and achieving stable and consistent solution even for an avalanche stopping on a sloping terrain.

1. Introduction

The growing concerns regarding natural hazards, particularly snow avalanches, have led to development of numerical models as support tools for the analysis of these hazards (Dent and Lang, 1980; Perla and others, 1980; Margreth and Funk, 1999; Blagovechshenskiy and others, 2002; Jamieson and others, 2008; Bartelt and others, 2015; Fischer and others, 2015; Eglit and others, 2020). From the viewpoint of their dynamics, avalanches have traditionally been classified into dense flow avalanches and powder snow avalanches (Gauer and others, 2008). However, large avalanches are usually classified as mixed-motion avalanches, which include a flowing, dense core and a less dense, fluidised transition zone that develops ahead and over the core, along with a very diluted, suspension cloud or aerosol (Sovilla and others, 2015; Issler and others, 2019). In some cases, the dense core is maintained by cohesion forces that provide additional flow resistance, thereby jointly retaining the snow particles. These facts highlight the complexity and challenges associated with numerical modelling of this physical process.

Up to date, most numerical simulation tools available for practitioners provide information about runout distance, flow depths, flow velocities and impact pressure for dense flow snow avalanches or for the dense core of mixed-motion avalanches. The first models were developed in the 1920s, which evolved into mass-point models (Ancey and others, 2005). Afterwards, 1D models with simple geometries, based on cross-sections or longitudinal profiles, were developed (Christen and others, 2001; Podolskiy and others, 2013; Eglit and others, 2020). Later came 2D models (Christen and others, 2010; Eglit and others, 2020), which are widely used today. Typically, 2D-models achieve a better discretisation of the spatial domain by using a mesh of elements, thereby improving the representation of the avalanche dynamics and the deposit distribution. Nowadays, the 3D smoothed-particle hydrodynamics (SPH) technique, where the particles are spheres (López and others, 2010; Salazar and others, 2016), is used in very few numerical models (Schraml and others, 2015), and mainly for research purposes due to the high computational cost. For the dense flow layer, Sampl and Granig (2009) use a 2D-SPH method that assimilates each particle to a cylinder with the height of the avalanche depth.

In both 1D and 2D models, the simulation of dense snow avalanches is based on the solution of the depth averaged mass and momentum conservation equations, which are also used for the simulation of free surface water flows. In the Shallow Water Equations (SWE) for water flows, the frictional terms are usually represented with the Manning formula; however, for any non–Newtonian shallow flow, such as snow avalanches, other specific rheological models are required.

The Voellmy friction model (Voellmy, 1955) is most often used to define the friction terms for granular flows (Pirulli and Sorbino, 2008; Schraml and others, 2015), but other approaches have also been proposed (Issler and Gauer, 2008; Faccanoni and Mangeney, 2013; Issler and others, 2018). This friction model expresses the total flow resistance as the sum of two contributions, namely turbulent friction resistance and dry-Coulomb friction, respectively. To help avalanche modellers with the selection of the appropriate values of these parameters, different guidelines and handbooks have been published (Buser and Frutiger, 1980; Bakkehøi and others, 1981; Brugnot, 2000; Bartelt and others, 2017). Despite that, the choice of the Voellmy parameters in practical applications has been the object of many analyses and discussions (Bartelt and others, 1999; Ancey and others, 2004; Gruber and Bartelt, 2007; Issler and Gauer, 2008; Keylock and Barbolini, 2011; Gauer, 2014; Fischer and others, 2015; Issler and others, 2018).

When solving the SWE, or any hyperbolic system of equations with source terms in general, the numerical scheme must achieve a proper balance between the homogeneous part of the equations and the source term (Bladé and others, 2008, 2012a); otherwise, oscillations leading to instabilities may occur. Alternative, for some schemes, the imbalance can make them incapable to simulate quiescent states (Capart and others, 2003; Hou and others, 2013). In the case of snow avalanches, the source term contains the friction model, and an improper balancing of it may cause the avalanche to not stop flowing even for very small slopes. For such cases, Bartelt and others (2017) proposed a stopping criterion based on controlling the momentum, where the avalanche is made to stop when its momentum is lower than a user-defined fraction of its maximum momentum. However, this criterion lacks a physical basis, as the maximum momentum depends on the avalanche's characteristics at very different location and time than those when it stops. To address this issue, Bartelt and others (2015) proposed the inclusion of an additional friction term related to snow cohesion, which is a real physical snow property that has the effect of retention and can stop the avalanche irrespective of the maximum momentum that is reached during the avalanche propagation. With this cohesion term, the stopping position crucially depends on the flow depth. Nevertheless, the fraction of momentum method introduced by Bartelt and others (2017) continues to be an extensively used method not only for practitioners but also for researchers (Schraml and others, 2015; Wever and others, 2018).

Within this framework, this study presents a specific numerical treatment for the non–velocity-dependent friction terms based on an upwind discretisation, which are counterbalanced with the pressure forces. The scheme ascertains that the avalanche stops with no need for any additional condition and, thus, the balance of the different terms of SWE is ensured. The developments presented in this study are based on an adaptation of the numerical scheme used in Iber (Bladé and others, 2014), which is a 2D hydraulic numerical modelling tool that uses the Roe scheme (Roe, 1986) and couples Godunov's method with Roe's approximate Riemann solver (Toro, 2009). The scheme was applied in equilibrium and quiescent states, the numerical results were compared with the laboratory experiments and the case study of the 2014 snow avalanche of Bonaigua. Additionally, the effect of considering nonhydrostatic pressure or anisotropic flow properties was analysed upon two practical cases, and demonstrated a good performance for all of them in terms of the avalanche dynamics and the deposit distribution.

2. Governing equations

2.1. 2d Shallow water equations for non–Newtonian flows

Most of the existing avalanche simulation tools are based on the solution of mass and momentum conservation equations, which are written similarly to the equations for free surface water flows but differ in the terms describing friction (i.e. the rheological model). These equations, when applied to water under a 2D framework, are named 2D Shallow Water Equations (2D-SWE). These equations are derived from the Navier–Stokes equations, through a time averaging to filter the turbulent fluctuations (Reynolds Averaged Navier–Stokes equations, or RANS) and a depth averaging to obtain the final 2D equations.

The 2D-SWE are a hyperbolic nonlinear system of three partial differential equations, which can be written in compact vectorial notation as follows:

(1)$$\displaystyle{\partial \over {\partial t}}{\boldsymbol U} + \displaystyle{\partial \over {\partial x}}{\boldsymbol F}( {\boldsymbol U} ) + \displaystyle{\partial \over {\partial y}}{\boldsymbol G}( {\boldsymbol U} ) = {\boldsymbol H}( {\boldsymbol U} ) $$

where U is the vector of the conservative variables, F and G are the x and y components of the flux vectors, respectively, and H is the source term. The momentum equations contain the gradients of the pressure and inertia terms (through the flow vectors F and G) as well as the bottom slope and friction terms (through the source term H).

While using 2D-SWE-based numerical models to simulate non–Newtonian shallow flows in a global coordinate system, particularly for dense snow avalanches, one must modify the SWEs because the bottom slope is usually not small and the cosine of the angle (θ) composed by the bed normal and the vertical axis cannot be approximated by 1. This has a direct effect not only on the pressure terms (Chow, 1959) but also on the bed friction and the bed slope.

Considering the 1D problem presented in Figure 1, the friction forces exerted over the bed and the pressure terms can be corrected by replacing the gravity acceleration g by g′ = gcos ²θ (Ni and others, 2019; Zugliani and Rosatti, 2021; Maranzoni and Tomirotti, 2022). Thus, assuming that the free surface flow is parallel to the terrain, hydrostatic pressure in the local axis z′ is produced such that $p{\rm ^{\prime}} = \rho g\cos \theta ( {h{\rm^{\prime}}-z{\rm^{\prime}}} )$. By accounting for the relation among the coordinate systems, where z′ = zcosθ and h′ = hcosθ, the pressure distribution can be expressed as p = ρgcos ²θ(h − z). Moreover, the velocity in the normal direction can be considered negligible (Savage and Hutter, 1989) and, thus, the velocity vector lies on a plane parallel to the terrain. Therefore, the depth-averaged velocity follows the relation v _x = v′_xcosθ in the global coordinate system with the same norm as for the local coordinate system (Gray and others, 1999; Iverson and Denlinger, 2001).

Figure 1. One-dimensional problem described by using a global coordinate system and a local coordinate system.

Furthermore, two hypotheses are usually considered: a monophasic fluid, in which the fluid is formed by a unique phase where all components are perfectly mixed, and shear stress grouping, in which the effect of different shear stresses are grouped as five components of a single term (Julien and León, 2000) as follows:

(2)$$\tau = \tau _d + \tau _t + \tau _v + \tau _{mc} + \tau _c$$

where τ _d represents the dispersive term, τ _t the turbulent term, τ _v the viscous term, τ _mc the Mohr–Coulomb terms and τ _c the cohesive term. In these components, the appropriate rheological model for the particular purpose of each work is obtained by selecting one or several components of Eqn (2).

For water flow, the shear terms due to friction are generally incorporated in the equations through of the concept of friction slope. It is part of the source term H of the momentum equations and expresses the contribution of the momentum change caused by the energy dissipation resulting from flow-boundary interactions and sub-grid turbulence (if, as is generally done, no turbulence model is used). By analogy, the rheological model for non–Newtonian fluids is sometimes expressed as a friction slope (S _rh), which comprises one or more of the aforementioned components of the total shear stress in Eqn (2).

As for free surface water flows, the hydrostatic and isotropic pressure distribution is assumed (Chow, 1959). This means that there is a linear variation of pressure, or a pressure gradient of the specific weight value of the fluid, in the direction normal to streamlines. For non–Newtonian flows, and depending on the nature of the fluid, this assumption may not be realistic and, in turn, lead to inaccurate results. In particular, the consideration of nonhydrostatic and anisotropic pressure distribution can help improve the representation of the avalanche dynamics (Hungr, 1995; Bartelt and others, 1999; Hungr and McDougall, 2009; Ruiz-Villanueva and others, 2019). Moreover, it can be especially relevant during the first moments after a snow avalanche release. This pressure correction is commonly made through of a factor K _p (Savage and Hutter, 1989), which multiplies the pressure terms in the momentum equations. A K _p value equal to 1 implies hydrostatic pressure distribution. Thus, the terms of Eqn (1), for a non–Newtonian shallow flow (e.g. a snow avalanche) can be written as follows:

(3)$$\matrix{ {\matrix{ {{\boldsymbol U} = \left[{\matrix{ h \cr {hv_x} \cr {hv_y} \cr } } \right]} & {{\boldsymbol F} = \left[{\matrix{ {hv_x} \cr {hv_x^2 + g{\rm^{\prime}}\displaystyle{{h^2} \over 2}K_{\,p, x}} \cr {hv_xv_y} \cr } } \right]} \cr } } & \cr{\matrix{ {{\boldsymbol G} = \left[{\matrix{ {hv_y} \cr {hv_xv_y} \cr {hv_x^2 + g{\rm^{\prime}}\displaystyle{{h^2} \over 2}K_{\,p, y}} \cr } } \right]} & {{\boldsymbol H} = {\boldsymbol \;}\left[{\matrix{ 0 \cr {g{\rm^{\prime}}h( {S_{o, x}-S_{rh, x}} ) } \cr {g{\rm^{\prime}}h( {S_{o, y}-S_{rh, y}} ) } \cr } } \right]} \cr } } \cr } $$

where h is the flow depth, v _x and v _y are the two depth-averaged velocity components, g′ is the projected gravitational acceleration, S _o,x and S _o,y are the two bottom slope components, and S _rh,x and S _rh,y are the two friction slope components of the rheological model. The bottom slope is computed as S_o = (∂z _b/∂x, ∂z _b/∂y)^T, where z _b is the bed elevation in the global coordinates.

Several rheological models for non–Newtonian fluid flows have been proposed in the literature. These models usually consist of two terms: one related to the yield stresses and the other to turbulent dissipation. In the case of snow avalanches, the Voellmy friction model is probably the most common approach. Therefore, it is the approach taken here as a reference, along with the Bartelt and others (2015) extension to include cohesion. Nevertheless, the developments that follow can be extended to other rheological models with the same yield-turbulent stress structure.

For the friction–cohesion (Voellmy–Bartelt) model, S_rh can be split into two as follows:

(4)$${\boldsymbol S}_{rh} = {\boldsymbol S}_{rh}^{\rm ^{\prime}} + {\boldsymbol S}_{rh}^{\rm \prime\prime } $$

where ${\boldsymbol S}_{rh}^{\rm ^{\prime}}$ is the friction slope contribution for the flow resistance forces while ${\boldsymbol S}_{rh}^{\rm \prime\prime }$ constitutes the cohesion forces.

The Voellmy friction model integrates the total flow resistance $( {\boldsymbol S}_{rh}^{\rm ^{\prime}} )$ as the sum of a solid phase, which is related to the Coulomb friction coefficient (μ), and a turbulent resistance coefficient (ξ):

(5)$${\boldsymbol S}_{rh}^{\rm ^{\prime}} = ( {S_{rh, x}^{\rm^{\prime}} , \;S_{rh, y}^{\rm^{\prime}} } ) = \left({\mu + \displaystyle{{v_x\sqrt {v_x^2 + v_y^2 + v_z^2 } } \over {\xi h{\cos }^2\theta }}, \;\;\mu + \displaystyle{{v_y\sqrt {v_x^2 + v_y^2 + v_z^2 } } \over {\xi h{\cos }^2\theta }}} \right)$$

where v _z = v _xtanθ _x + v _ytanθ _y according to Zugliani and Rosatti (2021).

According to Bartelt and others (2015), the cohesion forces $( {\boldsymbol S}_{rh}^{\rm \prime\prime } )$ can be defined as an additional flow resistance that depends on the cohesion (C) and the normal stress:

(6)$$\eqalignb{{\boldsymbol S}_{rh}^{\rm \prime\prime } & = ( {S_{rh, x}^{\rm \prime\prime } , \;S_{rh, y}^{\rm \prime\prime } } ) \cr & = \left({\displaystyle{1 \over {\rho g{\rm^{\prime}}h}}C( {1-\mu } ) \left({1-e^{-{{\rho {g}^{\prime}h} \over C}}} \right)\cos \theta_x, \;\;}\right. \cr & \qquad \left. \displaystyle{1 \over {\rho g{\rm^{\prime}}h}}C( {1-\mu } ) \left({1-e^{-{{\rho {g}^{\prime}h} \over C}}} \right)\cos \theta_y \right)}$$

where ρ is the flow density, ξ is the turbulent friction coefficient, μ is the Coulomb friction coefficient, C is the cohesion parameter, and θ _x and θ _x the two components of the angle composed by the bed normal and the vertical axis.

This rheological model has been implemented in a numerical model based on Iber. It is a 2D numerical tool that was initially developed for modelling hydrodynamic and sediment transport (Bladé and others, 2014, 2019) and solves the SWE on irregular geometries by using the Finite Volume Method (FVM). The tool has been continuously enhanced and, currently, includes a series of modules for different free surface flow processes, such as hydrological processes (Cea and Bladé, 2015; Sanz-Ramos and others, 2022), pollutant propagation (Cea and others, 2016), large-wood transport (Ruiz-Villanueva and others, 2014) and physical habitat suitability assessment (Sanz-Ramos and others, 2019, 2023). As for some other existing modules of Iber (Ruiz-Villanueva and others, 2014; Cea and Bladé, 2015; Cea and others, 2016), the computation of dense-snow avalanche dynamics requires a specific adaptation of the numerical scheme.

2.2. Numerical scheme

The final discretisation of the equations depends on the numerical method (i.e. finite differences, finite elements, finite volumes, etc.) as well as the particular numerical scheme and its characteristics (for fluid dynamics, either the centred or upwind scheme). Centred schemes use an approximation of the variables in each calculation point obtained by interpolating the surrounding point values with no consideration of the flow direction. By contrast, upwind schemes can be spatially biased and based on the characteristics theory for the hyperbolic systems of equations (Lax, 1973) and, thus, reproduce the fact that flow perturbations propagate along the characteristic lines in space and time. This results in the waves propagating in different directions and velocities in relation to the fluid velocity, which depends on the slope of the characteristic lines and, in their turn, on the Froude number (Fig. 2). The differences in the state variables between the neighbouring cells lead to the propagation of waves with characteristic speeds. Hence, an increase of flow depth in the downhill direction may lead to a wave propagating uphill in the numerical model, even if the terrain is steeper than flow the depth increases and the avalanche surface coordinate decreases in the downhill direction (Fig. 2).

Figure 2. A sketch of the wave propagation directions generated by a perturbation of the flow depending on the flow regime: (a) fluid at rest over horizontal terrain, (b) fluid in movement with subcritical flow (F _R < 1), and (c) fluid in movement with supercritical flow (F _R > 1). F _R is the Froude number, which is the quotient between the inertial and gravitational forces ($v/\sqrt {gh}$ for free surface flows, where v is the velocity, c is the celerity, h is the depth, and g is the gravitational acceleration).

Upwind schemes consider the flow velocity, the wave propagation direction, and the celerity in the discretisation of the differential equations. To achieve a well-balanced numerical scheme, the same upwind scheme has to be applied to both the homogeneous part of the equations as well as the source term, which is not trivial in some cases.

The proposed numerical model solves the 2D-SWE by using a conservative scheme based on the FVM on a structured mesh, unstructured mesh, and/or a combination of both kinds of meshes composed of triangles and/or quadrilaterals. For the convective fluxes, it uses an explicit first-order Godunov-type upwind scheme (Toro, 2009), in particular the Roe scheme (Roe, 1986) (under which name the Godunov Method with the approximate Riemann solver of Roe is commonly known), with an upwind discretisation for the geometric slope source terms (Vázquez-Cendón, 1999). A centred scheme is used for the other source terms.

When the 2D-SWE are extended to model non–Newtonian shallow flows, such as dense snow avalanches, the source term H is modified with the additional terms previously indicated in Eqn (3). As for the original SWE, the proper balance between the source term H and the gradients of the flux vectors F and G must be ensured (Bladé and others, 2012a). For snow avalanche modelling, some of the friction terms in H (i.e. those of the solid phase and cohesion) do not depend on the velocity. This is in contrast to what happens with water, where all friction terms vanish for the quiescent flow. This means that these terms must be properly balanced with the pressure terms at the left side of the equation, within the flux vector.

The FVM calculates the average values of the flow variables (depth and velocity) for each finite volume. However, these averages are updated at each time step with the flows (mass flow or discharge and momentum flow) across each element side. In the Roe scheme, the decomposition of the integral of the flow vectors on a finite volume, which results in the evaluation of the fluxes across the finite volume surface, is performed as a non–centred linear combination of the eigenvectors of the Jacobian matrix (J_n) of the flow vectors as follows:

(7)$${\boldsymbol J}_{\boldsymbol n} = \left({\matrix{ 0 & {n_x} & {n_y} \cr {( {{g}^{\prime}h-v_x} ) n_x-v_xv_yn_y} & {v_yn_y + 2v_xn_x} & {v_xn_y} \cr {( {{g}^{\prime}h-v_y} ) n_y-v_xv_yn_x} & {v_yn_x} & {v_xn_x + 2v_yn_y} \cr } } \right)$$

The discretisation of the proposed scheme is applied to mesh elements of any number of faces, but it might be equally applied to regular grids. A sketch of two neighbour cells and the geometric variables that will be used in the discretisation is presented in Figure 3. The geometric centre of the control volumes stores all the flow variables.

Figure 3. The geometric variables used to compute the flux between the elements of an unstructured finite-volume mesh composed of quadrilateral and triangle elements. V _i and V _j are the areas of the control volumes of the elements i and j, respectively, ${\boldsymbol n}_{i, w_l}$ is the exterior normal vector on the element edge w _l, and $l_{i, w_l}$ is the length of that edge (w _l).

According to Figure 3, any explicit finite volume scheme for the 2D-SWE, applied to an element (i.e. finite volume), can be expressed as follows:

(8)$${\boldsymbol U}_i^{n + 1} = {\boldsymbol U}_i^n -\displaystyle{{\Delta t} \over {V_i}}\left[{\mathop \sum \limits_{l = 1}^{N_i} ( {{\boldsymbol F}_{i, w_l}^\ast {\boldsymbol n}_{i, w_l}} ) l_{i, w_l}} \right] + \Delta t{\boldsymbol H}_i^{} $$

where Δt is the time step, V _i is the area of the control volume, ${\boldsymbol n}_{i, w_l}$ is the exterior normal vector of the element edge w _l, $l_{i, w_l}$ is the length of that edge (w _l), N _i is the number of edges of that element, F^* are the numerical fluxes, and ${\boldsymbol H}_i^{}$ is the integral of the source term over the element i, which includes friction and slope.

As Eqn (8) is explicit in time, the stability criterion of Courant–Friedrichs–Lewy (CFL) over the computational time step applies (Courant and others, 1967). The CFL condition establishes a relation between the flow velocity and the water depth, the element size and the maximum permissible computational time step. The presented scheme uses the computational time step proposed by Cea and Bladé (2015).

As stated before, to ensure the equilibrium between the pressure terms of the flow vector and the bottom slope (in the source term), which is crucial to avoid spurious oscillations or an unphysical movement of the free surface for quiescent water on irregular geometries, a decomposition identical to that of the pressure terms has to be performed for the bottom slope (Bermúdez and others, 1998; Vázquez-Cendón, 1999; LeVeque, 2002; Toro, 2009; Bladé and others, 2012a, 2012b). Therefore, an upwind discretisation for the slope source terms (S_o) must also be used. This results in the source term H in Eqn (3) being separated in its slope (H¹) and friction (H²) contributions as follows:

(9)$${\boldsymbol H} = {\boldsymbol H}^1 + {\boldsymbol H}^2$$

where

(10)$${\boldsymbol H}^1 = {\boldsymbol \;}\left[{\matrix{ 0 \hfill \cr {g{\rm^{\prime}}hS_{o, x}} \hfill \cr {g{\rm^{\prime}}hS_{o, y}} \hfill \cr } } \right]; \;\;{\boldsymbol H}^2 = \left[{\matrix{ 0 \hfill \cr {-g{\rm^{\prime}}hS_{rh, x}} \hfill \cr {-g{\rm^{\prime}}hS_{rh, y}} \hfill \cr } } \right]$$

Instead of a continuous sloping surface, the part corresponding to the terrain slope (H¹) is discretised through a series of ‘steps’ (Δz)_i,j between adjacent cells i and j that separate horizontal surfaces (Cea and Vázquez-Cendón, 2012). The notation $\widetilde{{}}$ used in the following equations represents the Roe averages of the variables at the element edges described below, H¹ being:

(11)$${\boldsymbol H}^1 = \left[{\matrix{ 0 \hfill \cr {-g{\rm^{\prime}}h{( {\Delta z} ) }_{i, j}n_x} \hfill \cr {-g{\rm^{\prime}}h{( {\Delta z} ) }_{i, j}n_y} \hfill \cr } } \right]$$

where n _x and n _y are the x and y components of the outer normal vector to the edge i,j, respectively, and (Δz)_i,j is the difference in the terrain elevation between adjacent elements.

Thus, Eqn (8) can be expressed as follows:

(12)$$\Delta t{\boldsymbol H}_i^{} = {-}\displaystyle{{\Delta t} \over {V_i}}\left[{\mathop \sum \limits_{l = 1}^{N_i} ( {{\boldsymbol H}_{i, w_l}^{{\ast} 1} } ) } \right] + \Delta t{\boldsymbol H}_i^2 $$

In this equation, the term containing the bottom slope, ${\boldsymbol H}_{i, w_l}^{{\ast} 1}$, is finally expressed as follows:

(13)$${\boldsymbol H}_{i, w_l}^{{\ast} 1} = \displaystyle{1 \over 2}l_{i, w_l}\left({\mathop \sum \limits_{k = 1}^3 {\tilde{\beta }}_k( {1-{\rm sign}( {{\tilde{\lambda }}_k} ) {\tilde{{\boldsymbol e}}}_k} ) } \right)_{i, j}$$

where $\tilde{\beta }_k$ are the coefficients that allow the decomposition of the slope source term as a linear combination of the eigenvectors of the Jacobian matrix in Eqn (7). Moreover, $\tilde{{\boldsymbol e}}_k$ and $\tilde{\lambda }_k$ are the eigenvectors and eigenvalues of the Jacobian matrix of the Roe approximate flux vector, which are as follows:

(14)$$\matrix{ {{\tilde{\lambda }}_1 = {\tilde{v}}_xn_x + {\tilde{v}}_yn_y + \tilde{c}} \cr {{\tilde{\lambda }}_2 = {\tilde{v}}_xn_x + {\tilde{v}}_yn_y} \cr {{\tilde{\lambda }}_3 = {\tilde{v}}_xn_x + {\tilde{v}}_yn_y-\tilde{c}} \cr } $$

(15)$$\tilde{e}_1 = \left[{\matrix{ 1 \hfill \cr {{\tilde{v}}_x + \tilde{c}n_x} \hfill \cr {{\tilde{v}}_y + \tilde{c}n_y} \hfill \cr } } \right]\;; \;\;\tilde{e}_2 = \left[{\matrix{ 0 \hfill \cr {-\tilde{c}n_y} \hfill \cr {\tilde{c}n_x} \hfill \cr } } \right]\;; \;\;\tilde{e}_3 = \left[{\matrix{ 1 \hfill \cr {{\tilde{v}}_x-\tilde{c}n_x} \hfill \cr {{\tilde{v}}_y-\tilde{c}n_y} \hfill \cr } } \right]$$

where n _x and n _y are the x and y components of the outer normal vector to the edge i,j, respectively, and c is the wave celerity. Moreover, the terms $\tilde{u}$, $\tilde{v}$ and $\tilde{c}$ are the Roe averaged states of the conserved variables (Toro, 2009). These are as follows:

(16)$$\tilde{u} = \displaystyle{{\sqrt {h_i} u_i + \sqrt {h_j} u_j} \over {\sqrt {h_i} + \sqrt {h_j} }}; \;\;\tilde{v} = \displaystyle{{\sqrt {h_i} v_i + \sqrt {h_j} v_j} \over {\sqrt {h_i} + \sqrt {h_j} }}; \;\;\tilde{c} = \sqrt {g^{\prime}\displaystyle{{h_i + h_j} \over 2}} $$

Upon the definition of the slope part of the source term, in Eqn (13) the expression for the $\tilde{\beta }_k$ coefficients must be the following:

(17)$$\eqalign{{\tilde{\beta }}_1 = & -\displaystyle{{\tilde{c}} \over 2}( {\Delta z} ) _{i, j} \cr {\tilde{\beta }}_2 = & 0 \cr {\tilde{\beta }}_3 = & \displaystyle{{\tilde{c}} \over 2}( {\Delta z} ) _{i, j}} $$

when modelling shallow flows of non–Newtonian fluids, S_rh represents the friction–cohesion model, which may be based on the Voellmy–Bartelt model in the case of snow avalanches. For its turbulent-drag part that is proportional to 1/ξ, which is velocity-dependent, a centred scheme can be used, as is done with the drag in water flows (Cea and Vázquez-Cendón, 2012; Bladé and others, 2012b). Conversely, in the flux vector discretisation, the pressure terms must be properly balanced with the Coulomb stresses (μ-dependent) and the cohesion stresses (C-dependent). Thus, these two friction terms must be treated with an upwind discretisation as is done with the bottom slope. The physical interpretation of this is that the three terms (i.e. bottom slope, Coulomb stresses and cohesion) represent the forces that are counterbalanced with the pressure forces for a flow at rest. All this results in the following scheme for the friction term H²:

(18)$$\Delta t{\boldsymbol H}_i^2 = {-}\displaystyle{{\Delta t} \over {V_i}}\left[{\mathop \sum \limits_{l = 1}^{N_i} ( {{\boldsymbol H}_{i, w_l}^{{\ast} 2( {\mu + C} ) } } ) } \right] + \Delta t{\boldsymbol H}_i^{2\xi } $$

where the ${\boldsymbol H}_{i, w_l}^{{\ast} 2( {\mu + C} ) }$ includes the Coulomb and cohesion stresses, which leads to the following final expression for the scheme:

(19)$${\boldsymbol U}_i^{n + 1} = {\boldsymbol U}_i^n -\displaystyle{{\Delta t} \over {V_i}}\left[{\mathop \sum \limits_{l = 1}^{N_i} ( {{\boldsymbol F}_{i, w_l}^\ast {\boldsymbol n}_{i, w_l}} ) l_{i, w_l} + \mathop \sum \limits_{l = 1}^{N_i} ( {{\boldsymbol H}_{i, w_l}^{{\ast} 1} + {\boldsymbol H}_{i, w_l}^{{\ast} 2( {\mu + C} ) } } ) } \right] + \Delta t{\boldsymbol H}_i^{2\xi } $$

In Eqn (19), the term containing the slope, ${\boldsymbol H}_{i, w_l}^{{\ast} 1}$, is discretised for water in the same way as for other rheological models with non–velocity-dependent terms, such as the Voellmy–Bartelt model used here. By analogy with the bottom slope treatment, the term containing the Coulomb stress and cohesion $( {\boldsymbol H}_{i, w_l}^{{\ast} 2( {\mu + C} ) } )$ is discretised in a similar manner as Eqn (13). This is expressed as follows:

(20)$${\boldsymbol H}_{i, w_l}^{{\ast} 2( {\mu + C} ) } = \displaystyle{1 \over 2}l_{i, w_l}\left({\mathop \sum \limits_{k = 1}^3 {\tilde{\beta }}_{( {{\boldsymbol \mu } + {\boldsymbol C}} ) }( {1-{\rm sign}( {{\tilde{\lambda }}_k} ) {\tilde{{\boldsymbol e}}}_k} ) } \right)_{i, j}$$

where $\tilde{\lambda }_k$ and $\tilde{{\boldsymbol e}}_k$ are the expressions shown in Eqns (14) and (15) respectively, and the $\tilde{\beta }_{( {{\boldsymbol \mu } + {\boldsymbol C}} ) }$ coefficients are as follows:

(21)$$\eqalign{{\tilde{\beta }}_{( {\mu + C} ) 1} = & -\displaystyle{{\tilde{c}} \over 2}( {\Delta z_{( {\mu + C} ) }} ) _{i, j} \cr {\tilde{\beta }}_{( {\mu + C} ) 2} = & 0 \cr {\tilde{\beta }}_{( {\mu + C} ) 3} = & \displaystyle{{\tilde{c}} \over 2}( {\Delta z_{( {\mu + C} ) }} ) _{i, j}} $$

Now, the term Δz _(μ+C) represents a fictitious bottom ‘step’, or friction step that would have the same effect on the flow motion as the actual Coulomb and cohesion stresses. It can be expressed in the following manner for the non–velocity-dependent part of the Voellmy–Bartelt friction–cohesion model:

(22)$$\Delta z_{( {{\boldsymbol \mu } + {\boldsymbol C}} ) } = \left({\tilde{\mu } + \displaystyle{1 \over {\rho g^{\prime}\tilde{h}}}\tilde{C}( {1-\tilde{\mu }} ) ( {1-e^{-( \rho {g}^{\prime}\tilde{h}/\check{C}) }} ) } \right)d_{i, j}$$

where $\tilde{\mu }\;$and $\tilde{C}$ are the arithmetic averages of μ and C at the elements i and j, and d _i,j is the distance between the centre of gravity of the two elements. The contribution of Δz _(μ+C) is limited to counterbalancing the pressure terms in a static situation; otherwise is limited to the momentum. It should be noted that if the cohesion contribution is neglected, this discretisation is also valid, and only the contribution of the Coulomb part will be considered.

Conceptually, this numerical treatment can be interpreted as the addition of an artificial elevation difference (friction step) to the actual elevation difference between the adjacent elements (geometric step) (Fig. 4). With this, the equilibrium of a stopped avalanche on a sloping terrain due to the balance of gravity forces and friction is properly represented. If the avalanche is still moving, the friction step Δz _(μ+C) acts against the fluid motion, whereas when it has stopped, the friction step Δz _(μ+C) completely counterbalances the geometric slope.

Figure 4. A sketch of the numerical treatment of Coulomb (μ) and cohesion (C) friction stresses as a ‘friction step’: (a) for a stopped avalanche, the geometric step (Δz) (upper left) is counterbalanced by the ‘friction step’ (Δz _(μ+C)) (upper middle), and thus the velocity is kept null (upper right) and (b) the friction step (lower middle) opposes the avalanche flow, the calculated velocity (lower right) is lower than the velocity in the case of no friction, and the avalanche keeps moving (lower figure left); moreover, h _i and h _j are the flow depth at the elements i and j respectively, and v is the flow velocity.

2.3. Wet-dry fronts

The interfaces between the elements belonging to the avalanche itself and the adjacent elements must also be treated properly to conserve mass. In the solution of the SWE, one commonly defines a fluid-depth threshold (ɛ_wd) below which a finite volume (i.e. mesh element) is considered to be dry. Proper numerical treatment is mandatory for avoiding numerical instabilities, especially for irregular geometries (Cea and others, 2007; Liang and Borthwick, 2009; Herty and Seaid, 2018).

In 2D-SWE-based modelling of non–Newtonian fluids as well as for water flows, the ɛ_wd parameter must be a value that allows for the proper representation of the extent of the area occupied by the fluid. However, in contrast to free surface water flows (some examples can be found in references (Cea and others, 2007; Bladé and others, 2012b; Cea and Bladé, 2015; Herty and Seaid, 2018; Sanz-Ramos and others, 2020)), no references have been found in the existing literature regarding the ɛ_wd values that are commonly used for dense snow modelling.

The numerical scheme presented here integrates different algorithms for the wet-dry front treatment based on the study of Cea and others (2007). In brief, to ensure mass conservation, the finite volume is considered to be dry when the fluid depth (h) in an element is lower than the threshold (h < ɛ_wd). In wet-dry fronts, the mass and momentum fluxes are controlled, i.e. mass and momentum fluxes are allowed from wet to dry elements; however, in case the Roe scheme produces fluxes from dry to wet elements, which is possible with the Roe averages, the fluxes are all set to zero. With this method, it is possible to ensure mass conservation without errors due to the discretisation of the wet-dry fronts.

3. Applications and discussion

This section presents the performance of the numerical numerical tool by simulating three idealised flow processes, a laboratory experiment, and a real snow avalanche event. Initially, to test the robustness and well-balancedness of the scheme, two classical 2D dam-break problems and a straight, inclined channel with an obstacle in the flow path were simulated (Test Case 1). Then, two test cases of laboratory experiments in a straight channel with cohesive fluid were numerically reproduced (Test Case 2). Afterwards, the nonhydrostatic pressure distribution effect during the first instants of snow avalanche release was analysed through a simple test and an experimental-based case (Test Case 3). Finally, the numerical scheme presented here was tested on a real snow avalanche that occurred in the Pyrenees in 2014 (Test Case 4).

3.1. Test case 1: well-balanced equilibrium states

The robustness of the numerical scheme was first tested through the classical 2D asymmetric dam-break presented by Chaudhry (2008), which has been widely used as a benchmark for water flow modelling (Sleigh and others, 1998; Zoppou and Roberts, 2000; Baghlani, 2011). The aim was not to compare the numerical results with the observed ones but to demonstrate how the scheme can achieve a stable solution with a nonhorizontal free surface within a 2D framework.

The analysed case consists of two horizontal domains of 95 × 200 m connected by a 75 m wide and 10 m long conduit, asymmetrically placed 30 m from one of the laterals. The domain was discretised using a structured mesh of square elements with side length of 2.5 m, and a wet-dry threshold of 0.001 m was considered. Considering an initial flow depth of 10 m imposed on one side of the model, two scenarios for the initial condition on the other side were simulated with a maximum time of 20 s: dry conditions and 2.5 m of depth. The fluid rheology was assumed to be purely Voellmy-type with μ = 0.25, ξ = 2000 m s⁻², and a fluid density of 300 kg m⁻³.

The fluid was released and followed a dam-break pattern, producing positive waves that travelled downstream and negative ones that travelled upstream; however, the propagation velocity and shape were characteristic of the fluid's rheology (Figs 5a1, b1). An equilibrium state was reached after 10 s for dry conditions (Fig. 5a1) and 5 s for an initial fluid level of 2.5 m on the other side (Fig. 5b1). A nonhorizontal free surface was obtained (Figs 5a2, b2), without numerical instabilities even when the fluid was stopped, which is a particular state for dense snow avalanches and also characteristic of several non–Newtonian fluid flows.

Figure 5. 2D asymmetric dam-break benchmark. Scenario with initial dry conditions: (a1) evolution of the free surface in cross-section y = 135.5 m, and (a2) 3D representation of the free surface (x5 distortion of the vertical scale) at the end of the simulation. Scenario with 2.5 m of initial conditions: (b1) evolution of the free surface in cross-section y = 135.5 m, and (b2) 3D representation of the free surface (x5 distortion of the vertical scale) at the end of the simulation.

Then, a dam break of a cylinder fluid that was 100 m in diameter, 10 m in height, and centred on a flat surface of 200 × 200 m² was modelled. This tested the treatment of the wet-dry front without considering the effect of the boundary conditions. The domain was discretised using triangle-shaped elements with a side length of 1 m, and a threshold of 0.001 m was chosen for the wet-dry transition. Further, a 20 s simulation was carried out considering a fluid with values of μ = 0.3 and ξ = 1250 m s⁻² for the parameters of the Voellmy friction model.

Figure 6a shows the evolution of the free surface at seven points 20–80 m from the centre of the cylinder (Fig. 6c). The fluid spread radially, preserving a circular shape, and stopped before 10 s of simulation with a nonhorizontal surface and without reaching the boundaries. Further, the slope of the final free surface of the mobilised area decreased radially from almost 0.3 up to 0 (Fig. 6b) due to the effect of the turbulent friction (ξ ≠ 0) during the dynamic stage. As expected, the maximum slope was obtained at the edge of the central area that remained at rest (r ≈ 22 m) because the influence of the turbulent phase was very low (v ≈ 0 m s⁻¹). The slope of the free surface during the static and the dynamic phase was less than μ, which agrees with the physical considerations. The resting free surface at t = 20 s is plotted in Figure 6c with an XZ plane view. In a flat, circular area that was 22 m in diameter, the pile retained its initial depth of 10 m. No numerical instabilities arose, neither at the flat area nor at the wet-dry boundaries, which attests to the robustness of the numerical scheme and the treatment of the wet-dry front.

Figure 6. Numerical results of the circular dam-break: (a) evolution of the free surface considering a radial distance from the centre of the circle of r = 20, 30, 40, 50, 60, 70 and 80 m, (b) slope of the free surface at t = 20 s and (c) XZ plane view of the free surface (vertical scale exaggerated 5 times) at t = 20 s.

The last test aimed to represent the triggering and propagation of an avalanche as well as its collision with an obstacle. It was inspired by observations of avalanches reaching constructions, partially accumulating behind them and partially flowing around them. The model consisted of a straight, inclined channel of 25 m width and 155 m length with two slopes, where the upper part was inclined at 36.87° (125 m) and the lower part was inclined at 1.15° (30 m). In the lower part, 1 m downstream of the slope change, there was a 10 m high and 5 m wide baffle that represents an obstacle (e.g. a building) centred in the deceleration area. An initial volume of 25 m³ of fluid, centred in the channel with a square-shaped area of 25 m² and located 5 m downstream of the upper boundary, was instantaneously released and, then, propagated towards the baffle. The fluid properties were as follows: μ = 0.2 and ξ = 1250 m s⁻² (Voellmy friction model) and a density of 300 kg m⁻³. A grid of square elements with a side length of 0.5 m was used to simulate 25 s of flow propagation with a wet-dry threshold of 0.01 m.

Figure 7 represents the avalanche evolution at 0, 5, 10, 15 and 20 s of the simulation. During the propagation stage, the avalanche spread laterally, nearly reaching the lateral boundaries (notably, no lateral slope was considered), and reached a maximum velocity of up to 12 m s⁻¹. When the avalanche arrived at the baffle, which was placed at the beginning of the deceleration area (with the slope being less than μ), a part of the fluid accumulated at its upstream part with a maximum depth of approximately 0.7 m. The rest of the avalanche passed by the obstacle and continued flowing for 5 m. The avalanche stopped 20 s after its triggering, reaching a static solution with a nonhorizontal free surface and resting partially in the two different slopes.

Figure 7. Map of fluid depth at (a) 0, (b) 5, (c) 10, (d) 15 and (e) 20 s of the simulation (XY plane view). Values above 1 m of depth are represented with the colour of the maximum (deep red).

3.2. Test case 2: a laboratory experiment with a cohesive flow

Dent and Lang (1980) built up a seminatural facility near a ski resort with the aim to analyse the bulk behaviour of dense snow avalanches during the deceleration phase. The facility consisted of a semicircular inclined (30^o) channel lined with a plastic sheet that aimed to achieve maximum flow velocities of up to 18 m s⁻¹ when the snow entered the decelerating zone. A 2.4 m-wide horizontal deceleration area followed immediately after the inclined part. The position of the leading-edge of the avalanche was filmed before being transformed into a function of time. Experiments 1 and 3 (Dent and Lang, 1980), which were associated with terminal velocities of 12 and 18 m s⁻¹, respectively, were analysed.

The computational domain comprised only the flat part of the chute that was discretised with 1D elements of 0.01 m of side. An initial snow volume of 3.36 m³ was introduced in the model with the terminal velocity as an initial condition. Since no data were available, the flow density was assumed to be 300 kg m⁻³. The simulation duration was 3 s and the writing time of the results was 0.1 s, considering a depth threshold of 0.01 m and no-wall friction. Forty-seven combinations of the parameters μ (0.1–0.3), ξ (5,500– 10 000 m s⁻²) and C (490–1060 Pa) were tested per each experiment.

Figure 8a shows the averaged value for all the simulations (black) of the leading-edge position vs time in comparison with the reported data (white) for Experiment 1 and 3. In general, an overestimation of the avalanche front at starting times was observed, but the final position at the end of each experiment (1.8 s for Experiment 1 and 2.4 s for Experiment 3) was well captured for all simulations.

Figure 8. Results of the simulation of the experiments performed by Dent and Lang (1980): (a) leading-edge position of the avalanche of all simulation (mean values) vs time (black) compared to Experiments 1 and 3 (white), (b) evolution of the flow depth and (c) inertial forces for μ = 0.2, ξ = 6500 m s⁻² and C = 815 Pa in Experiment 1 (values above 4000 Pa are coloured in deep red; XY plane view of the first 12.3 m of the channel).

More significant differences can be identified for Experiment 1 than for Experiment 3. A possible explanation for this can be the spreading of the leading-edge position (as in the experiments, the snow flowed as a block). This was probably produced by the initial condition used in the numerical model (which is analysed in detail in Section 3.3). In addition, the simulations of Experiment 1 required 0.3 s more for the avalanche to stop, resulting in a run-over that was approximately 0.25 m larger.

When the flow velocity and the bottom slope are equal to zero, the numerical model can reproduce the avalanche stopping, thereby achieving a stable solution at the end of the simulation. In these situations, the final avalanche depth is not constant along its longitudinal profile (Fig. 8b). The final shape of the avalanche depends on the previous dynamics, in contrast to what happens for water flows. This static final shape of the avalanche is achieved due to the balancing of the source term.

In this case, the analysis of the inertial forces can also help in understanding the evolution of the avalanche shape. As Figure 8b shows, after 1.5 s, the avalanche stops at the rear part, but the rear of the avalanche, still having high velocity, pushes against the front, which has decelerated already (Fig. 8c). This model configuration, with constant velocity as an initial condition, leads to high inertial terms on the rear avalanche part at initial time steps, but then the inertia evolves from the rear to the front part and the velocity decreases accordingly, whereas the depth increases.

Positive results in terms of the leading-edge evolution and final position were obtained for several combinations of the values of the three friction parameters within a relatively large range (μ, 0.1–0.3; ξ, 5,500– 10 000 m s⁻²; C, 490–1060 Pa). The Voellmy–Bartelt friction–cohesion parameters were parameterised through a dependency analysis for the combination of μ, ξ and C that best approximates the observed results of the final leading-edge position for the two experiments. A linear dependency of cohesion on μ and ξ (R ² = 0.97) was obtained for Experiment 1:

(23)$$\eqalign{C = & \alpha \cdot \xi + \beta \cr \alpha = & -0.025\cdot \mu + 0.0265 \cr \beta = & -1 506.4\cdot \mu + 980.66} $$

In the case of Experiment 3, cohesion was found to be log-dependent on the Voellmy parameters, with R² being larger than 0.99. This is shown in Eqn (24), which is as follows:

(24)$$\eqalign{C = & \gamma \cdot ln( \xi ) + \delta \cr \gamma = & -255.32\cdot \mu + 289 \cr \delta = & -216\cdot ln( \mu ) + 1 109.3} $$

Different combinations were used for both experiments, showing the possibilities provided by these parameters for model calibration. In some cases, high ξ values (i.e. up to 10 000 m s⁻²) were required to achieve good results in contrast to the previous test cases. These values were two or three times greater than the values generally found in the literature; however, they are in line with the values presented by others authors in several analyses (Gubler, 1987; Gauer, 2014; Fischer and others, 2015).

3.3. Test case 3: non-hydrostatic anisotropic pressure distribution

Two test cases were used to analyse the effect of nonhydrostatic and anisotropic pressure distribution, which is considered in Eqn (3) with the parameter K _p.

The first test considered a flat channel of 5 m width and 20 m length. A dam-break-like flow was generated by instantaneously releasing fluid from an area of 5 × 5 m², with an initial depth of 2 m at one end of the channel. Two scenarios were calculated after considering a different expression for the friction law: one scenario was calculated using the Manning formula (n ²v ²/h ^4/3), which is traditionally used for water fluid flow in 2D-SWE hydraulic models; whereas the other was calculated by using the velocity-dependent term of the Voellmy friction model (v ²/ξh) with a value of ξ = 1600 m s⁻². Both friction laws are equivalent through ξ = h ^1/3/n ². Particularly for the snow flow, two different K _p values (1 and 0.5) were tested.

This test case is oriented to show the effect of K _p factor in the fluid dynamics; therefore, additional frictional effects have been intentionally limited by omitting the Coulomb and cohesion contributions. This analysis allows comparing the effect of K _p in water flows, but especially in snow flow.

Figure 9a shows the free surface evolution during the first 2 s of both fluids with time increments of 0.5 s. As expected, the flow behaviour was nearly identical for water (blue lines) and snow with K _p = 1 (green dashed lines). For K _p = 0.5 (brown dotted lines), the snow flow was slower, but the free surface acquired similar shapes as for K _p = 1. The pivoting point of the free surface was the same for all simulations, maintaining the length and depth positions at approximately 5 m and 0.9 m, respectively. Further, the inertia terms, presented in Figure 9b, were very similar for water and snow with K _p = 1 and were approximately twice those for K _p = 0.5, which highlighted the effect of K _p on flow propagation.

Figure 9. Effect of the K _p factor on the flow behaviour of water (blue lines) and snow with K _p = 1 (green dashed lines) and K _p = 0.5 (brown dotted lines). Evolution of the (a) free surface and (b) the inertia during the first 2 s, with intervals of 0.5 s in a dam break.

In terms of inertia, lower K _p values imply a lower momentum and lower velocity. When K _p tends to 0, the fluid takes a longer time for the mass to move and deform. K _p = 0 means that the pressure is not transmitted through the mass. Thus, in non–Newtonian fluid flows, such as snow, K _p can be low but not zero.

When only the turbulent parameter of the Voellmy model (ξ) is used (i.e. no solid friction), the fluid continues flowing and does not stop. In this case, given that the friction slope consists only of the turbulent part, i.e. ${\boldsymbol S}_{rh}^{\rm ^{\prime}} ( \xi )$, a centred scheme is used for its representation (as for water) where the friction terms are modelled using the Manning formula, which is also a velocity-dependent equation.

Moreover, a laboratory experiment is used to show the effect of K _p on the avalanche dynamics modelling under controlled conditions. The triggering of avalanches is a complex phenomenon (Schweizer and others, 2003; Gaume and others, 2017) that is beyond the scope of this study; however, it is well known that in nature, during the first-time steps, slab avalanches tend to move like a block. In avalanche modelling, this flow pattern is usually not properly reproduced because it is assumed that the material is fluid from the beginning, which can be quite different from what happens in nature. To bypass this problem, many studies (e.g., Dent and Lang, 1980; Bartelt and others, 1999, 2015) considered an artificial initial condition, i.e. a snow mass with a constant velocity located at a short distance downstream of the release area. Nevertheless, a low value of K _p at the initial steps can help represent the solid-like behaviour of the avalanche near the triggering area with no need to alter the initial conditions. Besides, it is a more physical representation of the phenomenon.

To demonstrate this, Experiment 9 presented by Bartelt and others (2015) was simulated with two different values of the pressure correction: K _p = 1 and K _p = 0.1. An initial analysis was performed using the parameters proposed by Bartelt and others (μ = 0.55, ξ = 2000 m s⁻² and C = 396 Pa). Then, several tests were carried out by varying ξ while maintaining the same value for the others in all the simulations.

The experiment consisted of a 2.5 m-wide straight chute with three different slopes, which were 45, 32 and 1.5° from the upper to the lower part. An initial volume of 13 m³ was released from the first 5 m of the upper part of the chute (Platzer and others, 2007a, 2007b). The domain was discretised by a structured mesh of rectangular elements of 0.01 m in the flow direction with the full channel width in the perpendicular direction. Additionally, a wet-dry threshold of 0.001 m was used.

Figure 10a displays the evolution of the measured shear stress (black squared line), and that simulated by Bartelt and others (black dot line) compared to the results of the proposed scheme (red line). The Coulomb (μ) and turbulent (ξ) contributions in the simulation carried out with the presented scheme are also represented (coloured dashed lines). The results of the shear stress that best match the observations are for K _p = 0.1 compared to K _p = 1 (Fig. 10a). When K _p = 0.1, the snow arrives at the measuring point with a block-like shape, with higher depths and lower spreading (Fig. 10b2, red line). Hence, modifications on the K _p factor allow for controlling the dispersion of the avalanche, which are related to the speed of the wave (K_p < 1 implies less spread). This would also be valid at later stages if this parameter changes with time, the fluid properties or the avalanche dynamics (e.g. depth, velocity, etc.).

Figure 10. Effect of the K _p factor on the flow for Experiment 9 of Bartelt and others (2015). Observed and simulated results of (a) the shear stress and (b) flow depth by Bartelt and others (black dot-line) and with the proposed numerical scheme for (1) K _p = 1 and (2) K _p = 0.1 10 m downstrem of the release area.

Furthermore, different tests were performed after considering ξ values of 500 and 1000 m s⁻². In these cases, only the depth evolution is plotted (Fig. 10b). As a consequence of the velocity reduction due to lower values of ξ, the proposed numerical scheme exhibits small differences in the arrival time of the avalanche, particularly for K _p = 1 (<0.3 s). For the particular case of ξ = 500 m s⁻² and K _p = 1, part of the avalanche was stopped near the end of the flume (Fig. 10b1, orange dashed line). This behaviour was not produced for K _p = 0.1, as higher inertia values were achieved, such as for ξ = 2000 m s⁻² (Fig. 10a2).

Both numerical simulation results differ considerably from the observed data, especially when K _p = 1 for the proposed scheme. When nonhydrostatic conditions are applied (K _p = 0.1), the shear stress obtained with the proposed scheme reproduces the measurements more closely (Figu. 10a2). This is also observed in the depth profile (Fig. 10b2). The numerical results presented by Bartelt and others show that the bulk arrives at the measured point 0.4 s faster than the observations (Fig. 10a), probably due to the assumption of an initial condition. In the proposed scheme, this gap is reduced to 0.3 s for K _p = 1, while it is reduced to 0.1 s for K _p  = 0.1 (ξ = 2000 m s⁻²), which demonstrates the benefits of using nonhydrostatic anisotropic pressure distribution in this particular test. However, values lower than 1 could not be well-suited once the avalanche fluidified, which behaves more like an ideal fluid.

The results of the proposed numerical scheme were obtained by simulating the entire experiment, in contrast to the results presented by Bartelt and others, where an initial velocity was assumed. In this regard, assuming the flow depth and velocity at each time as a mean value for the 2.5 m-wide chute, the simulations of Bartelt and others (2015) were probably carried out with a higher volume (estimated at approximately 32 m³). Additionally, the avalanche simulated by Bartelt and others did not follow a ‘block-like’ behaviour, as the flow depth decreases smoothly (Fig. 10b), which again highlights the relevance of considering K _p values lower than 1, at least for the first stages of the avalanche releasing.

3.4. Test case 4: Bonaigua avalanche zone

The Bonaigua mountain pass, located on the southern side of the Pyrenees range (Catalonia, Spain), is one of the main access routes to the Val d'Aran and the Baqueira Beret ski resort (Fig. 11a). According to the Avalanche Database of Catalonia, BDAC (ICGC, 2020), there are 85 avalanche paths in the Bonaigua Valley (Fig. 11b, highlighted in yellow). Since 1986, more than 500 avalanches have been registered within Bonaigua, and 62 of them correspond to the path called BNG045 (Fig. 11c, black polygon).

Figure 11. The Bonaigua case study: (a) location of the study area (red point) (background image source: Copernicus Land Monitoring Service), (b) avalanche paths identified in Bonaigua Valley (yellow polygons), (c) BNG045 potential avalanche path (black polygon) and observed avalanches (blue polygons) and (d) observed avalanche of 26 January 2014 (green polygon), and estimated release area (red polygon) (background image source: Institut Cartogràfic i Geològic de Catalunya [CC by 4.0]).

On 26 January 2014 a slab avalanche, which was rather smooth, was observed within the BNG045 path and its extent was later mapped (Fig. 11d, green polygon). The avalanche crossed the road C-28 and stopped several metres below it after travelling a distance of approximately 650 m, with the deposition area being split into two branches. The return period of this event was estimated at 30 years, as this was the largest avalanche-related event in ~35 years. Neither the exact starting area nor the fracture depth could be measured in situ. According to this return period, the associated snow depth was estimated at 1.09 m. The fracture depth was approximated by evaluating the DH3dd in accordance with Burkard and Salm (1994), which resulted in 0.64 m for the release area. Despite the lack of some important field data, the performance of the model in this case can be investigated, and insights can be gained into its dynamics.

An unstructured mesh of triangles with a side length of 2 m was used to discretise the study area, which resulted in more than 185 000 calculation points or elements. An initial condition of 0.64 m of snow, which represents an initial volume of 16 230 m³, was imposed on the release area. The flow density was assumed to be 300 kg m⁻³ (for which no data were available), and a depth threshold limit of 0.01 m was chosen. The maximum simulation time was 180 s.

The performance and applicability of the numerical model were evaluated under a real snow avalanche in two phases. First, a performance test was conducted through the Monte Carlo (MC) method by varying the values of the friction–cohesion parameters, which are: μ (0.1–0.4), ξ (250–2000 m s⁻²), and C (0–500 Pa). The MC method uses a rather large number of simulations with all uncertain input parameters being sampled randomly. The parameters were assumed to be uniform in space and time in the 1296 scenarios simulated. On the other hand, accounting for the wide range of combinations of the friction–cohesion parameters that provide similar results, a sensitivity analysis of the runout area to the friction parameters was presented. This kind of approach can be useful to determine the areas where the flow concentrates and could be a key component in a tool for probabilistic hazard mapping.

The results of the performance test can be grouped into four situations depending on the location of the deposit: (1) the avalanche stopped before the road, (2) the avalanche stopped on the road, (3) the avalanche stopped after the road, and (4) the avalanche continued flowing. Approximately 44% of the simulated avalanches stopped either on the road or a few metres below of it. In scenarios with μ equal to or lower than 0.2, this percentage increased up to 72%. Further, this model configuration (μ ≤ 0.2) also generated all the scenarios where the avalanche continued flowing (8% of the total). The values of the parameter μ greater than 0.2 generated the detention of avalanche before the road (83%), whereas in this case, the avalanche reached and stopped on the road only in 17% of the simulated scenarios.

Table 1 exemplifies the main results obtained from the performance test (each scenario is labelled with the values of the parameter's combination as μ_ ξ_ C). The numerical results of situations 2, 3 and 4 show accumulations from a few decimetres to less than 3 m. Moreover, the scenario 0.1_1000_0, which stopped after the road (situation 3), showed the maximum snow depth on the road at 2.7 m. The scenarios with low resistance contribution such as 0.1_2000_0 allowed the avalanche to continue flowing. Furthermore, the maximum velocity (approximately 40 m s⁻¹) was obtained in scenario 0.1_2000_0 (situation 4), so the simulated velocities (<40 m s⁻¹) were reasonable. Besides, the accumulations of snow above 1 m were obtained on the road at the end of the simulation, which reinforces the good performance of the numerical model in front of topographical irregularities.

Table 1. Summary of the extreme values of each parameter's combination of the 1296 scenarios simulated

nd, no data.

Each scenario is labelled with the values of the parameters' combination as μ_ ξ_ C _.

The combination of parameters that provided a similar runout distance than the observations ranges widely, with values of μ [0.1–0.2], ξ[250–2000] and C[0–500]. To show the performance of the numerical model, particular results of the scenario μ = 0.125, ξ = 500 m s⁻², and C = 0 Pa are plotted in Figure 12. With this configuration, the model was able to simulate the following specific behaviour of the observed avalanche: the avalanche crossed the road and the deposition area was divided into two branches. The particular change of direction in the avalanche deposit observed between 40 and 100 s (white rows) was produced because the first part of the avalanche stopped and generated a huge deposit on the east side of the road (Fig. 12b). As the rear part of the avalanche continued flowing, when it reached this deposit the avalanche turned to the west generating a new deposit. The velocity rows of Figure 12 (bottom) show not only the direction of the velocity modulus but also how the avalanche is displacing. If the mass and energy of the rear part of the avalanche would be large enough, the existing deposit could be modified, allowing for the avalanche to continue flowing in this direction.

Figure 12. Maps of maximum snow depth (top) and velocity (bottom) for the parameters combination of μ = 0.125, ξ = 500 m s⁻², and C = 0 Pa. Evolution of the variables at: (a) 40 s, (b) 100 s and (c) 140 s. For the velocity maps (bottom), white rows represent the direction of the velocity modulus (background image source: Institut Cartogràfic i Geològic de Catalunya [CC by 4.0]). The transparent cyan polygon represents the observations.

The differences between the simulations and the observed data can be attributed to the uncertainties and assumptions regarding the release area (shape, extension and depth) as well as the topography used. The terrain elevation data corresponded to the summer topography (as the terrain elevation did not consider the snow pack during the event) and may differ in some areas from the one during the event, thereby producing different avalanche trajectories (Sanz-Ramos and others, 2021). The existence of two branches in the deposit area could be due to a single release area, but the triggering of two or more avalanches, whether simultaneously or not, cannot be discarded, especially due to the combination of the low values of the friction–cohesion terms in the simulations that can reproduce this behaviour. As previously stated, if a previous deposit exists, the dynamics of later avalanches could be conditioned. In any case, the scheme performed without numerical instabilities for all the simulations, and the avalanche flowed and stopped in accordance with the topographical data and the values of the rheological model used.

Figure 13a shows, though the overlap of the maximum extent of all simulated avalanches, the sensitivity of the runout area to the friction parameters. The bright cyan areas indicate the locations where the extent of the avalanche was achieved a high number of times, while light cyan (i.e. almost transparent) areas were obtained in fewer simulations. The first metres of the simulated avalanche paths agreed with the observations. Moreover, the narrowing in the avalanche path observed a few metres before the road (green polygon) presented worse results with a clear change of direction towards the south. Further, the surrounding of the road presented the lowest overlapped areas of the avalanches. Some simulations reached and overtopped the road (Fig. 13a, light cyan) due to low values of the resistance parameters (μ ↓ – ξ ↑ – C ↓). Only a handful of them were able to reproduce the particular deposit of the observed avalanche, which was split into two branches.

Figure 13. Sensitivity of the runout area to the friction parameters. (a) Overlapped runout areas of the 1296 simulations (cyan areas). Maps of conditional probability of snow depth greater than (b) 0.5 m, (c) 1 m, (d) 1.5 m and (e) 2.0 m for the indicated avalanche's triggering area and release volume (background image source: Institut Cartogràfic i Geològic de Catalunya [CC by 4.0]).

Assuming the previously indicated triggering area and release volume, and considering that the friction–cohesion parameters (μ, ξ and C) have uniform probability distribution functions, maps of conditional hit probability of avalanches can be calculated. In other words, the conditional probability of an avalanche reaching a point was evaluated whether a threshold of snow depth was exceeded in any element of the calculation domain and scenario or not. In the affirmative case, class 1 was assigned to the element. Accounting for the number of class 1 of each element in all the scenarios, the conditional probability of exceeding a particular threshold was directly obtained at each element.

A snow depth threshold equal to or larger than ɛ_wd will provide the conditional probability of an avalanche reaching a point, given there is an avalanche release under the previous assumptions. Figures 13b–d show the conditional probability of exceeding a snow depth of 0.5, 1, 1.5 and 2.0 m, respectively. These maps reveal the areas where the flow would concentrate during the dynamic phase and where the avalanche would potentially accumulate in the deposit. In this regard, a first potential stop area would be placed at the middle part of the avalanche path, before the narrowing zone of the observed avalanche, due to the low uncertainty of overcoming a snow depth of 2.0 m (Fig. 13e) within the stated intervals of the friction–cohesion parameters. This is consistent with the observed avalanches (Fig. 11c), where most of the triggered ones from the top of the release area stopped there. Avalanches that have enough energy to pass through this area (with a decreased slope) would potentially stop either at the east side of the road or a few metres below (Fig. 13d), as the observations show.

This analysis agrees with the observed slab avalanche of 26 January 2014, which was triggered from an unknown release area and produced two separate deposits on the east and the west side of the road. The generation of an intermediate branch that bypasses the narrowing zone by the west (Fig. 13c) is less probable and only possible for high-energy avalanches and if there is a previous deposit in the middle part of the avalanche path that pushes the avalanche to the west. A high degree of uncertainty is produced on the boundary areas, especially from the mid part of the avalanche path to the potential deposit area after the narrowing zone (Fig. 13b).

4. Conclusions

The numerical modelling of non–Newtonian fluid flows based on the solution of the 2D-SWE requires not only that the rheological model be changed but also that the numerical scheme used to solve the equations be adapted. These modifications must ensure a proper balance between the homogeneous part of the equations and the source term, which contains the rheological model. Although the Voellmy–Bartelt friction–cohesion model was chosen for this study as a reference, the presented developments are also valid for other models with a similar yield-turbulent stress structure (e.g. the Bingham fluid model).

The proposed numerical scheme for non–Newtonian fluid modelling with Iber is based on the consideration of the flow propagation properties (characteristics theory for hyperbolic systems) and ensures the balance among the different terms of SWE. On the one hand, only the parameters of the source term of the SWE that represent forces that are counterbalanced with the pressure forces must be treated with upwind discretisation, while on the other hand, a centred discretisation must be used for the velocity-dependent terms of the rheological model, in the same way as it is for water flows.

The performance of the numerical scheme was proved with a commonly used snow avalanche friction model (Voellmy) and, in some cases, together with a recently developed cohesion model (Bartelt). The test cases revealed that the numerical scheme can achieve stable and consistent solutions even when the avalanche stops on a sloping terrain, as the balance of gravity forces and friction is properly represented.

The presented numerical scheme is well-balanced, even when including the pressure factor parameter K _p, which facilitates the consideration of nonhydrostatic and anisotropic pressure distribution. Therefore, the results show that the consideration of nonhydrostatic pressure, at the first steps of the avalanche motion when it has a block-like behaviour, is crucial for a good representation of its dynamics. The inclusion of K _p modifies the SWE, which are derived by assuming hydrostatic pressure and, thus, predict a more fluid-like behaviour than expected for snow flows. As has been shown here, the hypothesis of nonhydrostatic and anisotropic pressure can improve this behaviour, provided that the K _p parameter is properly estimated.