Computational Snow Avalanche Simulation in Forested Terrain

Two-dimensional avalanche simulation software operating in three-dimensional terrain is widely used for hazard zoning and engineering to predict runout distances and impact pressures of snow avalanche events. Mountain forests are an effective biological protection measure against avalanches; however, the protective capacity of forests to decelerate or even to stop avalanches that start within forested areas or directly above the treeline is seldom considered in this context. In particular, runout distances of small-to medium-scale avalanches are strongly influenced by the structural conditions of forests in the avalanche path. We present an evaluation and operationalization of a novel de-trainment function implemented in the avalanche simulation software RAMMS for avalanche simulation in forested terrain. The new approach accounts for the effect of forests in the avalanche path by detraining mass, which leads to a deceleration and runout shortening of avalanches. The relationship is parameterized by the detrainment coefficient K [kg m −1 s −2 ] accounting for differing forest characteristics. We varied K when simulating 40 well-documented small-to medium-scale avalanches, which were released in and ran through forests of the Swiss Alps. Analyzing and comparing observed and simulated runout distances statistically revealed values for K suitable to simulate the combined influence of four forest characteristics on avalanche runout: forest type, crown closure, vertical structure and surface cover, for example, values for K were higher for dense spruce and mixed spruce-beech forests compared to open larch forests at the upper treeline. Considering forest structural conditions within avalanche simulations will improve current applications for avalanche simulation tools in mountain forest and natural hazard management.


Introduction
Avalanche dynamics models are widely used for hazard zoning and engineering to predict runout distances and impact pressures of snow avalanche events (Gruber and Margreth, 2001;Ancey et al., 2003;Gruber and Bartelt, 2007).The effect of mountain forests as an effective biological protection measure against avalanches has rarely been addressed in this context (Berger and Rey, 2004;Gruber and Bartelt, 2007;Teich and Bebi, 2009).Large destructive avalanches, 1 which often destroy the forest without a significant deceleration, are of major interest in hazard zoning (e.g., Gruber and Häfner, 1995;Fuchs et al., 2005).Yet, frequent smallto medium-scale avalanches are also often a threat to roads, railways and ski runs below the forest (Techel et al., 2013;Teich et al., 2013).Especially when it comes to decisions about the size and extent of avalanche defense measures (including afforestation) in potential starting zones in forested areas (e.g., in newly created forest openings due to wind disturbance), or directly above the treeline, forest and civil engineers could benefit from reliable avalanche simulation in forested terrain (e.g., Weir, 2002;Schönenberger et al., 2005;Bebi et al., 2009).
The avalanche flow is not only influenced by terrain characteristics, but also by vegetation in the avalanche path (McClung, 2003).A recent study showed that forest structural parameters (e.g., the type of forest and the stem density in avalanche starting zones) have a significant influence on runout distances of small-to medium-scale avalanches starting in forested areas (Teich et al., 2012a).For large avalanches released high above the treeline, this effect is negligible (de Quervain, 1979;Bartelt and Stöckli, 2001;Margreth, 2004;Schneebeli and Bebi, 2004;Christen et al., 2010b).The decreasing speeds and runout distances of largescale avalanches depend mainly on the topography and the distance an avalanche travels through open terrain before penetrating into forests (McClung, 2003;Takeuchi et al., 2011;Anderson and McClung, 2012;Teich et al., 2012a).Both cases have only rarely been implemented in avalanche models (Anderson and McClung, 2012).
Flow models used for avalanche simulation often employ Voellmy-type relations, splitting the total basal friction into a velocity-independent dry-Coulomb term and a velocitydependent "viscous" or "turbulent" friction (Voellmy, 1955).The friction approach has been applied by several authors to model the effect of forest on avalanche runout by increasing friction in forested areas compared to open unforested terrain (Gubler and Rychetnik, 1991;Bartelt and Stöckli, 2001;Gruber and Bartelt, 2007;Teich and Bebi, 2009), but has been verified for few real large-scale avalanche events (Casteller et al., 2008;Takeuchi et al., 2011).Avalancheforest interactions may be only poorly represented within the framework of this model (Teich et al., 2012b).Especially for small-scale avalanches, physical processes within the avalanche flow such as snow entrainment (mass uptake) and detrainment (mass extraction) along the avalanche path are important and are not included in the calibrated Voellmy friction coefficients (Maggioni et al., 2012).The local braking effect of forests on avalanche flow seems to be difficult to model with a frictional relationship at the grid scale (Feistl et al., 2014).
Instead of using higher friction values, Feistl et al. (2014) propose an additional one-parameter function (detrainment function) to account for avalanche-forest interactions.Based on field observations, they assume that trees stop fractions of the granular snow flow by a combination of impact, rubbing dissipation, deflection, cohesion and jamming.The stopped snow deposits behind trees, groups of trees or remnant stumps and, therefore, mass is directly extracted from the avalanche volume and the corresponding momentum is removed from the total momentum of the moving snow.This detrainment function accounts for the braking effect of forests on avalanche flow, and can be implemented in numerical avalanche dynamics models.The relationship is parameterized by the detrainment coefficient K, which is related to forest characteristics.Currently, values of K for forested areas have only roughly been estimated and tested for few real avalanche events (Feistl et al., 2014).
Detailed analyses of two-dimensional avalanche simulation software working in three-dimensional terrain objectively require a suitable data selection and a comprehensive and standardized way of processing multiple simulation results (Fischer, 2013).Automatically processing and analyzing large quantities of one-dimensional avalanche model outputs have been conducted in several studies (e.g., Ancey, 2005;Gauer et al., 2009;Eckert et al., 2009).In contrast, multidimensional simulation data have mainly been evaluated manually along predefined cross sections within the avalanche path (e.g., Christen et al., 2010b;Bühler et al., 2011).Manually comparing two-dimensional simulation results with field observations for a high number of avalanche events by visual (image) interpretation is time consuming and rather subjective.To overcome this weakness, a standardized evaluation and comparison method for models operating in three-dimensional terrain has been suggested by Fischer (2013).This approach is employed here to analyze avalanche simulation results automatically and objectively.
In this study, we apply a novel detrainment modeling approach in order to investigate the effect of different forest characteristics on small-to medium-scale avalanches.We compare simulation results of the avalanche simulation software RAMMS (RApid Mass Movement System; Christen et al., 2010a) with runout observations of 40 small-to mediumscale avalanches released in forests of the Swiss Alps in order to operationalize the detrainment function.We evaluate our model by systematically analyzing parameters, characterizing forest structural conditions and their effects on simulated compared to observed runout distances.The overall aim is to define combinations of forest characteristics corresponding to a specific value of the detrainment coefficient K to be applicable in practice.

