Snow avalanche friction relation based on extended kinetic theory

Rheological models for granular materials play an important role in the numerical simulation of dry dense snow avalanches. This article describes the application of a physically based model from the field of kinetic theory to snow avalanche simulations. The fundamental structure of the so-called extended kinetic theory is outlined and the decisive model behavior for avalanches is identified. A simplified relation, covering the basic features of the extended kinetic theory, is developed and implemented into an operational avalanche simulation software. To test the obtained friction relation, simulation results are compared to velocity and runout observations of avalanches, recorded from different field tests. As reference we utilize a classic phenomenological friction relation, which is commonly applied for hazard estimation. The quantitative comparison is based on the combination of normalized residuals of different observation variables in order to take into account the quality of the simulations in various regards. It is demonstrated that the extended kinetic theory provides a physically based explanation for the structure of phenomenological friction relations. The friction relation derived with the help of the extended kinetic theory shows advantages to the classic phenomenological friction, in particular when different events and various observation variables are investigated.


Introduction
Within the past few decades several software tools for the simulation of snow avalanches or, generally speaking, shallow granular flows have been developed, such as SamosAT (Sampl and Zwinger, 2004;Zwinger et al., 2003), TITAN2D (Pitman et al., 2003;Patra et al., 2005), RAMMS (Christen et al., 2010) or r.avaflow (Mergili et al., 2012).In this study the software SamosAT is utilized.The implemented flow model therein is the Savage-Hutter model (Savage andHutter, 1989, 1991), which is related to the shallow water or Saint-Venant equations (de Saint-Venant, 1871).These models idealize avalanches and other free surface flows as depthaveraged flows.The Saint-Venant equations are set up in a Cartesian coordinate system and the normal stresses are assumed to be hydrostatic.In contrast, the Savage-Hutter equations are set up in a curvilinear coordinate system (compare, e.g., Bouchut and Westdickenberg, 2004) and the hydrostatic pressure assumption is replaced by an assumption for the lateral active or passive earth pressure, common in soil mechanics.The density is assumed to be constant in both models.
Within this framework rheological models attract a significant portion of attention.A widespread, classic phenomenological rheological model still utilized in depth-averaged models is the Voellmy friction relation (Voellmy, 1955).A motivation for this phenomenological friction relation, based on a physical model following Bagnold (1954;1966), is presented by Salm (1993).The underlying model of Bagnold (1954) represents a specialization of kinetic theory (Mitarai and Nakanishi, 2005;Lee and Huang, 2010).Another approach is presented by Issler and Gauer (2008), who apply the Norem-Irgens-Schieldrop model (Norem et al., 1987) to establish a rheological model for snow avalanches, yielding a similar relation as in this work.More recently, Buser and Bartelt (2009) introduced the concept of random kinetic energy to snow avalanche simulations.The energy evolution is described by a transport equation and influences the frictional behavior along the avalanche path.This approach shows some similarities to the kinetic theory model used in this work.
The rheology of granular materials has been investigated in many scientific works within the framework of three-Published by Copernicus Publications on behalf of the European Geosciences Union.dimensional continuum mechanics.An important category of microrheological models, dealing with rapid granular flows, is the kinetic theory (Campbell, 1990;Goldhirsch, 2003).Standard kinetic theory struggles to describe the dense flow regime at low shear rates.Recently developed extensions aim to take into account the formation of clusters (Jenkins, 2006(Jenkins, , 2007) ) and enduring force chains (Berzi et al., 2011;Vescovi et al., 2013).The basis of the presented approach is the extended kinetic theory, as formulated by Vescovi et al. (2013).This microrheological model deals with both the quasi-static regime, described by the critical state theory, and the rapid, collisional flow, described by the kinetic theory of granular gases.
To implement the constitutive model into the depthaveraged dynamic model, several assumptions about the vertical structure of the flow are made, yielding a simplified friction relation.It is shown that the simplified expression is similar to classic friction relations.In a further step the obtained relation is compared to classic friction relations which are often applied in natural hazard estimation.For this task, back calculations of well-documented avalanches are conducted and minimal residuals in multiple observation variables, such as runout distance, affected area and velocity, are determined.It is shown that the relation obtained by kinetic theory allows to reduce the residuals for the presented events.

Constitutive relations in the Savage-Hutter model
The governing equations of the Savage-Hutter model, extended for entrainment, as implemented in SamosAT, for an incompressible, granular flow over a one-dimensional terrain can be expressed as These expressions can also be written for the twodimensional case (e.g., Sampl and Zwinger, 2004).Equation (1) describes the conservation of mass and Eq. ( 2) the conservation of momentum parallel to the slope surface.x is the curvilinear coordinate, z the coordinate perpendicular to x and t the time.ρ represents the flow density, assumed to be constant, h the slope perpendicular flow depth, u the slope parallel, depth-integrated velocity in x direction, 1 h h 0 u x dz and g x and g z are the gravitational accelerations in x and z direction, respectively.The entrainment rate q represents the mass entrained by the avalanche within a specific amount of time and area (with unit kg m −2 s −1 ).Note that Eqs. ( 1) and ( 2) are only valid for small curvatures.
The resistance of the material against its deformation is considered with the second and third term on the right-hand side of Eq. ( 2).The second term represents the slope parallel stress gradient ∂σ x ∂x , expressed by the basal normal stress σ b and the earth pressure coefficient K a/p .σ b is calculated with respect to the centripetal acceleration due to the basal curvature κ as (Sampl and Zwinger, 2004;Fischer et al., 2012) σ b = ρ h g z + κ u 2 . (3) The earth pressure coefficient K a/p is given as following a Mohr-Coulomb yield criterion (Savage andHutter, 1989, 1991).Another approach by Salm (1966Salm ( , 1993) ) is based on Rankine's earth pressure.Here, φ is the internal friction angle and δ is the bed friction angle.K a/p = 1 coincides with hydrostatic pressure and the Saint-Venant assumption and is a commonly used simplification (Christen et al., 2010).Bartelt et al. (1999) showed that the sensitivity of the relevant simulation results to the internal friction angle is rather small for snow avalanches.A detailed analysis is not in the scope of this work and a fixed internal friction angle, employing Eq. ( 4) is utilized.
The third term on the right-hand side of Eq. ( 2) describes the basal friction.Usually this term is a combination of a Coulomb type friction, σ b tan δ and a velocity-dependent drag term, f (σ b , u) u 2 (Hutter et al., 2005).