Avalanche modeling in forested terrain
In this contribution, avalanche flow is modeled using depthaveraged mass and momentum equations; for a detailed mathematical description see Christen et al. (2010a).To briefly summarize: avalanche flow is characterized by unsteady motion with varying flow depth and velocity.Therefore, avalanche flow depth H (x, y, t) and mean avalanche velocity U (x, y, t) are the unknown field variables.The depth-averaged field variables are a function of time (t) and space (x, y) and, thus, the equations to model avalanche flow, i.e., mass balance and momentum equations, are solved from avalanche release (t = 0) to avalanche deposition.
The mass balance in terms of the avalanche flow depth (H ) is given by where Q(x, y, t) denotes the mass source term with Q = Qe + Qd , the sum of the volumetric entrainment Qe and detrainment Qd rates.The mass source term specifies the mass uptake (entrainment) with Q > 0 or mass extraction (detrainment) Q ≤ 0 from the snow cover per unit area as a function of time t; U is the velocity in x and y direction.
The depth-averaged momentum balance is given by and where c x and c y are the velocity profile shape factors, k a/p is the earth pressure coefficient and S f = (S fx , S fy ) T is the total friction (for details on c and k a/p we refer to Christen et al., 2010a).The right-hand side terms of Eqs. ( 2) and ( 3) add up to the driving, gravitational acceleration g in x and y direction: S gx = g x H and S gy = g y H. (4) Avalanche flow resistance is implemented by a "Voellmyfluid" friction relation, assuming small shear strains in the flow body (Salm et al., 1990;Bartelt et al., 1999).The model splits the total basal friction S f into a velocity-independent dry-Coulomb term, which is proportional to the normal stress at the flow bottom (friction coefficient µ) and a velocitydependent "viscous" or "turbulent" friction (friction coefficient ξ ) (Salm, 1993): where g z is the surface normal component of the vector of gravitational acceleration g = (g x , g y , g z ) (see Fig. 1).
|U| is the magnitude of the mean flow velocity given by |U| = U 2 x + U 2 y .Based on observations, we assume that trees in the path of small-to medium-scale avalanches do not break, acting as rigid obstacles causing mass to stop (Faug et al., 2004).The mass removal behind trees starts immediately after the avalanche is released in forests, leading to a significant deceleration and a runout shortening (Teich et al., 2012a).When modeling avalanche flow in forested terrain, we assume that snow detrainment, (i.e., mass removal by trees, remnant stumps or dead wood) is predominant in relation to potential snow entrainment (mass uptake) and, therefore, snow entrainment is neglected, in other words, the mass source term (see Eq. 1) corresponds to Q ≤ 0 (as the sum of the volumetric entrainment rate Qe = 0 and the volumetric detrainment rate Qd < 0).The extracted mass stops promptly and, thus, is instantly subtracted from the avalanche volume (Eq. 1) and the associated momentum of the stopped mass The release area (A r ) as well as forested areas (A f ) have to be defined by the avalanche expert and assigned an appropriate K value dependent on specific forest characteristics which determine the detrainment rate ( Qd ).Avalanche flow in general is modeled by the velocities in x and y direction (U x and U y ) and by the friction S acting in the opposite direction than |U|, and the gravitational acceleration g.
is removed from the total momentum of the avalanche flow (Eqs. 2 and 3).The stopping process is immediate and can be associated with infinite friction.To account for the effect of differing forest conditions on avalanche flow, this relationship is parameterized with the forest detrainment coefficient K [kg m −1 s −2 ] according to with Ṁd as the mass lost by the avalanche in front of treestands.The density of the avalanche snow is denoted with ρ.
Parameter K accounts for the amount of mass detrained by different forest types per unit area and time (Fig. 1).Currently, two approaches exist to model the braking effect of forests on avalanches: -The friction approach increases the turbulent drag of the basal friction S f of the Voellmy-fluid model in forested areas compared to open unforested terrain (Eq.5), e.g., Bartelt and Stöckli (2001), Casteller et al. (2008), Takeuchi et al. (2011).
-The previously summarized detrainment approach, which is based on extracting the mass of snow from the avalanche volume caught behind trees, groups of trees or remnant stumps; for a detailed description see Feistl et al. (2014).
Both the friction and the detrainment modeling approach have the same goal: to explain and quantify the deceleration of the avalanche by the forest (for a detailed discussion see Feistl et al., 2014).Conceptually, the effect of forest on avalanche flow could also be described by an additional drag term interpreted as a retarding force that a tree exerts on the avalanche.The magnitude of such a drag force would be, however, similar to the increase in friction used to model avalanche-forest interactions with a Voellmy-fluid model.Applying the friction approach to small-to medium-scale forest avalanches, utilizing a turbulent drag coefficient ξ of 400 m s −2 for forested areas independent of the forest structure (Gruber and Bartelt, 2007), was not satisfying (Teich et al., 2012b;Feistl et al., 2014).The local retarding effect of forests modeled with the friction approach is significantly larger than the deceleration due to detrainment.On the avalanche path scale, however, this relationship may be reversed.That is, detrainment leads to reduced flow depths as soon as the avalanche penetrates into forests, which reduces the avalanche's potential energy.Furthermore, reduced flow depths also increase the magnitude of the retarding effect associated with the turbulent drag.This effect is especially important for small-to mediumscale avalanches since the ratio of the detrained flow depth to the total flow depth is high, leading to a significant runout shortening.
In this contribution, we evaluate and operationalize the detrainment approach for avalanche simulation in forested terrain.We specifically aim to quantify the effect of the detrainment function (Eq.6) and, therefore, keep the turbulent drag constant, that is, to link the magnitude of this effect to different types of forests growing in the avalanche path.However, a combination of both the friction and the detrainment modeling approaches might finally be suitable to fully account for the effect of forests on avalanche flow.