Constitutive model
In this section we outline the fundamental structure of the extended kinetic theory for simple shear conditions.Classic kinetic theories are based on the statistical description of binary particle collisions, which are the governing processes at low volume fractions.The extended kinetic theory used here includes extensions (Jenkins, 2006(Jenkins, , 2007;;Berzi et al., 2011;Vescovi et al., 2013) to make the model suitable for high volume fractions and therefore dense flows.In this section we employ material parameters for idealized quartz sand to identify the basic features of the extended kinetic theory.In Sect.6.5 we give estimations for material parameters suitable for snow, based on back calculations of real avalanche events.

Critical state theory
The critical state theory was developed in soil mechanics, where the deformation rates are very small and therefore no inertial forces occur.This theory can be outlined with the help of the simple shear test with constant normal stress σ = 0, which is similar to the situation sketched in Fig. 1.
A soil sample is sheared by an increasing shear angle γ .The sample has an initial density, and with that an initial concentration ν 0 = V p /V , with the volume of the particles V p and the total volume of the sample V .Note that in soil mechanics it is more common to use the void ratio e = (1 − ν)/ν.The shear stress increases due to the imposed shearing.An initially dense sample increases its total volume (see Fig. 1a) and an initially loose sample decreases its volume (see Fig. 1b).However, the concentration of both samples of the same material approach the same value, the so-called critical concentration ν ss and this concentration remains constant for further shearing (see Fig. 1c), i.e., ν = 0.At this stage the shear stress remains constant too, i.e., τ = 0.
The shear stress at critical state is related to the applied normal stress: τ = σ tan φ ss , with the critical (steady state) friction angle φ ss .The critical concentration depends on the applied normal stress σ , i.e., for tests with a higher normal stress a larger critical concentration (a higher critical density) will be reached, which constitutes a relation between normal stress and concentration at critical state.In turn the normal stress can be regarded as function of the critical concentration: σ (ν ss ).It is important to know one limiting value of this relation: the concentration of any granular media at rest cannot be less than a minimal concentration ν rlp (random loose package or maximal void ratio e max in soil mechanics), depending on the material.The critical concentration is bounded by this limit.Physical experiments to determine this limit are performed at very low stress levels, virtually zero stress level, and the relation between the normal stress and the concentration at critical state has to provide σ = 0 for ν ss = ν rlp .Note that τ = 0 in this limiting state, too.Since the shearing in any flowing mass is high, it can be assumed that the critical state is reached (compare Figs. 1c and 2).With this assumption, the concentration in the flow is ν = ν ss , and we will drop the index ss in the following.

Kinetic theory
The second constituent of the applied theory is the kinetic theory, which provides relations for the normal and shear stress based on statistical considerations of the collisions between particles in the volume V for high shear rates γ .It is reasonable to assume that the particle collisions will lead to high stresses for high shear rates and high concentrations.Appropriate relations can be established to provide both stresses σ c (ν, γ ) and τ c (ν, γ ) as functions of the concentration and the shear rate.These functions are very complex and will be given in the next section.So far, it is sufficient to keep the qualitative behavior in mind: σ c and τ c will approach zero for a vanishing shear rate and will decrease for a decreasing concentration.
The main idea of the extended kinetic theory by Vescovi et al. (2013) is to combine the quasi-static stresses for vanishing shear rates at critical state (denoted by q) with the collisional stresses (denoted by c) of the kinetic theory: It is important to note that the concentration is the critical concentration in both stress contributions.For concentrations lower than the loosest quasi-static limit ν rlp , the quasi-static stresses will vanish, and for a vanishing shear rate the collisional stresses will vanish.As such, critical state behavior will be ensured for low shear rates and purely collisional behavior will be ensured for densities lower than any statically possible density (which is not considered as dense flow herein).At the bottom, the normal stress σ (Eq.5) has to be in equilibrium with the basal normal stress σ b (Eq. 3) of the flow model.This can be used to calculate the basal shear stress τ b as shown in Sect. 4. Equation (3) can be extended by an additional term, accounting for vertical acceleration due to dilatancy, called dispersive pressure (Buser and Bartelt, 2015;Bartelt and Buser, 2016;Bartelt et al., 2016), but is neglected in this work.Another approach to include dilatancy in depthaveraged granular flow models can be found in Iverson and George (2014).

Model formulation
The simple shear setup, on which the presented model has been developed, is shown in Fig. 2. In this setup, the relevant stresses and deformations are constant over the whole domain.Consequently, only local field values have an influence on the constitutive behavior.This matches arbitrary fields when zeroth-order expansions of spatial and temporal gradients are employed (e.g., Garzó and Dufty, 1999).Moreover, no diffusive or convective transport is present in this setup.
Flows of granular material can display a large span of grain concentrations.Microscopic mechanical processes and consequently the macroscopic behavior of the material changes substantially with the concentration or solid fraction ν and the granular temperature T .The concentration (or solid volume fraction) is statistically defined as where d is the particle diameter and n the number of particles per unit volume at position x and time t, given as where f (v, x, t) is the single particle velocity distribution function, defined as the number of particles at position x with velocity v at time t.Consequently f (v, x, t)/n(x, t) is the probability density for a particle at position x and time t to have the velocity v (e.g., Goldhirsch, 2003;Jenkins and Berzi, 2010).Two different approaches are used to derive f (v, x, t), the Boltzmann equation and the Enskog equation (Garzó and Dufty, 1999;Vescovi, 2014).The granular tem-perature is associated with the fluctuation of the particle velocity, where the mean velocity u is defined as Although the rheological model assumes simple shear conditions, the fluctuating velocity is three-dimensional, leading to the factor 3/2 in Eq. ( 9).Concentration and granular temperature are determined by the hydrodynamic fields and their gradients.For the presented simple shear setup, these are the normal stress along the transversal direction σ and the shear rate γ , given as the derivative of the velocity in its perpendicular direction, ∂u ∂z (see Fig. 2).
To describe the whole range of flow configurations, multiple mechanical processes, described by different theories, need to be taken into account (Berzi et al., 2011;Vescovi et al., 2013).
On the one hand, the critical state theory (Roscoe et al., 1958;Schofield and Wroth, 1968) describes granular material at vanishingly small shear rates γ .This model is completely time independent and does not take into account the velocity of any process.The stresses in the material are based completely on enduring force chains between particles.Also the assumption of the incompressibility of granular flows follows from this theory: granular material under motion always reaches asymptotically a certain stress dependent concentration, the critical concentration.The contribution of quasistatic stresses, e.g., force chains, to the total stresses, as described by the critical state theory is given as with ( 13) Equations ( 11) and ( 12) can be considered as critical state line in the ν-σ -τ space.Material parameters are the tangent of the internal friction angle at the critical state tan φ ss , the Young's modulus of the particles E, the particle diameter d, the concentration at random loose packing ν rlp , the concentration at closest packing ν s and the dimensionless parameter a.
On the other hand, the kinetic theory of granular gases describes the granular material under the influence of high shear rates.The stresses in the material are based on short Nat.Hazards Earth Syst.Sci., 16,2016 www.nat-hazards-earth-syst-sci.net/16/2325/2016/contacts between particles, i.e., elastoplastic collisions.The following form of this theory was proposed by Garzó and Dufty (1999) and extended by Jenkins and Berzi (2010) and Vescovi et al. (2013).It is limited to homogeneous, steady, simple shear flows of identical, dry, spherical particles.The contribution from collisions to the total stresses is given as and the dissipation rate of the granular temperature is given as with Additional material parameters, introduced by the kinetic theory, are the particle density ρ p , the coefficient of restitution and the dimensionless parameter c.The parameter set is summarized in Table 1.The functions f 1 , f 2 and f 3 are solely dependent on the concentration ν and the material parameters; the function f 4 takes into account the particle stiffness.g 0 is the radial distribution function, following the approach of Carnahan and Starling (1969) for ν ≤ 0.49 and the approach of Garzó and Dufty (1999) for ν > 0.49.L is the correlation length, accounting for correlated motion of Parameter for correlated motionparticles at high concentrations and s is the mean separation distance among particles (i.e., the mean free path between collisions) (Vescovi et al., 2013).The total stresses can be expressed as sum of the quasistatic and the collisional stresses, combining critical state theory and kinetic theory: Coupling of quasi-static force chains and collisions is ensured by sharing the concentration (compare Eqs. 5 and 6).
The concentration determines the dominant process, e.g., quasi-static force chains at high concentrations and collisions at low concentrations.Blending between those limit cases is based on the well-known concentration dependence of both processes.The evolution of the granular temperature is described by the conservation equation, with production τ c γ , flux q c and dissipation c .The assumptions of steady state and simple shear imply that granular temperature is dissipated where it has been produced.This approach, called equilibrium assumption (van Wachem, 2000), is also applied to dense flows apart from steady, simple shear conditions (Syamlal et al., 1993;Boemer et al., 1997;van Wachem et al., 1998van Wachem et al., , 1999;;van Wachem, 2000) and is justified by the dominance of generation and dissipation terms over convection and diffusion terms.Other theories do not assume equilibrium; e.g., for a depth-averaged granular temperature transport equation, see Christen et al. (2010).With the equilibrium assumption, the transport equation can be reduced to an algebraic formulation of the equilibrium state, given as (Vescovi et al., 2013) This assumption allows to apply the kinetic theory to operational simulation software without introducing additional differential equations.By introducing Eqs. ( 17) and (18) into Eq.( 33), the granular temperature can be expressed as a function of the shear rate γ and the concentration ν: Introducing Eq. ( 34) in Eqs. ( 30) and ( 31) leads to an expression for the total stresses, only depending on γ and ν: According to Eqs. ( 35) and ( 36), it is possible to characterize the flow regime with only two state variables.In the case of known values for σ and γ , Eq. ( 35) can be used to solve for ν, using Newton-Raphson (e.g., Press et al., 1996) or another root-finding routine.τ can then be calculated with Eq. ( 36), T with Eq. ( 34), if required.Material parameters for snow are not available at the current stage.However, they can be estimated from the results of back calculations of real-scale avalanches, conducted in the second part of this paper.To qualitatively highlight the most important features of the constitutive model, it is analyzed for an idealized 1 mm quartz sand (d = 1 mm, ρ p = 2600 kg m −3 , K = 2.8 × 10 7 Pa m, = 0.6, c = 0.5, a = 1.8 × 10 −6 , tan φ ss = 0.5, ν s = 0.619, ν rlp = 0.55; Vescovi et al., 2013, see Figs. 3-10).Basic features identified in this section will be approximated with a simplified relation to reduce the amount of material parameters, increase the numerical stability and reduce the computational effort.
Figure 3 shows the dynamic critical state surface in the σγ -ν space.According to the presented theory, flow states are limited to this surface.The color in Fig. 3 shows the dominant source of stresses, which can be interpreted as flow regime.In yellow areas, enduring contacts, forming elastic networks between particles, are dominant (referred to as quasi-static regime in da Cruz et al.Vescovi et al., 2013, or elastic-quasi-static regime in Campbell, 2002, 2005, 2006).In blue areas, collisional stresses are dominant (referred to as collisional regime in da Cruz et al., 2005;Vescovi et al., 2013;Berzi et al., 2011, kinetic regime in da Cruz et al., 2005;Forterre and Pouliquen, 2008;Vescovi et al., 2013;Berzi et al., 2011, or inertial-collisional regime in Campbell, 2002, 2005, 2006).The flow is purely collisional for concentrations below the random loose package: ν < ν rlp .The red area represents a transitional zone between those cases, where both effects co-exist (referred to as dense regime in da Cruz et al., 2005, elastic-inertial regime in Campbell, 2002, 2005, 2006, or intermediate regime in Vescovi, 2014).
The granular temperature, which is not shown, is almost solely dependent on the shear rate γ .
The stress ratio τ/σ and the respective concentration ν for a given set of σ and γ are shown in Figs. 4 and 5.The constitutive model predicts a stress ratio τ/σ = tan φ ss for γ → 0, corresponding to critical state theory and a stress ratio corresponding to kinetic theory for high shear rates.The concentration is almost unaffected by changes in shear rate γ until it drops below ν rlp , e.g., at γ ≈ 5 × 10 2 s −1 for σ = 10 3 Pa.At this point the concentration decreases abruptly and the shear stress lowers with increasing shear rate, after reaching its peak.This behavior can be interpreted as a transition from a dense flow to a powder cloud like flow.This work is focused on dense flow -the post-peak behavior is not investigated.
A separation of quasi-static and collisional stresses is shown in Figs. 6 and 7.For small stress levels, increasing collisional stresses τ c can compensate the decreasing quasistatic stresses τ q .At high stress levels, σ > σ * (ν > ν * ), this is not the case and the stress ratio τ σ = τ q +τ c σ shows a nonmonotonic behavior before reaching the peak (Fig. 4).The thresholds for the non-monotonic behavior, σ * and ν * can be calculated as (Vescovi et al., 2013) where In the following, we focus on the monotonic dense flow regime, where It is not feasible to implement the complete extended kinetic theory model in an operational simulation tool.The main reason is that for certain combinations of normal stress σ and shear rate γ , no solution for the concentration ν can be obtained (e.g., white space in Fig. 10), which makes the numerical treatment of the whole flow model unstable.Moreover, the procedure to solve the constitutive equations requires an iterative root-finding method, which is computational expensive.Another difficulty arises because the velocity profile is required to link the shear-rate-dependent rheology model with the depth-integrated flow model.The velocity profile cannot be calculated analytically for the complex kinetic theory model.Numerical simulations of the velocity profile, employing the full extended kinetic theory can be found in Rauter (2015).In this work we stick to the analytical description of the velocity profile.Because of these reasons, approximations of Eqs. ( 35) and (36) are made.In analogy to other studies (see Ancey, 2007, for a review) two approaches, varying the separation of friction into two parts, are evaluated.The first approach separates friction by their source (quasi-static-collisional), while the second approach separates friction into a shear-rate-independent and shearrate-dependent part.The first approach is given as with Within this formulation, λ γ 2 represents solely collisional stresses.The decrease of quasi-static stresses is considered with the share rate dependence of µ.
The second approach is given as with Here, the term λ γ 2 also accounts for the decreasing quasistatic stress.Both λ and λ are approximately constant within a certain range of σ and γ .The second approach with constant µ = tan φ ss leads to a better approximation (see Fig. 8) and a simpler formulation with less parameters.However, the non-monotonic behavior at high stress levels cannot be reproduced with this approach.Values for λ are shown in Fig. 10.Up to a normal stress of 10 4 Pa and until the peak is reached, λ can be approximated as constant.The resulting simplified relation is shown and compared to the full kinetic theory in Fig. 9. Thus, Eq. ( 43) with a constant value for λ is employed in the following: 4 Velocity profile and kinematic relations The constitutive model obtained by the granular kinetic theory results in a relation depending on the shear rate, which does not explicitly appear in depth-averaged models.However, the equilibrium of stresses at the bottom of the avalanche requires that where σ b is the normal stress at the bottom and γb is the shear rate at the bottom.According to the Savage-Hutter model and related friction models, the friction may depend on velocity.To obtain an expression of the form τ b (u, h) we need to express σ b (see Eq. 3) and γb with known flow variables.Therefore a reconstruction of the velocity profile is required.Supposing that the avalanche has reached its steady state on a slope with constant inclination α and that ∂h ∂x is very  ) for the idealized quartz sand.The simplified relation for the stress ratio approximates the complex extended kinetic theory good for concentrations between ν rlp and ν * , e.g., in the monotonic dense flow regime (blue area, top-right).Between ν * and ν s the regime is dense but non-monotonic (red area, top-right), below ν rlp the regime is purely collisional (yellow area, top-right) and the concentration decreases rapidly with increasing shear rate.small, as for the middle part of an avalanche (referred to as equilibrium shape of the velocity profile in Issler and Gauer, 2008, or simple shear infinite landslide model in Dutto, 2014), all volume forces and stresses can be expressed with the differential equations The left-hand side in Eqs. ( 47) and ( 48) describes the change of stresses in z direction, which is caused by the gravitational volume force (right-hand side).Introducing the constitutive relation (Eq.45) in Eq. ( 47) leads to Figure 10.The factor λ as a function of σ and γ .In the white area no solutions could be obtained with the model.In the yellow area the value for the factor λ has its maximum value of approximately 3.5 × 10 −4 Pa s 2 .This value decreases at the borders of the yellow area.In the gray area, values for λ are negative, indicating a nonmonotonic behavior.For comparison, velocity measurements of a dry snow avalanche from a field test in Vallée de la Sionne is shown (red filled circle is the mean velocity; bars are fluctuations) (data from Kern et al., 2009;Sovilla et al., 2015).The measurement shows a more bulbous velocity profile, which indicates a plug flow regime.The bars show the high fluctuation of velocity of grains, which agrees with the assumptions of the kinetic theory.
an algebraic expression of the velocity profile, the so-called Bagnold profile (e.g., Pouliquen and Forterre, 2009): The depth-averaged velocity can be calculated with Molecular dynamic simulations of granular particles on an inclined plane result in a similar velocity profile, yielding an averaged velocity of u ∝ h 1.52±0.05(Silbert et al., 2001).Moreover, this correlation was observed in experiments by Pouliquen (1999).Also, a comparison with velocity profile measurements in real-scale test sites (Kern et al., 2009;Sovilla et al., 2015) shows resemblance in the middle part of the avalanche; see Fig. 12.The obtained velocity profile differs from the plug flow profile assumed by Savage and Hutter (1989).Consistently, the convective flux in the momentum Eq. ( 2) should be corrected with the shape factor (e.g., Baker et al., 2015).However, because of the small influence of this factor (Christen et al., 2010) and potential numerical problems (Hogg and Pritchard, 2004;Baker et al., 2015), the shape factor is set to 1.
Finally a relation between the depth-averaged velocity (Eq.53) and the shear rate at the bottom of the avalanche www.nat-hazards-earth-syst-sci.net/16/2325/2016/Nat.Hazards Earth Syst.Sci., 16, 2325-2345, 2016 γb = γ | z=0 (Eq.51) can be derived: Introducing Eq. ( 54) into the constitutive relation leads to the basal shear stress The factor 5/2 is directly related to the shape of the velocity profile and will change for other profiles, e.g., the velocity profile at the front of the avalanche.Moreover, a plug flow near the free surface is not reproduced by Eq. ( 55) but is visible in the measurement shown in Fig. 12. Dent et al. (1998) show velocity profiles where most shearing is concentrated at the ground.The appearance of a plug flow in measurements can be explained with cohesion (Norem et al., 1987), segregation effects or a flow in transitional state (see also Rauter, 2015, for velocity profiles in transitional states).Expression (55) appears similar to classic friction relations.For example, the relation predicted by Issler and Gauer (2008), who split friction into a quasi-static and dynamic part.In their model, quasi-static stresses are reduced for low concentrations, referring to states with ν < ν rlp as fluidized.Instead of employing a microrheological model, such as kinetic theory, their relation is based on the macroscopic Norem-Irgens-Schieldrop rheology (Norem et al., 1987).Also the similarity to the Voellmy friction relation, given as is remarkable.Because of the similarity and its extensive application in snow avalanche simulations, the Voellmy friction relation is used as a reference to test the performance of the one derived here.For a better comparison to the Voellmy friction parameter ξ the parameter is introduced.This leads to the expression where χ contains the velocity profile dependent factor 5/2.A constant χ indicates that the same shape of the velocity profile in the whole avalanche is assumed.
The difference between the obtained relation and the friction relation of Voellmy is the inverse quadratic dependency on the flow height.This leads to a lower friction for larger flow heights and therefore larger avalanches.This behavior is in line with observations, where the runout is often related to the volume of the avalanche (e.g., Scheidegger, 1973;Bovis and Mears, 1976;Körner, 1980;Alean, 1985).

Model test and parameter evaluation
To test the obtained friction relation, we employ a multivariate optimization method, based on the work of Fischer et al. (2015).This method takes different optimization variables into account, which represent the main avalanche characteristics, e.g., runout or velocity.These can be obtained from simulations and field observations and their residuals can be quantitatively evaluated.Low residuals indicate a good simulation-observation correspondence.The variation of input parameters is limited to the friction parameters, which allows a simple and clear comparison.By scanning the entire physically relevant parameter space, parameter sets, yielding minimal residuals between simulation and observation, are identified.The combination of two different avalanche events, which differ significantly in volume and velocity, is investigated, allowing us to unify parameter sets for avalanches of different types, which is usually not only a superposition of the single events (compare Issler et al., 2005).The two avalanche events are (compare Fischer et al., 2014) avalanche no.103 from 10 February 1999 at the Vallée de la Sionne (VdlS) test site with a deposition volume of approximately 500 000 m 3 and a velocity of up to 70 m s −1 (see Sovilla, 2004;Sovilla et al., 2006, for details) and avalanche from 17 April 1997 at the Ryggfonn (Rgf) test site with a deposition volume of approximately 40 000 m 3 and a velocity of up to 40 m s −1 (see Gauer et al., 2007, for details).
The simulations have been carried out using the SamosAT simulation software, including entrainment and the respective friction relation.To calculate the earth pressure coefficient K a/p , a value of 15 • for the internal friction angle φ is used in all simulations.Note that in SamosAT, φ is set equal to δ, when φ < δ, to solve Eq. ( 4).
Especially for the VdlS avalanche, entrainment appears important because of the high increase of volume during its descent.A simple approach for the entrainment rate q of the form where e b represents the specific erosion energy (compare Fischer et al., 2015) is employed.To estimate appropriate erosion energy coefficients we calculated growth indices, determining the quotient of the deposition mass and the initially released mass.This index is mainly influenced by the entrainment model, the available snow mass and the corresponding parameter.The field observations yield growth indices of 2.3 and 6.0 for Rgf and VdlS, respectively.To resemble values in this range, erosion energy coefficients of 10 4 J kg −1 for Rgf and 10 3 J kg −1 for VdlS were found to be appropriate.The Nat. Hazards Earth Syst.Sci., 16, 2325-2345, 2016 www.nat-hazards-earth-syst-sci.net/16/2325/2016/snow cover height h msc , which is used to limit the entrainment and determine the release volume, was calculated with regards to the elevation and slope inclination α: In order to investigate the range of possible simulation results, the friction parameter µ is varied uniformly between 0.1 and 0.5 for both friction models, the friction parameters ξ and χ are varied between 10 2 and 10 4 (units are m s −2 and m −1 s −2 , respectively).A logarithmic distribution for these parameters is chosen in order to account for the large associated uncertainty, i.e., order of magnitude.
To judge the quality of our simulations, they are compared with measurements of the following three observation variables.
The velocity in the avalanche track obtained by pulsed Doppler radar measurements.The radar measures the radial velocity and in combination with the elevation model, the surface parallel velocity can be calculated.After an appropriate coordinate transformation these values can then be compared with velocities obtained by simulations (Fischer et al., 2014).The radars settings allowed a distance between range gates of 50 m.This leads to a resolution of 14 values along the avalanche path for both events.
The affected area near the deposition area.The deposition area cannot be analyzed directly because the dynamic model does not simulate the deposition process explicitly.Therefore areas where the simulation results exceed a specific dynamic peak pressure, p lim = 1 kPa in our case, are compared.The dynamic pressure is calculated from primary flow variables as with ρ = 200 kg m −3 .Note that the simulation results are independent of the density ρ and the pressure limit p lim may equivalently be expressed in terms of peak velocities.However, defining affected areas and runout in terms of pressure is in accordance with different international hazard mapping guidelines (c.f.Jóhannesson et al., 2009).
The runout distance along the avalanche path.The runout length is measured as projected length in the natural coordinate system, defined by the avalanche track.Just like the affected area, the runout length is defined by the farthest point where the avalanche exceeds the pressure p lim along its cross section (Fischer, 2013).To quantify the quality of a simulation with the parameter set i, we used the residuals between values obtained by the respective simulation X i and the measurements X, calculated as where δX can be the residual of velocity δu or of runout length δr.The residual of the affected area δA is calculated in a similar manner but integrated over the investigated area Here, a i (x, y) denotes whether the pressure exceeded the threshold p lim at the respective position x, y in the simulation i or not: â(x, y) represents the documented affected area in the same manner.Therefore, δA i represents the area where simulation and documentation disagree.In this way we could take into account not only the runout distance from a single point but also the shape of the avalanche.The area where the affected area was analyzed (A oi , area of interest, in Eq. 63) is shown in Fig. 13.It contains the whole runout area of all simulations and the documented affected area.To combine residuals expressed by more than one value (like the velocity in the avalanche track, represented by a value for each range gate) we used a value related to the residual sum of squares of the form The division by the number of values N and taking the square root ensures that the resulting residual is of the same unit and of comparable size with respect to the single values.This eases the interpretation from an engineering point of view.We used the same concept to combine residuals from more events to obtain a single residual which would be obtained by simulating these events with the same parameters.The events in this paper are VdlS and Rgf: To combine residuals of different kinds, like runout and velocity, we normalized the respective residuals with the minimum and maximum residuals from all simulations to eliminate the specific scale: The normalization was always performed after the combination of values of the same kind and after combining two events.This is important because the normalization and combination of residuals is not commutative.This method does not require reference values like an acceptable error or a measurement error.A possible drawback is that larger events have a bigger impact on the results than smaller ones because of the larger absolute values of velocity and runout.If this is not suitable for the respective problem, one could also perform the normalization before combining the events and therefore lay weight on different events equally.The combination of events and measures leads to four possibilities to evaluate and compare model performance with respect to different regards and events (compare Figs. 15,16 and  lest to be fulfilled sufficiently with the simulation results.The last contains the most information ble. he evaluation of 1 600 simulation runs.This number results from two events, two friction models arameters µ and ξ or χ respectively.residuals from all simulations is summarized, separated by event and friction model.  .The last contains the most information lts from two events, two friction models by event and friction model.) The first evaluation is the simplest to be fulfilled sufficiently with the simulation results.The last contains the most information and is therefore the most valuable.
6 Results and Discussion