Avalanche data
Our evaluation and operationalization of the detrainment function were based on 40 small-to medium-scale avalanches released in forests, with runout distances ranging between 50 and 700 m.Within this data set, 38 wet and dry snow avalanches were observed during the winters 1986-1990 in the Swiss Alps (avalanches #1 to #38; Table A1).For these avalanches, the starting points were specified as x, y coordinates and runout distances were recorded from the starting point as the horizontal projection.Detailed data on avalanche characteristics and forest parameters were collected in the field close to the events (Schneebeli and Meyer-Grass, 1993).Since adequately detailed maps of release areas existed only for 26 of these avalanches, we reconstructed the release areas of the remaining 12 avalanches based on given avalanche starting points, maximum release widths, and field notes and photos taken shortly after the avalanche events combined with digital elevation model (DEM) analysis and visual orthophotograph interpretation (Vassella, 2012).In addition, two avalanches (#39 and #40; Table A1), which were released in forests near Davos, Switzerland in the winter 2011/2012, were mapped using a hand-held differential GPS device (for details see Feistl et al., 2014).Forest structural parameters (Table 1), terrain variables and avalanche characteristics such as the type of snow (dry or wet snow avalanche) or the distance an avalanche ran through the forest were assigned to all 40 avalanche events based on collected field data, orthophotograph interpretation and DEM analyses (Table A1).Release heights were measured in the field for 38 observations; for two avalanches (#39 and #40) release heights were estimated based on field visits in combination with measurements of nearby snow and weather stations.Avalanche release volumes (V r ) were calculated corresponding to mapped and reconstructed release areas and release heights.
In this contribution, forests are characterized by a maximum distance between trees of 25 m, a minimum canopy density of 20 %, and a dominant height above 3 m.We chose forest and terrain variables due to pretests of potentially relevant variables and their compatibility with existing assessment methods.Forests were classified into three types, depending on the main tree species: "beech forests", containing beech as well as mixed beech-spruce forests with the main tree species being European beech (Fagus silvatica L.); "spruce forests", consisting of evergreen coniferous forests dominated by Norway spruce (Picea abies (L.) H. Karst.); and "larch forests", defined as deciduous coniferous forests formed by European larch (Larix decidua Mill.) at the upper treeline.Forest density can be characterized by the degree of the crown coverage (Bebi, 1999).Based on the classification system of Bebi et al. (2001), the crown coverage was delineated and digitized in GIS by orthophotographic interpretation and aggregated in three classes, which are described by the variable crown closure (see Table 1).The stage of development indicates the mean stem diameter distribution, which is, for our data set, also represented by its vertical structure (Tables 1 and A2).The terrain variables overall mean slope angle, the cross-slope curvature and terrain roughness were determined from a high-resolution DEM, which was gained from airborne lidar (light detection and ranging) data with a spatial resolution of 2 m and a vertical accuracy of approximately 0.5 m.Cross-slope curvature was categorized as "gully" or concave slope, and "flat" terrain, i.e., almost no curvature; terrain roughness as "low" and "high" (Table 1).For a detailed methodological description, we refer to Teich et al. (2012a).In addition to the terrain roughness gained from the DEM, the small-scale surface roughness was also assigned to each avalanche by the variable surface cover.This variable was mapped in the field and describes the nature of the surface cover.Categories are "smooth", "knobby", "scree" and "stumps/shrubs/saplings" (Table 1).

Simulation software and setup
The detrainment function (Eq.6) was implemented in the avalanche simulation software RAMMS (RAMMS: AVALANCHE version 1.5.01 © WSL/SLF).Based on a two-dimensional depth-averaged flow model (Eqs.1-4), RAMMS calculates the development of avalanche flow depth H (x, y, t) and depth-averaged avalanche velocities U (x, y, t) as a function of time t (see Sect. 2.1); the system of partial differential equations is solved numerically using first-and second-order finite volume techniques (Christen et al., 2010a).The depth-averaged field variables H and U are used to predict avalanche runout distances or impact pressures in complex three-dimensional terrain.Three spatially explicit quantities are required to perform the numerical calculation: (1) a DEM, (2) release areas (A r ), and (3) model friction parameters (µ and ξ , Eq. 5).In addition, to run RAMMS including the detrainment function, forested areas (A f ) have to be defined in the model domain and assigned a K value corresponding to specific forest characteristics.We determined forested areas based on existing forest maps and orthophotographs.In order to focus the evaluation and operationalization on the detrainment function only, snow density was set to ρ = 300 kg m −3 and we kept the friction parameters constant at µ = 0.29 and ξ = 1500 m s −2 throughout this study.We chose this combination since the estimated release volumes of our avalanche data set range between 19 and 3398 m 3 , which corresponds to the avalanche size class "tiny" (< 5000 m 3 ), and is applied in practice to simulate frequent avalanches (10-year return period), in unchanneled terrain above 1500 m a.s.l.(Buser and Frutiger, 1980;Salm et al., 1990).The simulations are based on a DEM with a spatial resolution of 2 m and a vertical accuracy of approximately 0.5 m.The mapped release areas and release heights were used to specify the initial conditions for each simulation run.All simulations were accomplished without any pre-defined stopping criteria.
For each observed avalanche, a reference simulation was computed by running RAMMS without accounting for any forest influence in the avalanche path (K = 0).In order to find optimal values for K dependent on different forest characteristics, we then simulated each observed avalanche with varying values for K of 5,10,20,30,40,60,80,100,130,160,190 and 220 kg m −1 s −2 .These K values were chosen based on results of a computational experiment performed by Feistl et al. (2014).
The main simulation results are maxima over time t of the flow depth H (x, y, t) and the two-dimensional slope parallel velocities U (x, y, t) at a constant density ρ.As usually applied in hazard assessment (e.g., Eckert et al., 2010), the according peak pressure field can then be derived as where x, y denote the two-dimensional Cartesian coordinates.Here U peak corresponds to its maximum U value over the entire simulation time t: For our analyses, we exported the spatially explicit maximum pressure output.

Analyzing simulation results
To compare the two-dimensional model outputs with the one-dimensionally recorded avalanche runout distances, we applied the analysis method AIMEC (Automated Indicator based Model Evaluation and Comparison) presented by Fischer (2013).
The AIMEC approach allows for a standardized and objective evaluation of two-dimensional simulation results.The simulation results are transformed from Cartesian coordinates (x, y) to a coordinate system dependent on the specific avalanche path (s, l) (Fig. 2), here applied for the peak pressure: P (x, y) → P (s, l). (9) As a scalar metric, the runout indicator is defined based on the peak pressure (Eq.7), and evaluated for each simulation run.This runout indicator corresponds to the horizontal projection of length measured along the avalanche path coordinate s, where the cross-sectional maximum peak pressure value falls below a certain pressure limit P max cross (s) < P limit (Fig. 2).We tested pressure thresholds P limit of 1, 3, 5 and 10 kPa as well as 0.5 kPa for very small avalanches with release volumes V r < 100 m 3 .For such small avalanches, the differences between runout indicators determined with P limit = 3 kPa and P limit = 1 kPa for the reference simulations with K = 0 ranged between 1 and 66 % (mean = 22 %).When calculating the difference between both runout indicators for all avalanches of our data set, the mean difference was only 14 % (ranges between 0 and 67 %).For simulations performed with the detrainment function (K > 0), mean differences between the two runout indicators (P limit = 3 kPa and P limit = 1 kPa) decreased for very small avalanches (V r < 100 m 3 ) to 2 % and for all avalanches to 7 %.Due to such small differences, we applied a pressure threshold of P limit = 3 kPa throughout this study, which corresponds to a pressure threshold used for hazard zone mapping in Switzerland (BFF/SLF, 1984).That is, for avalanches with return periods ≤ 30 years an impact pressure > 3 kPa is assigned to have consequences regarding land-use planning (Jóhannesson et al., 2009).In order to measure the differences of simulated runout indicators (runout sim ) to observed runout distances (runout obs ), the relative runout difference ( runout in [%]) is introduced as where positive values indicate overestimated runout distances and negative values for runout reveal that runout distances were underestimated by the avalanche simulation software compared to the recorded ones.

Statistical analysis
For an evaluation of general dependencies between variables describing forest structure, topography and avalanche characteristics, and the response variable runout, we calculated Spearman's rank correlation coefficient (r S ) for categorical and continuous predictor variables, since it is known as nonparametric and does not assume a linear relationship.In addition, Pearson's correlation coefficient (r) was calculated for all continuous variables and runout to reveal potential linear dependencies and to measure their strengths.A correlation was assumed to be statistically significant if the respective p value was 0.01 < p ≤ 0.05 and highly significant for p ≤ 0.01.
Table 2. Significant a (0.01 < p ≤ 0.05) and highly significant b (p ≤ 0.01) Spearman rank correlation coefficients (r S ) between predictor variables and runout ref calculated for the reference simulation runs with K = 0, and between predictor variables and the assigned optimal value for K (K opt ).The evaluation and operationalization of the avalanche model included four steps: 1. We tested all variables (see Table A1) against runout (further referred to as runout ref ) for the reference simulations without any influence of forest (K = 0).
2. Based on the simulations, including the mass extracting effect of forests parameterized with the detrainment coefficient K, we determined an optimal K value for each avalanche event (K opt ).That is, one value for K was defined for each of the 40 avalanche events, which resembled the observed runout distances "best", where K approaches zero of runout, on the condition that runout ≥ 0. A conservative evaluation of simulation results leading to overestimated rather than to underestimated runout distances is preferred to reveal optimal K values which are applicable in practice.
3. We again calculated r S and r, and tested the forest parameters forest type, crown closure, vertical structure, stage of development and surface cover, as well as the release volume and the distance an avalanche ran through forest against the response variable K opt .
4. We defined K values based on specific forest characteristics and their combined effects to be applicable in practice for avalanche simulation in forested terrain.
We evaluated our derived K values by simulating two avalanche events additionally observed in 2012 in forested terrain in the Swiss and Bavarian Alps.These avalanches differed in forest conditions and the distance they ran through forest as well as in the snow type.To further test the practical applicability of the derived K values, we ran RAMMS using a default simulation setup and compared simulation results manually.

Avalanche simulation with K = 0
Runout distances were overestimated by RAMMS for 38 of 40 investigated avalanches in forested terrain when forest influence was not considered.The relative runout difference runout ref (Eq. 11) revealed overestimations by RAMMS up to 700 % for the chosen parameters K = 0, µ = 0.29, and ξ = 1500 m s −2 .The two avalanches with negative values for runout ref (−34 and −48 %) are of very small release volumes (V r < 50 m 3 ).
Variables which affected runout ref of our data set significantly are the release height, the snow type, the absolute as well as the relative distance an avalanche ran through forest, and the surface cover (Table 2).Dependencies between the continuous variables release height and absolute and relative distance through forest are not linear since no significant correlations were found when calculating Pearson's correlation coefficient (r).However, it could be assumed that increasing release heights, accompanied with increasing release volumes (see Table A2), are related to an increase in runout ref .
That is, the bigger an avalanche, the larger the difference between observed and simulated runout distances.Both correlations imply, that a loss of avalanche volume modeled for forested areas may lead to a significant runout shortening and a more realistic avalanche simulation which would match the observations.Differences between observations and simulations were significantly higher for dry snow avalanches compared to wet snow avalanches (Fig. 3).Thus, one can assume that the accompanying snow densities and thermal snow temperatures also determine the detraining effect of forests.Here, snow density was kept constant at ρ = 300 kg m −3 , which is often applied for dry snow avalanches.The snow type was also correlated with release volume and release height (Table A2), where the latter also influenced runout ref significantly (Table 2).The nature of the surface cover was correlated significantly with runout ref .That is, a scree slope and higher small obstacles such as stumps and shrubs in the avalanche path were related to larger differences between observed and simulated runout distances and, therefore, also determine the amount of snow deposited in the avalanche track.
Besides surface cover, distributions of runout ref suggest influences of other forest parameters on avalanche simulations (Fig. 3).In particular, runout indicators for avalanches that started in spruce forests were highly overestimated (median = 88 %, mean = 154 %), but less overestimated for avalanches which ran through beech forests (median = 52 %, mean = 79 %) or larch forests (median = 44 %, mean = 49 %).For simulations without any forest influence, runout ref was largest for avalanches which ran through dense evergreen forests with a more than two-layered vertical Mean slope angle, cross-slope curvature and terrain roughness in terms of local differences in elevation did not influence runout ref significantly.This strengthens the theory that avalanche-forest interactions need to be implemented by a function dependent on forest characteristics in combination with snow conditions only.

Avalanche simulation with varying K values
For the next step of our evaluation and further operationalization, we calculated runout for each simulation run with varying values for K and analyzed relationships between forest characteristics and runout.In general, increasing K values corresponded to decreasing runout indicators, where the strength of this effect seemed to decrease around K values of 150 kg m −1 s −2 and higher (Figs. 4 and 5).Very small avalanches with release volume V r < 100 m 3 showed diverging simulation results.For such avalanches, values of runout were often negative when applying the detrainment function; one avalanche simulation did not even start with the smallest chosen K value of 5 kg m −1 s −2 .However, differences between avalanche simulations in terrain covered with different forest types are visible, especially between larch forests and the two other forest types, spruce and beech forests, when calculating mean values of runout corresponding to each chosen K value for the three categories separately (Fig. 4).In addition, differences in the vertical structure of a forest stand as well as in crown closure had a higher influence on the amount of snow extracted from the avalanche volume compared to a differing stage of development (Fig. 5).The latter forest variable is, however, relatively well represented by the vertical structure (Table A2).The nature of the surface cover also influenced the amount of snow removed from the avalanche volume.The effect of differences in surface cover could have even been underestimated, since our simulation setup did not account for changes in surface cover in unforested areas.
In terms of the operationalization, optimal values for K (K opt ) were assigned to each observed avalanche based on the election rule that runout approaches zero on condition runout ≥ 0. A significant correlation was found between K opt and the forest variable forest type (Fig. 6) as well as for the release volume and the absolute distance an avalanche ran through forest (Table 2); the latter two were even linear with r = 0.35 and p = 0.028 * for release volume, and r = −0.44 and p = 0.005 * * for the distance through forest.Thus, the larger the release volume the higher the K opt , and the longer the distance an avalanche runs through the forest, the lower the corresponding K opt .According to theory, K should only account for forest characteristics.
Thus, we propose choosing a value of K to simulate avalanche runout in forested terrain mainly dependent on the forest type.Based on Fig. 6, possible values for K can be obtained, that is, K values of 5 kg m −1 s −2 may be assigned to areas covered by larch forests, 80 kg m −1 s −2 to forests dominated by spruce, and 100 kg m −1 s −2 to beech and mixed beech-spruce forests.These values should be adapted with K values corresponding to classes of the three forest characteristics crown closure, vertical structure and surface cover (see Fig. 5); the mean value of the respective K values for the four forest characteristics were calculated for our case studies (see Sect. 3.3).The influence of K values higher than approximately 150 kg m −1 s −2 on runout decreases (Fig. 5).Therefore, K values > 150 kg m −1 s −2 do not seem meaningful for modeling avalanche-forest interactions.

Case studies
In order to test the practical application of our results, we simulated two additionally observed avalanches with RAMMS including the detrainment function (Table 3).Therefore, we assigned a K value to forested areas characterized by the forest parameters forest type, crown closure, vertical structure and surface cover as specified in Table 1.
Values of K were estimated based on Figs.4-6.For forest type, crown closure, vertical structure and surface cover, K values close to runout = 0 were chosen and, then, the mean value of K was calculated (Table 3).We ran RAMMS with a default simulation setup, in other words, values for friction parameters µ and ξ were not kept constant but defined by an automatic procedure of RAMMS depending on terrain features such as gullies or flat slopes, elevation, the return period (set to 10 years) and the avalanche size class ("tiny").The simulations were based on a 2 m grid for the avalanche observed in Switzerland, and a 1 m grid for the one from Germany.Forested areas and forest characteristics were delineated based on pixel maps, orthophotographs, and photographs taken during field visits.Again, we ran the simulations until the final pressure patterns were reached.In practice a stopping criteria of 5 % of the total momentum is often applied, indicating that if the sum of all momenta of all grid cells is lower than 5 % of the maximum   momentum sum, the simulation is stopped (Christen et al., 2010a).However, test simulation runs applying this threshold have shown that runout distances of our case studies and, therefore, such small-scale avalanches were highly underestimated.In contrast, we ran our simulations without any stopping criteria and analyzed the simulation results by only displaying the grid cells of the runout area, which exceeded a pressure threshold of 3 kPa.This corresponds to our limit for the maximum peak pressure (P limit ) when defining runout distances by applying AIMEC (see Sect. 2.4) as well as to the impact pressure threshold with consequences for hazard zone mapping in Switzerland (BFF/SLF, 1984; Jóhannesson et al., 2009).
The simulation results showed a good agreement with the observed runout when applying the novel detrainment function (Fig. 7).Even if the runout area did not match the observed ones exactly, runout distances were predicted relatively well by the model for both avalanche events; simulated runout distances stopped within −6 to 3 m compared to the observed ones.

Discussion
In this study, we applied a novel detrainment modeling approach (Feistl et al., 2014) to account for avalanche-forest interactions within computational avalanche simulations.The aim was to evaluate and operationalize the detrainment function (Eq.6) and, therefore, to quantify the detrainment coefficient K, which is associated to the amount of snow caught behind trees in the avalanche path.
In general, immediate stopping and removal of a certain amount of mass by trees has a greater influence on small-to medium-scale avalanches than on larger avalanches (Feistl et al., 2014).Large-scale avalanches are able to break and uproot trees linked to a low energy consumption which increases avalanche mass and, therefore, flow energy (Bartelt and Stöckli, 2001).When applying a Voellmy-type relation, which is often employed by avalanche flow models, the effect of forests on such avalanches can be modeled by increasing friction compared to unforested terrain (Bartelt and Stöckli, 2001).This is not valid for modeling small-scale avalanches in forested terrain: previous simulations of our data set with RAMMS with alternating ξ values for forested areas (100-1000 m s −2 ) showed that runout distances of 31 out of the 40 avalanches were still overestimated when applying the smallest chosen ξ value of 100 m s −2 (Teich et al., 2012b).Moreover, simulating small-scale avalanches with a model based on Voellmy-type frictional relationships only is generally questionable (Sailer et al., 2008), since the avalanche will not stop as long as the slope angle is larger than the friction angle, i.e., tan φ > µ.Therefore, including physical processes within the avalanche flow such as snow entrainment (mass uptake) and detrainment (mass extraction) is important, for example, modeling the mass removal by trees, remnant stumps or dead wood, as realized in this study.We hence kept friction parameters constant to focus the evaluation and operationalization on the detrainment function only.
The results gained from analyzing reference simulations accomplished without any forest influence (K = 0) highlight the importance of modeling the mass loss induced by forests growing in the avalanche path.Significant correlations between the predictor variables release height and the distance an avalanche ran through the forest with the response variable runout ref suggest that a loss of avalanche volume modeled for forested areas will lead to shorter runout distances.In addition, the surface cover in terms of stumps and shrubs or scree slopes also affected runout ref significantly.This effect should also be considered for small-to medium-scale avalanches' simulations in unforested areas, such as large forest openings caused by natural disturbances, which are often interspersed with shrubs, fallen logs, remnant stumps and root plates of upturned trees (Fig. 8).Remaining dead wood is able to increase the surface roughness at least over the first 10-20 years after the dieback (Brown et al., 1998;Rammig et al., 2007).Indeed, the effective heights and interacting avalanche flow depths also determine the mass deposited behind obstacles (Faug et al., 2004;Naaim et al., 2004).Based on sporadic field samples we can assume effective heights of approximately 0-30 cm for "smooth" slopes, 30-50 cm for "knobby" terrain, and 30-150 cm for "scree" slopes as well as for terrain interspersed with stumps, shrubs and/or saplings.The significant correlation between the snow type and runout ref indicates that the effectiveness of the mass removal by forests is also determined by snow densities as well as thermal snow temperatures, for instance, the wetter and more viscous the snow, the slower the avalanche; in future such processes need to be incorporated when modeling small-to medium-scale avalanches (Vera Valero et al., 2012).
In the next step, we simulated each avalanche with varying K values (between 5 and 220 kg m −1 s −2 ) and assigned an optimal value for K (K opt ) to each avalanche event.In general, runout distances decreased with increasing K values while this effect decreased around K = 150 kg m −1 s −2 .However, some of the 40 observed avalanches were still overestimated by RAMMS when simulating with the highest chosen K value of 220 kg m −1 s −2 .On the one hand, partially misinterpreting the orthophotographs and DEMs when reconstructing 12 release areas could have affected the simulation results (Vassella, 2012).On the other hand, other processes such as the influence of thermal snow temperature on the avalanche flow (see above) and the effect of different topographic features (usually modeled by varying friction parameters µ and ξ ) could have also influenced the simulations.In order to reduce uncertainties related to the avalanche modeling process and to account for effects of varying K values on the simulations only, we used constant values for µ and ξ throughout this study (see Sect. 2.3).
The statistical analyses between predictor variables and the response variable K opt revealed that the forest type in which an avalanche was released and ran through had an influence on runout.Thus, the forest type mainly determines the K value to be chosen for avalanche simulation in forested terrain in combination with crown closure, vertical structure, and surface cover since: -clear differences of mean runout between the categories of these forest parameters are visible (Fig. 5) -these variables can be largely derived from remotesensing-based data (orthophotographs, lidar-data) combined with sporadic field samples, but no extensive measurements are required -other studies on the effect of forest structural parameters on observed runout distances emphasize the relevance of these forest characteristics (e.g., McClung, 2003;Teich et al., 2012a).
The case studies performed by simulating two additional avalanches emphasized this argumentation (Table 3 and Fig. 7): the good agreement of the simulated and observed runout distances when applying K values based on the four suggested forest characteristics encourages the applicability of the detrainment function for hazard analyses and, therefore, for a practical natural hazard and protection forest management.
For these two avalanches we applied a default simulation setup and interpreted the simulation results visually, but based on an avalanche pressure threshold of > 3 kPa used for hazard mapping in Switzerland (BFF/SLF, 1984).Impact or peak pressure results are generally of high interest in snow avalanche modeling to estimate the avalanches' destructive potential, and are utilized for hazard zoning and engineering affecting land-use planning in many countries (Jóhannesson et al., 2009).
We also chose the threshold of P limit = 3 kPa when analyzing our simulation results automatically by applying AIMEC (Fischer, 2013).That is, we ran the simulations without any pre-defined stopping criteria such as for the flow momentum or flow depth and used a pressure-based runout indicator to determine simulated runout distances.In the case of very small avalanches, the pressure threshold P limit has to be defined carefully, since pre-defined pressure limits could be too high (i.e., never be exceeded).Defining too-low pressure limits could, however, lead to a misinterpretation of the simulation results, for example, when accounting for runout which is attributed to non-realistic stopping in flat natural terrain due to a diffusive runout behavior arising from the flow model (Fischer, 2013).On the contrary, defining runout distance based on thresholds for the maximum flow momentum or the minimum flow depth could also lead to misinterpretations, which would considerably influence further analyses.The P limit = 3 kPa yielded reliable runout indicators and did not differ considerably from runout indicators determined with lower values.In contrast, a P limit > 3 kPa is not appropriate to determine runout indicators of small-scale avalanches since tested values of 5 and 10 kPa were not exceeded for many simulated avalanches of our data set.A verification of the results received with AIMEC is, however, still necessary since numerical solutions can include singularities, especially when simulating small-scale avalanches.
In this study, we could only compare observed and simulated avalanche runout distances.Reliable observations, as well as measurements and experiments on the effect of forests on the avalanche flow (which also contain more avalanche characteristics such as avalanche velocity and avalanche mass balance), are rare.In addition, more welldocumented avalanches in forested terrain have to be analyzed in the way we did to establish better grounded results on the role of forest type, crown closure, vertical structure and surface cover in avalanche simulation to further improve the new detrainment function, in particular in forested areas with varying decelerating effects.The presented findings are nonetheless a valuable first step to simulating small-to medium-scale avalanches in forested terrain for applicability in hazard analyses.

Conclusions and outlook
The applied detrainment function, which can be implemented in numerical avalanche dynamics models, will improve the simulation of small-to medium-scale avalanches in forested terrain considerably.A value for the detrainment coefficient K can now be defined mainly dependent on the forest parameter forest type in combination with crown closure, vertical structure and surface cover.As the suggested forest characteristics can be largely derived from remotesensing-based data (orthophotographs, lidar data) combined with sporadic field samples, there is a high potential for practical implementations.In addition, we demonstrated that applying a standardized method to analyze a high number of two-dimensional avalanche simulation results automatically increases the reliability of an objective software evaluation; the employed method AIMEC provided accurate runout indicators as the basis for further analyses.
Implementing avalanche-forest interactions in avalanche simulations will facilitate current applications for such software, for instance, by better accounting for the protective effects of forests in natural hazard mapping (Berger and Rey, 2004;Gruber and Bartelt, 2007), for managing mountain forests efficiently (Weir, 2002;Brang et al., 2006;Teich and Bebi, 2009) or to value "avalanche protection by forests" as a key ecosystem service in mountainous regions (Grêt-Regamey et al., 2013).The detrainment function will be implemented in the next version of RAMMS and tested by practitioners based on the findings gained in this study.

Figure 1 .
Figure1.Schematic illustration of avalanche modeling in forested terrain.The release area (A r ) as well as forested areas (A f ) have to be defined by the avalanche expert and assigned an appropriate K value dependent on specific forest characteristics which determine the detrainment rate ( Qd ).Avalanche flow in general is modeled by the velocities in x and y direction (U x and U y ) and by the friction S acting in the opposite direction than |U|, and the gravitational acceleration g.

Figure 2 .
Figure 2. Schematic avalanche simulation result (see Fig.1).Red areas correspond to forests with specific forest characteristics (i.e., tree density illustrated by green dots); displayed is the outline of the peak pressure field with a new coordinate system along the central flow line z(x, y) (in bold).

Figure 3 .
Figure 3. Difference between simulated and observed runout distances ( runout ref ) calculated for the reference simulation runs without any forest influence (K = 0) shown for the subsets of variables snow type and small-scale surface cover which are statistically significant (first row) and for the subsets of four other forest structural parameters (no statistically significant relationships).Boxplots show minimum values, the lower quantile (Q 0.25), the median (Q 0.5), the upper quantile (Q 0.75) and maximum values of runout ref .Points are relative positions of extreme values.

Figure 4 .
Figure 4. Mean values of runout for each applied K value calculated separately for the three forest type categories.The dashed line corresponds to runout = 0 indicating the potential mean optimal K value for each category.

Figure 5 .
Figure 5. Mean values of runout for each applied K value calculated separately for the corresponding categories of four forest variables.The dashed line corresponds to runout = 0 indicating the potential mean optimal K value for each category of the respective forest variable.

Figure 6 .
Figure 6.Optimal K values (K opt ) assigned to each observed avalanche based on simulations with varying values of K shown for subsets of different forest types.Boxplots show minimum values, the lower quantile (Q 0.25), the median (Q 0.5), the upper quantile (Q 0.75) and maximum values of K opt .Point is relative position of extreme value.

Figure 7 .
Figure 7. Simulation results gained with RAMMS including the detrainment function in comparison to the observed runout distances of the two case studies "Dischma valley" (left panel) and "Brecherspitz" (right panel).Contour lines are displayed in 10 m steps.

Figure 8 .
Figure 8. Snow detrained by a stump, highlighting the significant effect of surface cover on small-scale avalanches which should be considered in avalanche simulations.

Table 1 .
Forest parameters and corresponding categories assigned to each avalanche.
* Mean diameter at breast height: outside bark diameter measured 1.37 m above the forest floor on the uphill side of the tree.

Table 3 .
Characteristics and K values corresponding to selected forest parameters of two avalanches which were not included in previous analyses to test the results of the operationalization.