Simulation results
The following section shows the evaluation of 1600 simulation runs.This number results from two events, two friction models and 20 values for the friction parameters µ and ξ or χ, respectively.In Fig. 15 the evaluation of residuals from all simulations is summarized, separated by event and friction model.Figure 16 shows the same result combined for both events.Additionally the combined residuals in dependence of the respective friction parameters are highlighted.

Predictive power of employed friction relations
To gain a first insight concerning the validity of the backcalculated model parameters for the two investigated friction relations, we determine the optimal parameter set for the Rgf event and obtain µ = 0.437, ξ = 10 000 m s −2 for the Voellmy relation and µ = 0.395, χ = 7848 m −1 s −2 for the kinetic theory relation.Prediction of the VdlS event with these parameters yields a global residual of 0.27 (δr = 300 m, δA = 140 000 m 2 , δ u = 18 m s −1 ) for the Voellmy relation and 0.15 (δr = 180 m, δA = 98 000 m 2 , δ u = 13 m s −1 ) for the kinetic theory relation.Switching the events and performing the same procedure (i.e., utilizing the parameters optimized for the VdlS event for a Rgf prediction) leads to a global residual of 0.26 (δr = 120 m, δA = 35 000 m 2 , δ u = 9 m s −1 ) for the Voellmy relation and 0.19 (δr = 110 m, δA = 24 000 m 2 , δ u = 9 m s −1 ) for the kinetic theory relation.These residuals indicate an enhanced prediction with the presented friction relation in comparison to the Voellmy reference for both cases.

Single and combined parameter optimization
The runout distance represents a point in the avalanche path.Simulations with high friction (high values for µ, low values for ξ and χ ) may not reach this point, simulations with low friction may exceed this point.Between those limiting cases, a simulation with optimal parameters fits the documented runout almost perfectly.In order to satisfy two observations from two events, each simulation demands its own optimal parameter set (Fig. 15).Therefore, simulations with the same parameter set for two events will lead to a residual between predicted and documented runout (Fig. 16).The residual for the runout length is 22.9 m for the kinetic theory and 27.4 m for the Voellmy relation at best when optimizing both events together.
The affected area is also a measure related to runout.However, it provides an additional important information on the lateral extend and spatial distribution of the avalanche.The correlation to runout is clearly visible in Figs. 15 and 16, as the respective areas of low residuals overlap.Figure 13 shows the documented affected area alongside with the affected area obtained from simulations with the smallest residuals in this respect (δA min ).A perfect correspondence between the documented area and the affected area in the simulation does not appear in any of our simulation runs.This can also be seen in Fig. 15.The smallest residual is about 7700 m 2 for Rgf and 65 000 m 2 for VdlS.
In case of the Rgf avalanche, it is observed that the agreement of documented and simulated affected area is limited.This can be attributed to a large amount of deposited snow in the runout, which is not considered in the digital elevation model, leading to an upstream spreading of the avalanche (see Fig. 13a).All simulations are affected equally by this effect, which leads to the big red areas in Fig. 15a and b.
For the VdlS event, the delineation of the documented affected area is accompanied by high uncertainties due to the large powder cloud of this avalanche.The applied documented affected area represents areas with clearly visible snow depth variations (deposition and erosion) caused by the avalanche (Vallet et al., 2001).Figure 13b shows that simula-  tions with an overall good accordance in the runout fail to reproduce the high climb on the counter slope (the two humps opposite to the two main avalanche tracks).This is the exact area where one expects the powder snow layer to detach from the dense flow layer (the dense flow layer follows the terrain more strictly than the powder snow layer).
Another interesting detail can be observed in Fig. 13b.In contradiction to the Voellmy relation, the kinetic theory relation predicts a separation of two branches in the runout zone, which matches the observed behavior.This is an indication for a proper description of important physical processes in the avalanche.
www.nat-hazards-earth-syst-sci.net/16/2325/2016/Nat.Hazards Earth Syst.Sci., 16, 2325-2345, 2016 The velocity along the radar path is visualized in Fig. 14.The residuals are shown in Fig. 15.For the Rgf event the smallest residual among all simulations is 2.4 m s −1 and for the event in the VdlS the smallest residual is 6.2 m s −1 .
In Fig. 14a the dashed yellow line, representing the Voellmy simulation with the smallest residual in velocity for both events combined, matches the observed velocity quite well.However, the runout prediction of the respective simulation appears to be inaccurate.This example shows the importance of the evaluation of different observation variables.
The dynamic pressure can be calculated from velocity with Eq. ( 61).The residual follows from Eq. ( 62) as for other observation variables.A direct calculation of the residual in pressure using the residual in velocity is not possible because of their nonlinear correlation.Minimal residuals in velocity yield minimal residuals in pressure of 28 and 118 kPa for Rgf and VdlS, respectively.Minimal residuals differ only slightly between the two investigated friction relations, when optimizing on single events.However, when using the same set of material parameters for both events, the residual gained with the kinetic theory model is 88 kPa compared to 116 kPa from calculation with the Voellmy model (about 25 %).

Combinations and parameter matching
Possible best-fit parameters can be obtained by analyzing the overlapping areas in Fig. 15.This areas represent parameter combinations which yield relatively small residuals.Figure 16a and b show the same kind of areas for a combination of the two events.The bigger influence of the event in VdlS is clearly visible, as Fig. 16a and b are quite similar to Fig. 15c  and d.The white circles in Fig. 16 mark the positions of the smallest combined residual.In both cases it is located in the overlapping area of relatively small individual residuals.
The combined residual of velocity, runout length and affected area is shown in Fig. 16c and d for both friction relations.The position of the minimal combined residuals in the parameter space is marked with white circles.The form of the isolines in Fig. 16c and d matches qualitatively the overlapping areas of small residuals in Fig. 16a and b.The combined normalized residual as calculated here seems an appropriate method for the determination of optimized parameter sets.
The parameters for the kinetic theory model, which yield the smallest combined residual for both events (e.g., residual in velocity, runout and affected area, white circles in Fig. 16b and d), are µ = 0.311, χ = 3793 m −1 s −2 .For the Voellmy model, these parameters are µ = 0.268, ξ = 2976 m s −2 (white circles in Fig. 16a and c).
For comparison, practical guidelines propose the variation of Voellmy friction parameters with different avalanche characteristics; e.g., µ is varied between 0.155 for big avalanches and 0.3 for small avalanches (e.g., Salm et al., 1990;Gruber et al., 1999).This tendency is in accordance to the presented results (µ = 0.437 for Rgf and µ = 0.268 for VdlS).How-ever, values for µ found in this work are higher compared to the ones utilized in practice.For the turbulent friction parameter ξ , guidelines propose a relation to terrain features and suggest values between 500 m s −2 (high roughness, channeled avalanche path) and 1000 m s −2 (small roughness).It appears that the newly obtained values for the turbulent friction coefficient ξ are about 1 magnitude larger, which is in accordance to previous observations (Fischer et al., 2015, and references therein).However a direct comparison between single model parameters is limited, since different model implementations (e.g., entrainment) lead to different optimal parameter settings.

Back calculation of microscopic material parameters
Since a direct measurement is hardly possible, optimized values for the parameters µ and χ can be combined with observable properties to estimate microscopic parameters of the extended kinetic theory in the spirit of Salm (1993).
The tangent of the critical friction angle tan φ ss is equal to the Coulomb friction coefficient µ.Moreover, tan φ ss can be used to estimate limit values for the concentration ν.The correlation between tan φ ss and ν s has been investigated by Chialvo et al. (2012) using 3D DEM simulations.For tan φ ss = 0.405, matching approximately the value obtained from the real-scale avalanche simulations (Rgf), Chialvo et al. (2012) obtain ν s = 0.581.With a similar approach, Silbert (2010) obtains ν rlp = 0.556 for tan φ ss = 0.405.The particle density ρ p is assumed to be 300 kg m −3 , following field measurements by McClung and Schaerer (1985) and experiments by Steinkogler et al. (2015a).This leads to a maximal bulk density of 174.3 kg m −3 , which is reasonable for snow avalanches (e.g., Gauer et al., 2008) and close to the assumed value of 200 kg m −3 employed in this work to estimate the dynamic pressure.The Young's modulus E is estimated using the empirical relation of Scapozza (2004), where ρ p is the particle density with unit kg m −3 and E has the unit Pa.For ρ p = 300 kg m −3 , E = 1.64 × 10 7 Pa and The coefficient of restitution is assumed to be 0.1, following Steinkogler et al. (2015a).The sensitivity on the parameter c is rather small and the value of 0.5 is not changed.The macroscopic friction parameter λ is calculated from χ = 7848 m −1 s −2 (Rgf), using the correlation given by the velocity profile, Eq. ( 57), leading to λ = 0.04 Pa s 2 .The unknown microscopical material parameters d and a can be estimated by fitting Eq. ( 44) to λ = 0.04 Pa s 2 for a specific range of σ = [10 2 Pa; 10 5 Pa] and γ = [10 0 s −1 ; 10 4 s −1 ] or until ν < ν rlp .This leads to a particle diameter d = 0.02 m and a dimensionless coefficient for the critical state line a = 0. .Areas in the parameter space with relatively small residuals (less than 10 % on the normalized scale) for the runout length (blue), the affected area (red) and the velocity (yellow) for the different events and friction models.The triangle in the respective color marks the simulation with the smallest residual δX min .These residuals correspond to the evaluations of a single variable of a single event (type a).The white triangle marks the smallest combined residual, which corresponds to an evaluation of combined residuals of a single event (type b).
as always performed after the combination of values of the same kind and after combining two events.
cause the normalization and combination of residuals is not commutative.
not require reference values like an acceptable error or a measurement error.A possible drawback is that bigger impact on the results than smaller ones because of the larger absolute values of velocity and runout.
e for the respective problem, one could also perform the normalization before combining the events and on different events equally.The combination of events and measures leads to four possibilities to evaluate performance with respect to different regards and events (compare Figs. 15, 16 and Table 2): ent with respect to all investigated observation variables (δr ∧ δA ∧ δu, marked by ) s with respect to a single observation variable (δr VdlS+Rgf , δA VdlS+Rgf and δu VdlS+Rgf , marked by ) s with respect to all investigated observation variables (δr VdlS+Rgf ∧ δA VdlS+Rgf ∧ δu VdlS+Rgf , marked by ) is the simplest to be fulfilled sufficiently with the simulation results.The last contains the most information most valuable.cussion sults n shows the evaluation of 1 600 simulation runs.This number results from two events, two friction models e friction parameters µ and ξ or χ respectively.luation of residuals from all simulations is summarized, separated by event and friction model.
The normalization was always performed after the combination of values of the same kind and after combining two events.This is important because the normalization and combination of residuals is not commutative.
This method does not require reference values like an acceptable error or a measurement error.A possible drawback is that 5 larger events have a bigger impact on the results than smaller ones because of the larger absolute values of velocity and runout.
If this is not suitable for the respective problem, one could also perform the normalization before combining the events and therefore lay weight on different events equally.The combination of events and measures leads to four possibilities to evaluate and compare model performance with respect to different regards and events (compare Figs. 15,16 and  The first evaluation is the simplest to be fulfilled sufficiently with the simulation results.The last contains the most information and is therefore the most valuable.
15 6 Results and Discussion

Simulation results
The following section shows the evaluation of 1 600 simulation runs.This number results from two events, two friction models and 20 values for the friction parameters µ and ξ or χ respectively.
In Fig. 15 the evaluation of residuals from all simulations is summarized, separated by event and friction model.Figure 16 20 shows the same result combined for both events.Additionally the combined residuals in dependence of the respective friction parameters are highlighted.

Predictive power of employed friction relations
To gain a first insight concerning the validity of the back-calculated model parameters for the two investigated friction relations, we determine the optimal parameter set for the Rgf event and obtain: µ = 0.437, ξ = 10 000 m s −2 for the Voellmy as always performed after the combination of values of the same kind and after combining two events.
cause the normalization and combination of residuals is not commutative.
not require reference values like an acceptable error or a measurement error.A possible drawback is that bigger impact on the results than smaller ones because of the larger absolute values of velocity and runout.
e for the respective problem, one could also perform the normalization before combining the events and on different events equally.The combination of events and measures leads to four possibilities to evaluate performance with respect to different regards and events (compare Figs. 15,16 and  is the simplest to be fulfilled sufficiently with the simulation results.The last contains the most information most valuable.cussion sults n shows the evaluation of 1 600 simulation runs.This number results from two events, two friction models e friction parameters µ and ξ or χ respectively.luation of residuals from all simulations is summarized, separated by event and friction model.
The normalization was always performed after the combination of values of the same kind and after combining two events.This is important because the normalization and combination of residuals is not commutative.
This method does not require reference values like an acceptable error or a measurement error.A possible drawback is that 5 larger events have a bigger impact on the results than smaller ones because of the larger absolute values of velocity and runout.
If this is not suitable for the respective problem, one could also perform the normalization before combining the events and therefore lay weight on different events equally.The combination of events and measures leads to four possibilities to evaluate and compare model performance with respect to different regards and events (compare Figs. 15,16 and  The first evaluation is the simplest to be fulfilled sufficiently with the simulation results.The last contains the most information and is therefore the most valuable.
15 6 Results and Discussion

Simulation results
The following section shows the evaluation of 1 600 simulation runs.This number results from two events, two friction models and 20 values for the friction parameters µ and ξ or χ respectively.
In Fig. 15 the evaluation of residuals from all simulations is summarized, separated by event and friction model.Figure 16 20 shows the same result combined for both events.Additionally the combined residuals in dependence of the respective friction parameters are highlighted.

Predictive power of employed friction relations
To gain a first insight concerning the validity of the back-calculated model parameters for the two investigated friction relations, we determine the optimal parameter set for the Rgf event and obtain: µ = 0.437, ξ = 10 000 m s −2 for the Voellmy  ciently with the simulation results.The last contains the most information simulation runs.This number results from two events, two friction models χ respectively.
ulations is summarized, separated by event and friction model.It appears reasonable that the particle diameter significantly varies in flow and deposition and increases during descent, since the temperature is rising (Steinkogler et al., 2015b;Vera Valero et al., 2015).Consequently, one expects to find larger particles in the deposition area as in the flowing avalanche.In snow chute experiments Rognon et al. (2008) mostly observed particle diameters on a centimeter scale following a power law and Tiefenbacher and Kern (2004) used a sub-centimeter setup to resolve snow particle and aggregate motion.Considering these experimental observations and earlier theoretical assumptions (e.g., Salm, 1993), a particle diameter of 0.02 m seems reasonable during avalanche motion.

Conclusions
From Fig. 15, one can see that both rheological models can be fitted almost equally well to single observations from field tests.The smallest residuals differ only slightly for the single cases.When different observation variables are combined, the kinetic theory approach allows a better fit to the observed data.This is indicated by the larger overlapping area of the three relatively small residuals in Fig. 16a compared to the areas in Fig. 16b.It stands to reason that the modification of the friction with the flow height can help to represent different stages or flow regimes of the avalanche better.This leads to a more realistic dynamic description in different parts of the avalanche, namely the avalanche track and the runout area.This tendency increases with the number of observations combined.Table 2 shows an overview over possible evaluations and values obtained for both investigated models, where this trend is clearly visible.The difference between the Voellmy relation and the kinetic theory approach increases with the number of combinations in the optimization process.Figure 16 shows residuals of the combination of events.The smallest combined residual for simulations with the Voellmy relation is 0.057 (combination of δr = 64 m, δA = 58 000 m 2 , δu = 8 m s −1 and δp = 159 kPa).The kinetic theory approach reduces this value to 0.020 (combination of δr = 24 m, δA = 50 000 m 2 , δu = 7 m s −1 and δp = 132 kPa).This corresponds to a reduction of the residual in runout by about 60 %, alongside with a reduction of the residual in the pressure along the avalanche track by about 20 %.
This improvement can be obtained with very little modification to current models and simulation tools; i.e., no additional transport equation needs to be solved.Convection and diffusion are neglected and a local equilibrium of the granular temperature is assumed to get an analytical expression for the shear stress.An additional improvement with a more accurate description of the velocity profile is expected.A more realistic velocity profile should also lead to different friction Nat.Hazards Earth Syst.Sci., 16,2016 www.nat-hazards-earth-syst-sci.net/16/2325/2016/ in head and tail of the avalanche like proposed by Buser and Bartelt (2009).
Overall, velocities predicted by the presented models can match the observations quite well with an optimized set of parameters.However, this may also be attributed to the considered entrainment process since the analysis of similar friction approaches showed less agreement of the velocities, disregarding entrainment (Fischer et al., 2014).This highlights the importance of considering friction and entrainment equally in a process-orientated approach and the respective impact on avalanche velocities along the track.
The evaluation of the affected area is accompanied by large uncertainties in the documentations.Therefore, assumptions about the quality of the model in this regard are limited.
Another problem of observations in the runout zone is the rising temperature of the avalanche with its descent.The temperature increases because of dissipating kinetic energy and entrainment of warm snow (e.g., Vera Valero et al., 2015).This is not considered in the presented rheological model.The assumptions of the applied model, e.g., incompressibility and equilibrium of the granular temperature in the dense flow, are encouraged by the obtained results.Applying the microrheological parameters for snow to the mean free path between collisions, we obtain s = 2 × 10 −4 m.In combination with the low coefficient of restitution = 0.1 the granular temperature will dissipate rapidly when no additional energy is supplied, leading quickly to the assumed equilibrium.Moreover, the difference between limit states of density in the dense flow, calculated using ν rlp and ν s found for snow, is about 4 %, which indicates an almost incompressible flow.The validity of a steady state velocity profile cannot be verified.However, applying prescribed velocity and pressure profiles is the core feature of all depth-averaged models.

Outlook
For future works we suggest to include the energy conservation into the dynamic flow model to track the evolution of the thermodynamic temperature.Consequently, the particle diameter can be calculated (e.g., with the relation proposed by Steinkogler et al., 2015a) and used to model the dynamic friction with respect to the diameter dependency as predicted by the extended kinetic theory, (70) This will lead to a significantly increased dynamic friction for wet snow avalanches.Assuming a particle diameter of 0.2 m for wet snow (Bartelt and McArdell, 2009;Steinkogler et al., 2015a), χ will decrease down to 80 m −1 s −1 .The possibility of a negative coefficient χ at high pressures, as proposed by the extended kinetic theory, should be further investigated.This effect cannot be described by the simplified relation and does therefore not appear in the extended kinetic theory relation with the back-calculated microscopic parameters.However, this effect may help explaining the low friction of catastrophic ice and snow avalanches (e.g., Alean, 1985).
The extended kinetic theory predicts a relation between classic flow variables (e.g., σ , γ ) and the bulk density ρ = ν ρ p , matching the critical state line for small shear rates.Although the dense flow regime is nearly incompressible (e.g., for the snow parameters obtained in this work, the limit states for bulk density in the dense flow regime differ by about 4 %), this relation can be exploited to extended the Savage-Hutter model to compressible flows, based on microrheological effects.
Moreover, the transition point between dense flow regime and purely collisional regime, where the density decreases rapidly is a meaningful result for the simulation of powder snow avalanches.We found a simplified empirical description for the transition point, when applying the microscopic material parameters for snow: γ ≈ 3.16 σ 0.368 . (71) This relation can be used to develop novel coupling models of dense flow and powder cloud.There are several other models simulating the coupling between dense flow and powder cloud (e.g., Sampl and Zwinger, 2004;Bartelt et al., 2016).

Summary
In summary this paper highlights the application of a rheological model based on kinetic theory to depth-averaged www.nat-hazards-earth-syst-sci.net/16/2325/2016/Nat.Hazards Earth Syst.Sci., 16, 2325-2345, 2016 snow avalanche simulations.To combine both frameworks we employed the commonly accepted assumption of a constant velocity profile along the avalanche and during its decent.The resulting relation shows similarities to classic friction relations.The employed comparison method allowed to evaluate the different basal friction models with respect to different observation variables.Here the residual sum of squares in combination with a normalization, such that values with different physical units and orders of magnitude can be combined, allowed the comparison of the presented friction relation to the widespread Voellmy friction relation.
Utilizing the new relation shows some improvements, particularly when evaluating different observation variables and multiple events.
10 Data availability The underlying data are intellectual property of the BFW and its scientific partners and are not available to the public.For scientific collaboration and data usage, interested researchers are invited to get in contact with the authors.

Figure 1 .
Figure 1.Simple shear test: C denotes the center of gravity with an arrow showing the direction of its displacement during shearing.

Figure 3 .Figure 4 .
Figure 3. Critical state surface in the σγ -ν space.The color marks the origin of stresses: in yellow areas frictional stresses are dominant and in blue areas collisional stresses.

Figure 9 .
Figure 9. Concentration ν (top) and stress ratio τ/σ (bottom) as proposed by the extended kinetic theory (colored solid lines) in comparison to the simplified relation (black dashed lines in bottom figure)for the idealized quartz sand.The simplified relation for the stress ratio approximates the complex extended kinetic theory good for concentrations between ν rlp and ν * , e.g., in the monotonic dense flow regime (blue area, top-right).Between ν * and ν s the regime is dense but non-monotonic (red area, top-right), below ν rlp the regime is purely collisional (yellow area, top-right) and the concentration decreases rapidly with increasing shear rate.

Figure 12 .
Figure11.Orientation of the coordinate system and stresses in the slope parallel direction on an infinitesimal small control volume.Slope parallel normal stresses (K a/p σ ) cancel each other and are not shown.

where h 0
represents the snow cover height at the elevation z ref and h represents its increase with elevation.This approach ensures a smooth initial snow distribution.The snow cover parameters (z ref = 2400 m, h 0 = 1 m, h = 10 −4 for VdlS and z ref = 1500 m, h 0 = 2 m, h = 6 × 10 −5 for Rgf) are chosen to match field observations of release volume and snow depth estimates.
Figure 16 ed for both events.Additionally the combined residuals in dependence of the respective friction ployed friction relations ning the validity of the back-calculated model parameters for the two investigated friction reimal parameter set for the Rgf event and obtain: µ = 0.437, ξ = 10 000 m s −2 for the Voellmy 848 m −1 s −2 for the kinetic theory relation.Prediction of the VdlS event with these parameters (δr = 300 m, δA = 140 000 m 2 , δ u = 18 m s −1 ) for the Voellmy relation and 0.15 (δr = 180 m, 1 ) for the kinetic theory relation.Switching the events and performing the same procedure (i.e.21 spect to a single observation varia ) b.
Figure 16 in dependence of the respective friction rs for the two investigated friction re-437, ξ = 10 000 m s −2 for the Voellmy f the VdlS event with these parameters Voellmy relation and 0.15 (δr = 180 m, nd performing the same procedure (i.e.

Figure 13 .
Figure 13.Outlines of the numerical simulations in comparison with the documented affected area for the avalanche event in Ryggfonn (a) and Vallée de la Sionne (b).The red and blue lines show the p lim isobars for the simulation with the smallest residuals in affected area (δA min ) for the Voellmy friction model and the new friction model, respectively.The yellow filled areas show the documented affected areas and the orange filled areas show the release areas.The evaluation of the affected area was limited to the area within the black polygons.In the figure showing the event in Rgf one can see that the avalanche (yellow area) stopped and spread apart about 50 m before the dam.This results from a large snow deposit uphill of the dam.Because the digital elevation model does not take into account the snow numerical simulations show the same behavior of the avalanche with an offset of about 50 m.The smallest residual is achieved by simulations with high friction, which leads to a stopping of the avalanche before reaching the dam and an untypical form.

Figure 14 .
Figure 14.Velocity measurements compared with simulation outcomes for the avalanche event in Ryggfonn (a) and Vallée de la Sionne (b).The flow direction is right to left.The x axis shows the distance from the radar station and the left y axis shows the velocity.For a better orientation the elevation profile is shown in black and the right y axis shows the elevation.The red line shows the velocity obtained by the pulsed Doppler radar measurements with an estimated observational error.The yellow and blue lines show the velocity along the radar path in the simulation with the smallest residual for the Voellmy friction model and the new friction model, respectively.The dashed lines show the best simulation when using the same material parameter for both events.For the kinetic theory friction model the optimized parameter set for VdlS and both events combined coincide.The background shows the distribution of velocities obtained by all simulations.
Figure15.Areas in the parameter space with relatively small residuals (less than 10 % on the normalized scale) for the runout length (blue), the affected area (red) and the velocity (yellow) for the different events and friction models.The triangle in the respective color marks the simulation with the smallest residual δX min .These residuals correspond to the evaluations of a single variable of a single event (type a).The white triangle marks the smallest combined residual, which corresponds to an evaluation of combined residuals of a single event (type b).
Figure 16 lt combined for both events.Additionally the combined residuals in dependence of the respective friction lighted.er of employed friction relations ht concerning the validity of the back-calculated model parameters for the two investigated friction ree the optimal parameter set for the Rgf event and obtain: µ = 0.437, ξ = 10 000 m s −2 for the Voellmy 95, χ = 7 848 m −1 s −2 for the kinetic theory relation.Prediction of the VdlS event with these parameters al of 0.27 (δr = 300 m, δA = 140 000 m 2 , δ u = 18 m s −1 ) for the Voellmy relation and 0.15 (δr = 180 m, = 13 m s −1 ) for the kinetic theory relation.Switching the events and performing the same procedure (i.e.
Figure 16 lt combined for both events.Additionally the combined residuals in dependence of the respective friction lighted.er of employed friction relations ht concerning the validity of the back-calculated model parameters for the two investigated friction ree the optimal parameter set for the Rgf event and obtain: µ = 0.437, ξ = 10 000 m s −2 for the Voellmy 95, χ = 7 848 m −1 s −2 for the kinetic theory relation.Prediction of the VdlS event with these parameters al of 0.27 (δr = 300 m, δA = 140 000 m 2 , δ u = 18 m s −1 ) for the Voellmy relation and 0.15 (δr = 180 m, = 13 m s −1) for the kinetic theory relation.Switching the events and performing the same procedure (i.e.

25
relation and µ = 0.395, χ = 7 848 m −1 s −2 for the kinetic theory relation.Prediction of the VdlS event with these parameters yield a global residual of 0.27 (δr = 300 m, δA = 140 000 m 2 , δ u = 18 m s −1 ) for the Voellmy relation and 0.15 (δr = 180 m, δA = 98 000 m 2 , δ u = 13 m s −1 ) for the kinetic theory relation.Switching the events and performing the same procedure (i.e.21 10 (a) To a single event with respect to a single observation varia Figure 16 ditionally the combined residuals in dependence of the respective friction ons e back-calculated model parameters for the two investigated friction rer the Rgf event and obtain: µ = 0.437, ξ = 10 000 m s −2 for the Voellmy inetic theory relation.Prediction of the VdlS event with these parameters 40 000 m 2 , δ u = 18 m s −1 ) for the Voellmy relation and 0.15 (δr = 180 m, ry relation.Switching the events and performing the same procedure (i.e. are shown in Fig. 17 . The full and simplified relation are in agreement for the dense flow regime, i.e., ν > ν rlp .At ν = ν rlp the stress ratio shows a sharp bend and the friction is overestimated by the simplified relation for ν < ν rlp .The value for λ, obtained from the real-scale simulations, cannot be reconstructed with any particle diameter significantly higher then 0.02 m because λ reacts much more sensitively to changes of d than to any other parameter modification.Particle diameters in the deposition of avalanches vary between 0.035 m (smaller particle aggregates) and 0.335 m (large snow boulders)(Bartelt and McArdell, 2009).Steinkogler et al. (2015a) obtain diameters between 10 −3 m (T < −1 • C), 0.02-0.1 m (−1 • C < T < 0 • C) and 0.1-0.25 m (T > 0 • C) in DEM simulations and laboratory experiments.
Figure 16.(a, b)  Areas in the parameter space with relatively small residuals (less than 10 % on the normalized scale) for the runout length (blue), the affected area (red) and the (yellow) both events combined.The circle in the respective color marks the with the smallest residual δX min (optimization to a single variable both events, c).The white circle shows the minimal combined residual of runout, affected area and velocity (optimization, type d).(c, d) The combined residual in the parameter space.This surface has a clearly visible local minimum within the physical relevant area.Their position matches with the white circles in the graphs above.The minimal residual is in both cases located within the intersecting areas.
Simple shear flow setup with linear velocity profile u(z).Flow variables included into the constitutive relation, e.g., the shear rate γ , normal stress σ and the shear stress τ , are consequently constant.

Table 1 .
Material parameters of the extended kinetic theory.

Table 2 .
Obtained residuals for all possible result evaluations.Each symbol marks a set of parameters which can also be seen in the parameter space in Figs. 15 and 16.