Diachronic quantitative snow avalanche risk assessment as a function of forest cover changes

Abstract

This work proposes a holistic quantitative snow avalanche risk assessment that evaluates, at reasonable computational costs and for various types of buildings, the impact of forest cover changes on the probability distribution of runout distances, impact pressures and subsequent risk estimates. A typical case study of the French Alps shows that, from a completely deforested to a completely forested path, avalanche risk for a building located downslope decreases by 53–99%, depending on how forest cover is accounted for in avalanche statistical–dynamical modeling. Local forest cover data inferred from old maps and photographs further demonstrates that a 20–60% risk reduction actually occurred between 1825 and 2017 at the site because of the afforestation dynamics, with significant modulations according to the considered building technology. These results (1) assert the protective role of forests against snow avalanches, (2) highlight the potential of combining nature-based solutions with traditional structural measures to reduce risk to acceptable levels at reasonable costs, (3) suggest a significant decrease in risk to settlements in areas that encountered similar forest cover changes and (4) open the door to the quantification of long-term avalanche risk changes as a function of changes of all its hazard, vulnerability and exposure drivers in various mountain context.

Introduction

Over the years, classical structural protection measures against natural hazards in mountains (e.g. dams, or additionally reinforced walls) have been increasingly criticized due to their negative environmental and aesthetic impact and, sometimes, their potential to rise exposure to risk by stimulating the development of new settlements, dwellings, infrastructures, etc. (Delage, 2003; Moos and others, 2018). Some studies even consider that traditional protection structures, sometimes referred to as ‘gray infrastructure’, are unable to adapt to changes in hazards driven by climate change (Kumar and others, 2020; Poratelli and others, 2020a). These arguments, in addition to the large construction and maintenance costs of structural protection measures have turned the focus toward more sustainable and cost effective ecosystem-based solutions for disaster risk reduction (Eco-DRR, a part of nature-based solutions for climate change adaptation and disaster risk reduction approaches).

In mountain areas, the most well-known example of Eco-DRR solutions such as forests, protecting people and their assets mostly against gravity-driven natural hazards (Brang and others, 2001, 2006). The protective effects of forest stands against, e.g. rockfalls (Dupire and others, 2016) and snow avalanches (Bebi and others, 2009), is of a great importance in alpine regions, preventing human deaths and destruction of buildings and infrastructures (Getzner and others, 2017). In the case of snow avalanches, the primary effect of forests is to hinder avalanche formation and prevent avalanche release by stabilizing the snowpack in release areas, which would, in terms of forest management, also be the primary target (De Quervain, 1978; Salm, 1978; Gubler and Rychetnik, 1991; Viglietti and others, 2010; Teich and others, 2012a). They also have the capacity to decelerate flowing avalanches (Anderson and McClung, 2012; Takeuchi and others, 2018). Their protective potential is directly linked to forest structural parameters, e.g. the stem density for small-medium avalanches (Teich and others, 2012a, 2014), and the distance traveled before penetrating into the forest for large avalanches triggered above the timberline (McClung, 2003; Takeuchi and others, 2011; Teich and others, 2012a).

Forest–avalanche interaction is an old topic in snow avalanche modeling. In earliest studies, friction was locally increased within the model to mimic the decelerating impact of forests leading to shorter runout distances (Salm, 1978; Gubler and Rychetnik, 1991; Bartelt and Stöckli, 2001; Takeuchi and others, 2011). In such so-called frictional approaches, avalanche models often use a Voellmy friction law that considers the total friction as the sum of a dry-Coulomb coefficient μ and a turbulent term depending on the squared velocity and on the inverse of a coefficient ξ (Voellmy, 1955). The static friction coefficient μ generally varies within the 0.1–0.7 range, and is often thought to summarize snow properties (Salm and others, 1990). The turbulent friction coefficient ξ generally varies within the 1000–10 000 m s⁻² range for forest-free terrain, and aims at representing the roughness of the path potentially related to land cover properties (Salm and others, 1990; Ancey and others, 2003). This explains why, in frictional approaches, more often than not, ξ is lowered to values in the order of 400 m s⁻² in forests. This leads to an increase of the velocity-dependent friction that some authors related to the entrainment of heavy trees and stems by the flow (Bartelt and Stöckli, 2001). However, some research also considered a concomitant slight increase in μ, with a Δμ ranging between +0.02 and +0.05 from forest-free to forested terrain (Gruber and Bartelt, 2007; Christen and others, 2010). More recently, as an alternative to frictional approaches, Feistl and others (2014) proposed the detrainment approach in which the forest–avalanche interaction is modeled through a single-parameter detrainment function. This approach accounts for the braking effect of forests on avalanche flows, and has been recently implemented in most up-to-date snow avalanche simulation models (Védrine and others, 2021).

From a different perspective, research that aimed at better mitigating avalanche risk historically focused on the sole hazard component of the risk, and mostly using deterministic approaches relying on physics (Harbitz and others, 1998). However, even the best numerical avalanche models cannot, without further probabilistic and vulnerability considerations, evaluate the risk levels and related uncertainties that are required, e.g. for land-use planning and the design of defense structures. As a consequence, dealing with snow avalanche risk has recently changed from sole hazard prevention to risk management, which includes explicit consideration of vulnerability and exposure (Bründl and others, 2009). First implementations relied on simple scenarios representing extreme avalanches (Fuchs and others, 2004), but, quickly, extreme-value based and/or statistical–dynamical models representing the full variability of avalanche events likely to occur where put at play (Keylock and others, 1999; Barbolini and others, 2004; Cappabianca and others, 2008; Eckert and others, 2008, 2009; Favier and others, 2014b, 2016). This was clearly a step toward more integrated avalanche risk mitigation solutions potentially accounting, e.g.for acceptability thresholds (Arnalds and others, 2004) and/or cost–benefit constraints faced by stake-holders. However, to our knowledge, few studies so far tried to take into account the protective effects of forests within quantitative avalanche risk assessments, and, when this was attempted, with simple statistical relationships (Grêt-Regamey and Straub, 2006) or scenarios (Teich and Bebi, 2009) as hazard model only. By contrast, comprehensive combinations of numerical probabilistic hazard models taking into account changes in land cover, elements at risk and their vulnerability have already been achieved, e.g. for rockfall risk (Moos and others, 2018; Farvacque and others, 2019) and it is a rather common strategy (although still difficult to properly achieve) for flood risk (Rogger and others, 2017; Bathurst and others, 2020).

That the evolution of avalanche disaster risk has rarely focused on the protective effects of forests is especially surprising in light of the strong reforestation that occurred in the European Alps (and numerous other mountain areas) over the last ~150–200 years (Mather and others, 1999; Bebi and others, 2017). This time frame roughly corresponds to the classical 100–300 year reference periods that define legal thresholds in land-use planning (Eckert and others, 2018). This pleads for a dynamical assessment of the impact of reforestation on avalanche risk accounting for potentially quick changes in forest extension and structure. Already existing examples provide insights about how reforestation can affect avalanche risk. Giacona and others (2018) demonstrated that, in medium-high mountains, interactions between avalanche activity, forest stands, social practices and climate result in strong temporal modulations of mountain landscapes and avalanche risk. García-Hernández and others (2017) highlighted strong decrease in snow avalanche damage in the Asturian massif (Spain) linked to reforestation following the end of industrial activities. Similarly, Mainieri and others (2020) and Zgheib and others (2022) showed a decrease in avalanche hazard and risk, respectively, in the Queyras massif (France), mostly due to agricultural abandonment. However, in these various studies, the diachronic methodology used to asses avalanche risk changes was mostly qualitative, which is not enough for supplying decision makers with diagnoses that are immediately usable.

On this basis, in this paper, we (1) include forest cover changes within a holistic quantitative avalanche risk assessment approach and (2) demonstrate on a typical case study of the French Alps how strongly the changes in forest cover that occurred over the 1825–2017 period may have affected avalanche hazard and subsequent risk for building-like elements at risk. By holistic, we mean an explicit modeling of forest cover changes effects on the full probability distribution of runout distances, impact pressures and on avalanche risk estimates for various types of buildings. Only the effect of forest cover change on avalanche release probability is neglected. To reach our goal, the Bayesian statistical–dynamical model of Eckert and others (2010b) is first expanded to account for multiple release areas within the same path and it is calibrated using on-site data. The local distribution of avalanche hazard is then tuned according to observed changes in the forest fraction (the latter being defined as the aerial percentage of the terrain covered by forests within the extension of the avalanche path). This is achieved by (1) evaluating forest cover changes within the case study extension from a combination of aerial photographs and ancient maps (Zgheib and others, 2020) and relating these changes to the friction coefficients μ and ξ of the Voellmy friction law. The resulting diachronic hazard distributions are eventually combined with fragility curves for different types of reinforced concrete (RC) buildings (Favier and others, 2014b), so as to produce first ever avalanche risk estimates that account for land cover changes in a comprehensive way. Note, however, that our aim is neither the in-depth modeling of the forest–avalanche interaction, nor of the damage due to avalanches for concrete, which would require more advanced numerical modeling techniques than those we use. Indeed, we rely on a rather simple frictional-like approach for avalanche–forest interactions, and on a large set of fragility curves evaluated under quasi-static assumptions. This allows our holistic risk assessment approach to be implemented at reasonable computational costs, so as to (1) grasp the on-site evolution of avalanche risk for buildings and (2) show how the protective effect of forests and building technology shall combine to limit the risk, which may ultimately allow the proper implementation of efficient gray-green protection measures against snow avalanches.

Case study

Our study site, the Ravin de Côte-Belle avalanche path (44°48′ N, 6°55′ E), is located in the Hautes Alpes department (Provence-Alpes-Côte d'Azur region) in the French Alps (Fig. 1a), on the northern slopes of the municipality of Abriés (Queyras massif) (Fig. 1b). Snow avalanches are mostly triggered from two distinct well localized release zones: zone I between 2500 and 2200 m a.s.l., and zone II between 1900 and 1800 m a.s.l. (Fig. 1c). Most of the avalanches stop in the runout zone ~1600 m a.s.l. The departmental road D411 linking le Roux village to Abriés cuts the runout area at the abscissa x = 1840 m (Figs 1c, 2a). Within the analysis, this position is taken as the location of a potential new building, so as to assess how risk to settlements has evolved as a function of forest cover extension within the path extension.

Fig. 1. Case study: Ravin de Côte-Belle avalanche path. (a) Location within the French Alps, (b) municipality of Abriés, Queyras massif, extension of avalanche paths from the EPA record, the Ravin de Côte-Belle avalanche path is highlighted in red, aerial photograph from 2017 (© IGN), (c) 2-D topography with historical data from the EPA record (the 17 fully documented events used for magnitude model calibration, Appendix E), and avalanche release zones I and II. Within the analysis, the road position at x = 1840 m is taken as the location of a hypothetical building, so as to assess how risk to settlements has evolved as a function of forest cover changes.

Fig. 2. Forest cover evolution at Ravin de Côte-Belle avalanche path. (a) Release zones I (2500–2200 m a.s.l.) and II (1900–1800 m a.s.l.) and diachronic map of forest cover extensions in 1825, 1948, 1980 and 2017. Elevation distribution of forest pixels in (b) 1825 (forest fraction f _k = 0.16), (c) 1948 (forest fraction f _k = 0.24), (d) 1980 (forest fraction f _k = 0.35) and (e) 2017 (forest fraction f _k = 0.46). Pixel size is 0.5 × 0.5 m².

Winter climate in the Queyras massif is cold, with relatively low precipitation in comparison with the rest of the French Alps (Durand and others, 2009a, 2009b). However, local avalanche activity is significant (Corona and others, 2013; de Bouchard d'Aubeterre and others, 2019), notably during ‘easterly returns’. The latter are atmospheric flows coming from the Mediterranean Sea that are responsible for heavy snowfall in the eastern part of the French Alps (Le Roux and others, 2021), leading to marked avalanche cycles (Eckert and others, 2010). As a consequence, the Ravin de Côte-Belle path is well documented in the French Avalanche Permanent Survey (referred to as EPA: Enquête Permanante des Avalanches (Bourova and others, 2016)), with 21 avalanches recorded between 1934 and 2018.

For the calibration of the avalanche statistical–dynamical model (Methodology section), all the 21 events from the EPA record were used to evaluate local avalanche occurrence frequency. By contrast, only 17 events (Appendix E) that occurred between 1960 and 2018 were sufficiently documented (release elevation, runout elevation, snow deposit), to be used for the calibration of the magnitude component of the model. Among these, 13 released from zone I and four released from zone II. Meteorological and snow conditions corresponding to the dates at which the events occurred were collected, notably snow depths at the elevation of the starting zones. These were taken from available reanalyses provided by the SAFRAN-Crocus model chain (Vernay and others, 2019).

Land cover changes from 1825 to 2017 within the Ravin de Côte-Belle were assessed using available historical maps and aerial photographs (Table 2, Appendix A) following the optimal combination approach proposed by Zgheib and others (2020). Historical maps of 1825 were georeferenced, and forest extensions were manually digitized. Regarding aerial photographs, pre-processed images were obtained for 2017 and 1948, whereas an orthorectification had to be applied to the 1980 image, and forest cover was digitized manually. For each date, the forest fraction, i.e. the aerial percentage of the terrain covered by forests within the extension of the avalanche path, could then be evaluated, as well as the more comprehensive elevation distribution of forest pixels. This shows that the forest fraction in the Ravin de Côte-Belle path increased all over the study period, and at accelerated pace over the last few decades, from 0.16 in 1825 to 0.24 in 1948, 0.35 in 1980 and to 0.46 in 2017 (Fig. 2a). In more detail, reforestation from 1825 to 2017 occurred mostly at high elevations of the path, between 2000 and 2200 m a.s.l. (Figs 2b–e). In 2017, the forest reaches the highest of the two avalanche release area (zone I), and completely covers the lowest one (zone II). Reforestation in the Queyras massif is a combined effect of the heavy depopulation and abandonment phase that began during the early 19th century and the development of the forest policy in 1860 that played an important role in forest management (Zgheib and others, 2022).

Methodology

A Bayesian statistical–dynamical model expanded to multiple release areas

The Bayesian statistical–dynamical model used in this study was developed by Eckert and others (2010b). According to Naaim and others (2004), it is based on the depth-averaged Saint Venant equations solved along a curvilinear profile z = f(x), where z is the elevation and x is the horizontal distance measured from the top of the avalanche path. Within the model, the following Saint Venant mass and momentum conservation laws represent a 1-D flow on the curvilinear profile solved using a finite volume scheme (Naaim, 1998). To facilitate the specification of the input conditions corresponding to each avalanche simulation and to reduce computation times, snow incorporation, deposition, entrainment and detrainment are ignored:

(1)$${\partial h\over \partial t} + {\partial ( hv) \over \partial x} = 0,\; $$

(2)$${\partial ( hv) \over \partial t} + {\partial ( hv^{2} + gr( {h^{2}}/{2}) ) \over \partial x} = h( gr\sin\phi-F) ,\; $$

where v is the flow velocity, h is the flow depth, ϕ is the local slope, t is the time, gr is the gravity acceleration and F is the total friction.

The total friction F considered is the classical Voellmy friction law (Voellmy, 1955):

(3)$$F = gr\mu \cos\phi + {gr\over \xi h}v^{2}.$$

Total friction is thus the sum of a dry-Coulomb coefficient μ and a turbulent term depending on the squared velocity and on the inverse of a coefficient ξ (Voellmy, 1955).

The full stochastic model (Fig. 3) representing the variability of avalanche events at the study site is noted p(y,  a|θ _M,  λ). Avalanche magnitude y is a random vector including all the correlated multivariate quantitative characteristics that vary from one event to another, namely runout distance, velocity and pressure profiles, or snow volume. Avalanche frequency a is a scalar discrete random variable corresponding to the number of avalanches recorded each winter. Avalanche magnitude and frequency are classically modeled as two independent random processes so that their joint distribution writes as a product and related parameters can be inferred separately (Eckert and others, 2010b):

(4)$$p( y,\; a\vert \theta_{M},\; \lambda) = p( y\vert \theta_{M}) p( a\vert \lambda) .$$

Fig. 3. Quantitative framework developed in this work to evaluate snow avalanche risk as a function of forest cover and building technology. The three parts of the analysis are as follows: (1) calibration of the statistical–dynamical model. The parameters of the model are: θ _M = (α ₁,  α ₂,  β ₁,  β ₂,  p,  b ₁,  b ₂,  σ _h,  c,  d,  e,  σ,  ξ) and θ _F = λ. (2) Simulation of the forest integration. Here three cases are considered: case I representing the dependency on μ, case II on ξ and case III on both μ and ξ. (3) Risk assessment and (4) risk calculation as a function of forest cover and building technology.

Avalanche frequency is modeled as a Poisson-distributed process, with a scalar parameter θ _F = λ representing the mean annual avalanche number i.e. a|λ ~ P(λ), necessary for the computation of return periods. The magnitude model specified below is more complex, with 13 parameters θ _M = (α ₁,  α ₂,  β ₁,  β ₂,  p,  b ₁,  b ₂,  σ _h,  c,  d,  e,  σ,  ξ). Also, the magnitude model evaluates, for each avalanche, the latent friction μ and the computed runout distance x _stop.

The studied path is characterized by the presence of two distinct starting zones: zone I between 2500 and 2200 m a.s.l. and zone II between 1900 and 1800 m a.s.l. (Figs 2a, b). To include this information into the analysis, rather than by the original single Beta distribution (Eckert and others, 2010b), we herein model x _start as a binomial mixture of two Beta distributions. This could be easily, in the future, generalized to even more complex cases using a multinomial mixture (e.g. Lavigne and others, 2012):

(5)$$ \eqalign{ \matrix { x_{{\rm start}_{i}} & \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!= P ( x_{{\rm start}_{i}} \in [ x_{\min_{1}},\; x_{\max_{1}}] ) x_{{\rm start}_{1}} + \cr & \,\,\,( 1-P ( x_{{\rm start}_{i}} \in [ x_{\min_{2}},\; x_{\max_{2}}] ) ) x_{{\rm start}_{2}},\; \cr & \!\!\!\!\!\!\! P ( x_{{\rm start}_{i}} \in [ x_{\min_{1}},\; x_{\max_{1}}] ) \sim B (p)}}$$

Since normalization is a requirement when dealing with the Beta distribution, normalized release abscissas are calculated and considered in the model as follows:

(6)$$\matrix{&x_{{\rm start}_{norm_{1}}} = {x_{{\rm start}_{1}}-x_{\min_{1}}\over x_{\max_{1}}- x_{\min_{1}}} \mathbbm{1}( x_{{\rm start}_{1}} \in [ x_{\min_{1}},\; x_{\max_{1}}] ) \sim Beta( \alpha_{1},\; \beta_{1}) ,\; \cr &x_{{\rm start}_{norm_{2}}} = {x_{{\rm start}_{2}}-x_{\min_{2}}\over x_{\max_{2}}- x_{\min_{2}}} \mathbbm{1}( x_{{\rm start}_{2}} \in [ x_{\min_{2}},\; x_{\max_{2}}] ) \sim Beta( \alpha_{2},\; \beta_{2}) ,\; }$$

where $x_{\min _{1}}$, $x_{\max _{1}}$, $x_{\min _{2}}$ and $x_{\max _{2}}$ are the minimal and maximal abscissas delimiting release zones I and II estimated for the case study using topographical thresholds (Figs 1, 2), and α ₁, α ₂, β ₁ and β ₂ the parameters of the corresponding Beta distributions, and $\mathbbm {1}(.)$ the indicator function.

The mean release depth h _start is assumed to be Gamma-distributed with parameters b ₁ and b ₂ reflecting its dependency on release abscissa and a constant dispersion around the mean σ _h (Eckert and others, 2010b). Here and in what follows the conditioning by $x_{\min _{1}}$, $x_{\min _{2}}$, $x_{\max _{1}}$ and $x_{\max _{2}}$ is dropped for simplicity:

(7)$$ \eqalign{ & h_{\rm start}\vert b_{1},\; b_{2},\; \sigma_{h},\; x_{\rm start} \sim \cr & Gamma\left({1\over \sigma_{h^{2}}}( b_{1} + b_{2}x_{{\rm start}_{\rm norm}}) ^{2},\; {1\over \sigma_{h^{2}}}( b_{1} + b_{2}x_{{\rm start}_{\rm norm}}) \right),\;}$$

where $x_{{\rm start}_{\rm norm}} = ( {x_{\rm start}-x_{\min _{1}}}) /( {x_{\max _{2}}-x_{\min _{1}}})$ is the normalized release abscissa evaluated all over the different release areas.

Given the normalized release abscissa $x_{{\rm start}_{\rm norm}}$ and the mean release depth h _start, the friction coefficient μ is modeled as a latent variable describing the random effects from one avalanche to another and it is assumed normally distributed (Eckert and others, 2010b):

(8)$$\mu\vert c,\; d,\; e,\; \sigma,\; x_{\rm start},\; h_{\rm start} \sim N( c + {\rm d} x_{{\rm start}_{\rm norm}} + eh_{\rm start},\; \sigma) .$$

The parameters c, d and e represent the dependency of μ on the release abscissa and mean release depth, with a constant dispersion σ around the mean. Parameters b ₂, d and e indirectly translate the impact of the altitude, and hence, of prevailing climate and snow conditions on the release depth and μ. Therefore the distribution of the μ and h _start varies according to the release altitude. Alternatively, the velocity-dependent friction coefficient ξ is modeled as a parameter in the strict statistical sense of the term. Both ξ and μ _i, i ∈ [1,  n] where n is the data sample size of fully documented events are estimated from the data.

Small Gaussian differences between the observed runout distances $x_{\rm stop_{\rm data}}$ and computed runout distances x _stop are postulated:

(9)$$x_{\rm stop_{\rm data}}\vert \sigma_{\rm num},\; x_{\rm stop} \sim N( x_{\rm stop},\; \sigma_{\rm num}) .$$

These differences can result from numerical errors due to the imperfection of the propagation model, and/or from observation errors. The SD of these numerical errors σ _num is set to 15 m to grant model identifiability.

Inference of the full model is a difficult problem which is solved by splitting it in simpler tasks. As stated before, the frequency and magnitude models are inferred separately. For the frequency model inference, all 21 avalanche events recorded since 1934 are used. For the magnitude model inference, only the n = 17 fully documented events are used. Avalanche events are assumed mutually independent. This implies that the joint likelihood of all events is the product of their individual marginal likelihood.

For the frequency model, Bayesian inference is analytical. With the chosen Gamma(a _λ,  b _λ) prior for the parameter λ, the posterior distribution of λ remains Gamma-distributed: Gamma(a _λ + T _obs,  b _λ + N _obs), where the pair (a _λ,  b _λ) represents the prior knowledge concerning the distribution of avalanche occurrences. N _obs is the total number of avalanches recorded on the study site (i.e. 21 events) during T _obs years of observation.

The parameters of the binomial Beta mixture model of x _start, α ₁, α ₂, β ₁ and β ₂, are estimated using the method of moments, a frequentist approach for parameter estimation (Appendix B).

The rest of the magnitude model is inferred using Bayesian methods, resulting in the joint posterior distribution of remaining parameters and latent variables:

(10)$$ \eqalign{ \matrix{\,p & \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! ( \theta,\; \mu,\; x_{\rm stop}\vert {data},\; \sigma_{\rm num}) \propto \pi( \theta) \cr & \times\, p( h_{\rm start},\; x_{\rm stop_{\rm data}}\vert \theta,\; \mu,\; x_{\rm start},\; x_{\rm stop},\; \sigma_{\rm num}) \cr & \times\, p( \mu,\; x_{\rm stop}\vert \theta,\; h_{\rm start},\; x_{\rm start},\; x_{\rm stop_{\rm data}},\; \sigma_{\rm num}) }}$$

where θ = (b ₁, b ₂, σ _h, c, d, e, σ, ξ), π(θ) is the joint prior distribution for the related parameters and data represents all observations. Numerical computation was achieved using the Metropolis Hasting algorithm within a Markov chain Monte Carlo scheme detailed in Eckert and others (2010b). Convergence was checked using several chains starting at different points of the space parameter. For each unknown various point estimates can be computed from the joint posterior distribution. We classically chose the posterior mean and summed-up the related uncertainty with the posterior SD and the 95% credible interval (Table 1).

Table 1. Parameter estimates of the statistical dynamical model

When Bayesian inference is used, for each parameter, the marginal prior distribution used is given, as well as summary statistics of the posterior distribution (the posterior mean is retained as point estimate). The mixture model for the release position (Eqn (5)) is calibrated separately using the method of moments (Appendix B).

Simulation of avalanche hazard conditional to local calibration

To evaluate avalanche hazard conditional to local calibration, 10 000 predictive simulations (Eckert and others, 2007, 2010b; Fischer and others, 2015, 2020) were performed with all parameters set to their Bayesian or frequentist point estimates (Table 1). In detail, for each simulation, x _start is evaluated according to the mixture model (Eqn (5)) that allows reconstructing a binomial mixture of two Beta distributions. Then, the normalized release abscissa injected in the simulation model is $x_{{\rm start}_{\rm norm}} = ( {x_{\rm start}-x_{\min {1}}}) /( {x_{\max _{2}}-x_{\min _{1}}})$, where x _start is the simulated release abscissa depending on its location on the path, $x_{\max _{2}}$ the maximal abscissa of zone II and $x_{\min _{1}}$ the minimal abscissa of zone I. A statistical–dynamical Monte Carlo approach is necessary to obtain the full joint distribution of the outputs of the numerical avalanche model knowing the distribution of the inputs. It involves integration over the distribution of the latent friction coefficient μ and writes as follows:

(11)$$ \eqalign{ & p( y \vert \,{ \hat {\! \theta} _{M}}) = \cr& \int p( x_{\rm start}\vert \hat {\alpha_{1}},\; \hat{\alpha_{2}},\; \hat{\beta_{1}},\; \hat{\beta_{2}},\; \hat{\,p}) \cr& p( h_{\rm start}\vert {\hat b_{1}},\; { \hat b_{2}},\; { \hat \sigma_{h}},\; x_{\rm start}) p( x_{\rm stop}\vert x_{\rm start},\; h_{\rm start},\; \mu,\; \hat{\xi}) {\rm d}\mu,\;}$$

where ‘$\hat {.}$’ classically denotes a statistical estimate. This simulation strategy leads, for instance, the joint distribution of the input variables $( x_{\rm start},\; \, h_{\rm start},\; \, \mu \vert { { \hat \theta} _{M}})$ and of the output variables $p( x_{\rm stop},\; \, v_{xt},\; \, h_{xt}\vert { \hat \theta _{M}})$ of the avalanche propagation model, where v _xt and h _xt represent, for each avalanche simulation, the velocity and flow depth computed for each abscissa and time step (Fig. 4).

Fig. 4. Avalanche hazard at Ravin de Côte-Belle: statistical–dynamical simulations according to local calibration and mean forest fraction f _k = 0.35 (as in 1980). Distribution of (a) release abscissa, (b) release depth, (c) runout abscissa, (d) friction coefficient μ, (e) maximal flow depth at abscissa position corresponding to a return period of T = 10 years and (f) maximal velocity at abscissa position corresponding to a return period of T = 10 years.

Subsequently, the return period $T_{x_{\rm stop}}$ associated with the runout distance x _stop, is estimated by combining $\hat {\lambda }$ and the cumulative distribution function (cdf) of runout distances $\hat {F}( x_{\rm stop}) = P( X_{\rm stop}\leq x_{\rm stop})$ as follows:

(12)$$T_{x_{\rm stop}} = {1\over \hat{\lambda}( 1-\hat {F}( x_{\rm stop}) ) }.$$

The inverse problem is then solved, to evaluate the runout distance quantile corresponding to the return period T as:

(13)$$x_{{\rm stop}_{T}} = \hat {F}^{-1}_{x_{\rm stop}}\left(1-{1\over \hat{\lambda}T}\right).$$

A Monte Carlo confidence interval is computed for the non-exceedance probability associated with a given runout distance. It allows checking if the sample size is large enough to give reliable estimates:

(14)$$CI_{\alpha} = \hat{F}( x_{\rm stop}) \pm q_{N_{\alpha_{c}}} \sqrt{{\hat {F}( x_{\rm stop}) ( 1-\hat {F}( x_{\rm stop}) ) \over n}},\; $$

where n is the sample size (in our case 10 000) and $q_{N_{\alpha _c}}$ is the quantile of the standard normal distribution corresponding to the desired confidence level α _c.

The runout return period obtained were used to extract the distribution of maximal velocities $v_{x}^{\max }$ (Fig. 4f) and maximal flow depths $h_{x}^{\max }$ (Fig. 4e) at abscissas corresponding to different return periods, notably 10 years.

Eventually, the distribution of impact pressures is classically calculated by transforming velocities into pressures as follows:

(15)$$P = C_{x}{1\over 2}\rho v^{2},\; $$

where C _x is the dimensionless coefficient of resistance i.e. drag coefficient, ρ is the snow density and v is the flow velocity. Studies like Sovilla and others (2008) and Naaim and others (2008) link velocity to pressure in a comprehensive but complex way that involves semi-empirical relationships between the drag coefficient C _x and the Froude number. For simplicity, we stick here to a value of C _x = 2, an approximation usable for wide, wall like structures, and ρ = 300 kg m⁻³ all along the analysis. Associated Monte Carlo confidence intervals are computed similar to Eqn (14).

Integration of forest cover changes within avalanche hazard assessment

The model calibrated using the above considers avalanche events as mutually independent, i.e. the result of the calibration is independent of the order of the events (except that each event is associated with the snow depth data in the release zone that prevailed at the date at which it occurred). Therefore, it is implicitly assumed that all events have occurred for the same path configuration and forest cover. This configuration corresponds to the average forest fraction over the period during which the events used for the calibration of the magnitude model occurred, namely 1950–2018. According (1) to the quasi-linear increase of forest fraction between 1948 and 2017 inferred from the aerial photographs, and (2) the slightly uneven temporal distribution of avalanche events within the EPA record over the 1950–2018 period, this mean forest fraction $\bar {f}$ arguably coincides with the forest fraction digitized from the 1980 aerial photograph, namely $\bar {f} = 0.35$.

Here, we seek at showing how the distribution of hazard and subsequent risk levels can be impacted by forest cover changes by introducing the forest fraction into the modeling in a robust way. Six values of f _k are considered, among which four represent the actual chronic of forest fractions digitized for the period 1825–2017 at the study site, and two extreme cases. In detail, f _k values considered are 0 (no forest, deforestation), 0.16 (forest fraction in 1825), 0.24 (forest fraction in 1948), $\bar {f} = f_{k} = 0.35$ (equal to the forest fraction in 1980), 0.46 (forest fraction in 2017) and 1 (the path is fully covered by forests).

As explained in the Introduction, in frictional approaches, forest–avalanche interaction is often modeled by decreasing the value of ξ significantly, and increasing the value of μ only slightly in the forested section of the path. However, it is known for long that the runout distance is mostly controlled by the static friction coefficient μ (e.g. Dent and Lang, 1980; Borstad and McClung, 2009), which has been recently confirmed by more systematic sensitivity analyses by Heredia and others (2021, 2022). Similarly, Teich and others (2012b) showed that focusing on the turbulent friction ξ only when investigating how forest cover impacts avalanche dynamics is not enough since the effects on runout distances of sole slight decrease of ξ remain too limited. Hence, we consider that changes in land cover can affect both ξ and μ and consider the three following cases, corresponding to the different possible relationships between forest fraction and friction coefficients: Case I: It is assumed that the dependability of the friction coefficient μ on land cover, particularly forests, can be modeled by adding a term $g( f_{k}-\bar {f})$ that increases/decreases the mean of the distribution of the static friction coefficient μ based on the value of the forest fraction f _k relative to $\bar {f}$. This uses the fact that μ is a latent variable in our statistical–dynamical model. A linear dependency between the mean of the distribution of μ and f _k is chosen, similar to the one between the mean of μ, release abscissa and mean release depth. To avoid the small number of unphysical negative values of μ that may result from this formulation with low values of f _k, we simply add a positivity constraint and sample the static friction coefficient μ from the truncated normal distribution (see Section ‘Dependence of the Voellmy friction coefficients on the forest fraction’ for discussion and Appendix D for the influence of the truncation on the case study):

(16)$$\mu \sim N( c + {\rm d} x_{{\rm start}_{\rm norm}} + eh + g( f_{k}-\bar{f}) ,\; \sigma) \mathbbm{1}( \mu > 0) ,\; $$

where the parameters c, d and e remain set to their posterior estimates (Table 1) resulting from the model calibration. Choosing a suitable value for g is a difficult task. To this aim, we use the fact that, when the effects of the variables x _start and h are neglected, the minimal value of the mean of μ, obtained for f _k = 0, is ${\bar \mu} _{\min } = c-g\bar {f}$. Physical knowledge and repeated calibrations in various avalanche terrains using deterministic or stochastic methods (Salm and others, 1990; Ancey and others, 2004; Naaim and others, 2013) indicate that classical very low values of μ are ~0.15. Within the equation $\bar {\rm \mu _{\min }} = c-g\bar {f}$, this corresponds to g = 0.6. Given that this remains, however, a strong assumption, we performed a sensitivity analysis with values of g spanning the full g = 0.4–0.8 range. Case II: It is assumed that the turbulent coefficient ξ decays as a power law with an increasing forest fraction:

(17)$$\xi = \hat{\xi}b^{( f_{k}-\bar{\,f}) }$$

where $\hat {\xi }$ is the posterior estimate resulting from the calibration step. This relationship uses the fact that ξ is a parameter in the strict statistical sense of the term. It is considered that, moving from deforested to forested terrain, ξ decreases by ~60% (~1000 to 400 m s⁻²) (Feistl, 2015). This corresponds to b = 0.5. However, just like in the previous case, the sensitivity of the results to this parameter is tested by spanning the range between b = 0.2 and b = 0.8. Eventually, the choice of a power law relationship stems from the rather limited protection offered by forests against natural hazard during the primary steps of forest colonization (Wohlgemuth and others, 2017) versus the maximum protection offered when reforestation of the avalanche path is complete. Case III: It is assumed that both μ and ξ depend on the forest fraction, combining the models corresponding to cases I and II. Results are analyzed and discussed for g = 0.6 and b = 0.5.

Note that, due to the specification of cases I–III, for $f_{k} = \bar {f} = 0.35$, friction coefficient models, and, hence, results of the previous section remain unchanged, so that no further simulation campaign is necessary. For all other forest fractions and in each case, 10 000 avalanches were simulated to reconstruct the joint distribution of model inputs–outputs according to Eqn (11). Subsequently, return periods, confidence intervals and impact pressures were computed based on Eqns (12–15).

Avalanche risk evaluation

Risk in the natural hazard field is defined as the expected damage resulting from the interaction between a damageable phenomenon and a vulnerable exposed entity. According to Eckert and others (2012), the specific (dimensionless) avalanche risk r _z for an element at risk z is:

(18)$$r_{z} = \lambda\int p( y) V_{z}( y) {\rm d}y,\; $$

where V _z(y) is the vulnerability of the element z to the avalanche magnitude y. Favier and others (2014b) further show that this generic expression can be rewritten at abscissa x _b as:

(19)$$r_{z}( x_{b}) = \lambda\int p( P\vert x_{b}\leq x_{\rm stop}) p( x_{b}\leq x_{\rm stop}) \times V_{z}( P) {\rm d}P.$$

Here, the avalanche magnitude distribution corresponds to the joint distribution of runout distances x _stop and pressures P i.e. p(P,  x _stop) = p(P|x _b ≤ x _stop)p(x _b ≤ x _stop) where p(P|x _b ≤ x _stop) is the pressure distribution at abscissa x _b knowing that x _b has been reached by an avalanche and p(x _b ≤ x _stop) is the probability for an element at x _b to be reached by an avalanche.

To assess the evolution of avalanche risk to settlements, we consider at x = 1840 m a typical mountainous dwelling house. Its overall vulnerability is determined by the failure probability of its most vulnerable part, i.e. the wall facing the avalanche (Favier and others, 2014a). Within a reliability framework, Favier and others (2014a) evaluated vulnerability curves for such typical buildings impacted by snow avalanches. This was done by modeling the response of an RC wall impacted by a uniform dense avalanche flow under the assumption of a quasi-static loading. These curves were obtained for various building types and for different limit state definitions: the elastic limit state (ELS), the ultimate limit state (ULS), the accidental limit state (ALS) and collapse (YLT: yield line theory). The ELS represents the upper limit of the elastic phase beyond which cracks begin to form and the concrete develops a non-linear behavior (Bertrand and others, 2010). At this point, the structure can still carry loads and the damage is considered low. Under continuously increasing pressure, the tensile crack will grow until the concrete or the steel reach respectively the ultimate compression strain and ultimate tensile strain (Favier and others, 2014a), thus announcing the reach of the ULSand the onset of steel yield (plastic (permanent) deformation inside the steel). The ALS is considered to ensure that the structure can withstand accidental events (statistically less likely to occur) e.g. explosions. Finally, yield lines or macro-crack form through the member leading to the collapse of the structure (YLT). For more details refer to Appendix C and Favier and others (2014a). In this study, we consider ten building types and the four limit state definitions. This results in a large set of 40 vulnerability relationships (Appendix C) providing the probability for the considered building type to surpass the considered limit state (i.e. the failure probability) when subjected to a given avalanche impact pressure. This large set allows studying in a robust way how risk to buildings varies with forest cover within the path, and, notably, how different combinations of building types and forest fractions can result in given risk levels.

Results

Impact of forest fraction on runout distances

The first part of the analysis examines the relationship between a varying forest fraction and runout distances. From the data in Fig. 5, it is apparent that the annual probability of an avalanche exceeding the house abscissa x _house = 1840 m significantly decreases with an increasing forest fraction. The largest variation is in case III (the forest fraction was introduced as a parameter impacting both the static friction μ and the turbulent friction ξ), followed very closely by case I (forest fraction acts on μ), whereas modification is much lower in case II (the forest fraction is linked to the velocity-dependent friction ξ). Considering the two extreme situations f _k = 0 (deforestation) and f _k = 1 (complete reforestation), this corresponds to a decrease from a 17.6% annual exceedance probability (f _k = 0) to 0.25% (f _k = 1) (Fig. 5c) in case III, from 17% (f _k = 0) to 0.36% (f _k = 1) (Fig. 5a) in case I, and from 10% (f _k = 0) to 7% (f _k = 1) only in case II (Fig. 5b). From a temporal point of view, focusing on changes that actually occurred within the path, for cases I and III, the probability of an avalanche exceeding the house abscissa x _house = 1840 m decreased from 14% in 1825 (f _k = 0.16) to 6% in 2017 (f _k = 0.46) (Figs 5a, c), whereas in 1980 (f _k = 0.35, mean forest fraction for the study period) it was ~9%. Eventually, note that, along the path, the runout exceedance probability varies more with forest fraction for abscissas $\eqslantless \!1840$ m than for abscissas >1840 m (Fig. 5).

Fig. 5. Annual exceedance probability of runout distances for: (a) case I (f _k acts only on the static friction coefficient μ), (b) case II (f _k acts only on the turbulent friction coefficient ξ) and (c) case III (f _k acts on both the static friction coefficient μ and the turbulent friction coefficient ξ). Results are shown for the six forest fractions: f _k = 0 (deforestation), 0.16 (as in 1825), 0.24 (as in 1948), 0.35 (mean forest fraction, as in 1980), 0.46 (as in 2017), 1 (complete reforestation). Due to the forest fraction integration model specification (Eqn (16–17)), by definition, for $f_{k} = \bar {f} = 0.35$ (as in 1980), results are the same in the three cases and correspond to those obtained conditional to local calibration (Fig. 4c).

With the mean forest fraction corresponding to 1980 ($f_{k} = \bar {f} = 0.35$), the return period is 11 years for events reaching the x _house = 1840 m abscissa. Logically, given the changes observed for exceedance probabilities, cases III (Fig. 6f) and I (Fig. 6d), have the largest impact on the return period at x _house = 1840, whereas case II (Fig. 6e) leads to the lowest variations. For the two extreme situations f _k = 0 (deforestation) and f _k = 1 (complete reforestation), T rises from 5.7 years (case III, f _k = 0) to T = 393 years (case III, f _k = 1), from 5.9 years (case I, f _k = 0) to T = 273 years (case I, f _k = 1) and from T = 10 years (case II, f _k = 0) to T = 14 years (case II, f _k = 1) only (Fig. 6e). From a temporal point of view, for cases I and III, events reaching X _house had a return period of ~7 years in 1825 which rose to 15.5 years in 2017 (Figs 6a, c, d, f). In general, largest variations in the runout distance–return period relationship with the forest fraction are for runout abscissas ≥1600 m, namely the complete lower part of the path where the return period increases strongly with abscissa (Figs 6a–c). All Monte Carlo confidence intervals are very small (Figs 6d–f), showing that the number of simulations performed is large enough to highlight significant changes of runout distance distributions with forest fractions.

Fig. 6. One-to-one relationship between runout distances and return periods for (a) case I (f _k acts only on the static friction coefficient μ), (b) case II (f _k acts only on the turbulent friction coefficient ξ) and (c) case III (f _k acts on the static and turbulent friction coefficient μ and ξ respectively). (d, e and f) Close up on the one-to-one relationship between runout distances >1600 m and return periods, and the associated 95% confidence intervals for the return period (Eqn (14)). Results are shown for six forest fractions: f _k = 0 (deforestation), 0.16 (as in 1825), 0.24 (as in 1948), 0.35 (mean forest fraction, as in 1980), 0.46 (as in 2017), 1 (reforestation). Due to the forest fraction integration model specification (Eqn (16–17)), by definition, for f _k= $\bar {f} = 0.35$ (as in 1980), results are the same in the three cases and correspond to those obtained conditional to local calibration (Fig. 4c).

Impact of forest fraction on pressures

Figure 7 shows the annual probability of occurrence of pressure values at three distinct positions of the avalanche path, for all cases and forest fractions studied. In general, results show that, for all three positions, the annual probability of having high pressures decreases with higher forest fractions, with cases III and I leading to the largest changes with forest fraction. In detail, the annual probability for an impact pressure ⩾30 kPa at x _house varies between 5 and 20% for all cases, except for complete reforestation (f _k = 1) in cases I and III where it drops to nearly 0 (Figs 7a, d). Considering the two extreme situations f _k = 0 (deforestation) and f _k = 1 (complete reforestation), the annual exceedance probability p(Pressure_house⩾30 kPa) at x _house = 1840 m varies from 16% (f _k = 0) to 0.2% (f _k = 1) for case I, and from 18% (f _k = 0) to 0.1% (f _k = 1) for case III (Figs 7a, d). Here also, case II has no large impact on the annual exceedance probability compared to the other two cases since p(Pressure_house ⩾30 kPa) varies from 9% (f _k = 0) to 5% (f _k = 1) only. From the temporal point of view, for cases I and III, the annual probability of impact pressures ⩾30 kPa at x _house decreased from ~13% in 1825 (f _k = 0.16) to 5% in 2017 (f _k = 0.46) (Figs 7a, d), whereas in 1980 ($f_{k} = \bar {f} = 0.35$) it was 7.6%. Again, 95% confidence intervals demonstrate the significance of changes according to the simulation sample size (Fig. 7d). Regarding the location within the path, high impact pressures are more likely to occur in the propagation zone i.e. at 1000 m (Fig. 7b), followed by the runout zone, here represented by the house abscissa x _house = 1840 m (Fig. 7c) and lastly in the starting zone at x = 400 m (Fig. 7a). For example, in 1980 (f _k = 0.35, i.e. the mean forest fraction), the annual probability of pressures exceeding 250 kPa was 4% at x = 1000, 2% at x = 1840 m and 1% only at x = 400 m (Fig. 7).

Fig. 7. Annual exceedance probability of impact pressure modeled for the three cases I (f _k acts only on the static friction coefficient μ), II (f _k acts only on the turbulent friction coefficient ξ) and III (f _k acts on both the static and turbulent friction. Results are provided for three abscissa positions within the path: (a) at x _house = 1840 m located in the runout zone, (b) at x = 1000 m within the propagation zone and (c) at x = 400 m located in zone I (release area) (Fig. 1d). (d) Close up on the annual exceedance probability of impact pressure at X _house = 1840 m, and the associated 95% confidence intervals. Results are shown for the six forest fractions: f _k = 0 (deforestation), 0.16 (as in 1825), 0.24 (as in 1948), 0.35 (mean forest fraction, as in 1980), 0.46 (as in 2017), 1 (complete reforestation). Due to the forest fraction integration model specification (Eqns (16–17)), by definition, for $f_{k} = \bar {f} = 0.35$ (as in 1980), results are the same in the three cases and correspond to those obtained conditional to local calibration.

Impact of forest fraction on risk for buildings

Figure 8 depicts the risk for the house located at x = 1840 m (T = 11 years in 1980, $f_{k} = \bar {f} = 0.35$), for all building types and the four limit states. At this location, in 1980, avalanche risk (i.e. the annual probability to reach a given critical state) was between 0.02 and 0.09 depending on the considered limit state and building type. Unsurprisingly given the definition of the limit states, for all building types and for a given forest fraction, the ELS (Fig. 8a) based risk > ULS risk (Fig. 8b) > ALS risk (Fig. 8c) > YLT-based risk (Fig. 8d). The strongest building type is the building VIII with four clamped edges and the weakest is configuration IV i.e. one free and three simply supported edges (Fig. 8). Overall, the weakest buildings include at least one free edge, whereas the strongest have at least two clamped and no free edges.

Fig. 8. Risk i.e. annual probability for a building, located at x = 1840 m within the Ravin de Côte-Belle avalanche path to reach the limit state considered (Fig. 11, Appendix C). The analysis considers ten building types (Table 3, Appendix C), four limit states and the three forest integration cases I (f _k acts only on the static friction coefficient μ), II (f _k acts only on the turbulent friction coefficient ξ) and III (f _k acts on both the static and turbulent friction coefficient μ and ξ respectively). All forest fractions are considered. Due to the forest fraction integration model specification (Eqns (16)–(17)), by definition, for $f_{k} = \bar {f} = 0.35$ (as in 1980), results are the same in the three cases and correspond to those obtained conditional to local calibration.

From the results, it is immediately clear that, for a given building type and limit state, the risk decreases with an increasing forest fraction. Yet, whereas the ELS-based risk to all building types varies almost similarly with forest fraction in all cases (Fig. 8a), for the other failure states (Figs 8b–d), risk varies according to the building type, the considered case and forest fraction in a more complex way. Considering the two extreme situations f _k = 0 (deforestation) and f _k = 1 (complete reforestation), cases III and I show the largest risk variations according to forest fraction changes for all limit states and building types. For example, the annual probability for a building located at x = 1840 m to reach the ELS drops for cases I and III from >0.16 (f _k = 0) to 0.00 (f _k = 1), but by no more than 0.03 for case II (Fig. 8a). From a temporal point of view, between 1825 (f _k = 0.16) and 2017 (f _k = 0.46), the risk roughly averaged over all building types and limit states decreased by almost 60% for cases I and III and by ~20% for case II. Given the very slight Monte Carlo uncertainty regarding pressure distributions, even small differences in risk according to building type and/or failure state are significant (Fig. 7).

Sensitivity to the forest fraction integration model

Sensitivity of our results to the forest fraction integration model was investigated by varying the parameters g and b controlling the degree of impact of the forest fraction on the annual exceedance probability of runout distances (Fig. 9) and, hence, on avalanche risk for a building located at x = 1840 m (T = 11 years in 1980 for $f_{k} = \bar {f} = 0.35$) (Fig. 10). The main result is that, over the large considered variation ranges, forest fraction changes affect the annual probability of exceedance of runout distances and the risk for buildings in a rather similar way regardless of the values chosen for g (case I) and b (case II). Obviously, in detail, results are slightly shifted when g and b are modified. For example, the higher g, the faster μ values increase with forest fraction, and the faster exceedance probabilities and ultimately risks for building decrease with forest fraction. More importantly, both in cases I and II, exceedance probability bands for f _k = 0 and f _k = 1 (in green and blue respectively) never intersect each other (Fig. 9), which suggest that the forest fraction value has a much more decisive effect on our results than the choice of given g and b value (as long as a reasonable range is considered).

Fig. 9. Sensitivity of the annual exceedance probability of runout distances to the forest fraction integration model in (a) case I (f _k acts only on the static friction coefficient μ) and (b) case II (f _k acts only on the turbulent friction coefficient ξ). g and b are the parameters representing the dependency of μ and ξ on the forest fraction f _k, respectively (Eqns (16)–(17)). represents the annual exceedance probability band delimited by g = 0.8 and g = 0.4. represents the annual exceedance probability band delimited by b = 0.8 and b = 0.2. Only the two extreme forest fractions are considered, i.e. deforestation (f _k = 0) and complete reforestation (f _k = 1).

Fig. 10. Risk sensitivity to the forest fraction integration model. Risk is the annual probability for a building located at x = 1840 m within the Ravin de Côte-Belle avalanche path to reach the limit state considered. The analysis considers ten building types, four limit states and the two forest integration cases I (f _k acts only on the static friction coefficient μ) and II (f _k acts only on the turbulent friction coefficient ξ). g and b are the parameters representing the dependency of μ and ξ coefficients on the forest fraction f _k, respectively. Only the two extreme forest fractions are considered, i.e. deforestation (f _k=0) and complete reforestation (f _k = 1).

Discussion

Avalanche risk changes with changes in forest cover

This study aimed at assessing the impact of the changes in the forest fraction on avalanche hazard and risk to a hypothetical building located at x = 1840 m in the Ravin de Côte-Belle path (Queyras massif). Unsurprisingly, our findings showed, for an increasing forest fraction, a decrease in p(X _stop ≥ x _house = 1840 m) (increase in return period) (Fig. 5), in the probability of high impact pressures at the house abscissa (Fig. 7a) and eventually in avalanche risk (Fig. 8). All these results originated from the combined effect of the forest fraction on runout distances and velocities, i.e. increasing f _k reduces both reduces the annual runout probability and the velocity conditional to reach, which results in decreased annual pressure probabilities, and, ultimately, lower risk levels. These results are in line with previous studies notably Teich and Bebi (2009) whom also highlighted a decrease in avalanche risk linked to the spatial variability and extent of the forest cover. Hence, using a more systematic risk-based methodology, we confirm and quantify the importance of the protective effect of forests against snow avalanches.

From a temporal perspective, we showed that the probability for a snow avalanche reaching the building decreased by almost half between 1825 and 2017 as a result of path's reforestation (forest fraction tripled between 1825 and 2017). Consequently, avalanche risk for buildings also decreased by more than a half (Fig. 8). This is consistent with the more qualitative findings of Mainieri and others (2020) and Zgheib and others (2022) regarding the evolution of hazard and risk in the Queyras massif. These authors suggested that both avalanche hazard and risk decreased sharply after 1950, linked to reforestation of the avalanche paths, the latter being a result of the socio-economic transitions (i.e. abandonment, tourism and changes in forest policies) that took place in the Queyras massif (Granet-Abisset, 1991) and in the European Alps in general (MacDonald and others, 2000) since around mid-19th century. It is therefore likely that, in the Queyras massif, avalanche risk trends similar to the one we highlight exist in many paths that went through the same intense reforestation process. Extrapolation beyond the Queyras massif is more difficult seeing that, despite common global socio-economic and environmental transitions in the mountains of the world, local disparities are always present and they highly impact the trajectory evolution of hazard and risk (Zgheib and others, 2022). Hence, similar studies should now be conducted in other massifs to indicate to which extent our findings apply to wide mountain areas or not. Given (1) the very strong decrease in risk our integrated quantitative methodology was able to identify, (2) the current lack of consideration of temporal changes in risk levels with land use and climate (Eckert and others, 2018), such studies would be of uttermost importance for improving both safety and sustainability of mountain communities.

A potential for combined nature-based and structural protection measures

Avalanche risk was assessed for ten different building types and four limit states for RC. It was quite a surprise to find that avalanche risk to all building types is relatively similar if the ELS is considered as the failure mode. This result suggests that preventing the elastic failure of the wall is independent of the building configuration. However, this is only valid for avalanches with an impact pressure ≤36 kPa, i.e. the maximum impact pressure needed for the elastic failure of the wall (all configurations considered, Appendix C, Fig. 12). Alternatively, for large avalanches exerting higher impact pressures, the building type is a decisive factor and stronger construction types logically provide additional safety.

Also, our results show that different forest fraction/building configurations lead to the same risk level downslope (Fig. 8). Hence, the same risk reduction can be ensured through an array of possible combinations of building reinforcement and forest management. This shows the potential for combining nature-based and traditional engineering solutions in snow avalanche risk mitigation. However, it must be kept in mind that forest disturbances (e.g. windthrow, bark beetle outbreaks, forest fires, etc.) challenge ecosystem services of mountain forests such as avalanche protection (Teich and others, 2019; Stritih and others, 2021; Caduff and others, 2022). In turn, forest disturbances are sensitive to climate changes (Maroschek and others, 2015; Seidl and others, 2017) that alters their frequency, intensity and duration (Dale and others, 2001). Hence, complex interactions and factors should be kept in mind to fairly assess the protection capacity of forests and taking it into account in the selection of the best risk management strategy.

Even more broadly, avalanche protection largely relies on structures capable of withholding the flow or inducing runout shortening, and different approaches exist to optimize their efficiency (Faug and others, 2008; Eckert and others, 2012; Favier and others, 2022). However, considering the large temporal evolution of forest cover in several areas, our results suggest that they could be, in the future, more systematically combined to greener solutions. Indeed, Ecosystem-based solutions for Disaster Risk Reduction (Eco-DRR) such as protective forests are already known to be very efficient, when combined with structural avalanche protection techniques (i.e. snow sheds, snow fences, snow bridges, etc.) to reduce collective risk in densely populated areas exposed to snow avalanches (Teich and Bebi, 2009). Of course, this may not always be possible, as, e.g. above the treeline, supporting structures made of steel remain the only option to prevent avalanche release. Also, intrinsic limitations of Eco-DRR solution such as their sensitivity to disturbances should be kept in mind. Yet, our results confirm that they may be, in some cases, an effective, economically viable (Poratelli and others, 2020) option to be considered for avalanche risk management. For instance, below the treeline, permanent avalanche supporting fences combined with afforestation can probably be a sensible solution in several avalanche paths, as an alternative to larger ‘gray’ protection structures.

Dependence of the Voellmy friction coefficients on the forest fraction

To reach our conclusions, a key step of the approach was the representation of the dependency of the Voellmy friction coefficients on the forest fraction. Notably, changing the static friction coefficient μ based on the forest fraction f _k, i.e. case I, had the largest impact on all the studied variables (runout return period, impact pressure and risk) compared to varying the friction coefficient ξ (case II). This relates to the high sensitivity to changes in μ of runout distances in accordance with previous studies (Barbolini and others, 2000; Heredia and others, 2021, 2022). To which extent our results are direct consequences of partially subjective model characteristics is thus debatable.

First argument in favor of our choices for the forest/friction dependency is that they are consistent with our statistical–dynamical model structure, i.e. addition of a fixed effect in the linear model for the variable μ and modification of the parameter ξ with an exponential dependency on f _k, as well as with physical and empirical knowledge regarding current values of μ and ξ and their variations with terrain properties. More pragmatically, the sensitivity analysis highlighted that most of the changes in the annual probability of exceedance of x _house and in related avalanche risk to buildings are explained by variations in the forest fraction and not by parametric choices. Thus, our results, notably the sharp risk decrease from 1825 to 2017 appear as robust.

Similarly, for μ, our modeling choice made that we had to include a positivity constraint within the simulation set-up (Eqn (16)). Appendix D shows that the fraction of μ values rejected by the truncation remains, however, low, except for the two lowest forest fractions considered. In addition, it is clear that, even with the truncation, the distribution of generated μ values well decreases with forest fraction (notably the mean of the truncated distribution, Appendix D). This suggests that our simple numerical trick is sufficient to support our conclusions.

However, in the future, more refined modeling of the functional relation between μ and ξ and forest cover could be envisaged. Notably, b and g could be considered within the calibration of the statistical–dynamical model. Also, additional information regarding forest structure, if available, could be included by proposing explicit functional relationships with, e.g. forest fraction, stem number and density instead of having within the model the forest fraction only. At the cost of additional inference difficulties (e.g. a potentially high number of parameters to calibrate with a limited amount of data and a numerically intensive procedure), this would allow a more in-depth analysis of the relationship between forest cover characteristics and friction parameters, and, ultimately, avalanche risk.

Other pros and cons of the modeling strategy

In addition to potentially debatable statistical modeling assumptions, this work also suffers from a number of other limitations. Notably, the proposed approach does not consider the altitudinal distribution of the forest (Figs 2b–e). Yet, depending on its position, a protective forest can either stabilize the snow and prevent avalanching (Salm, 1978; Viglietti and others, 2010) (location in the starting zone of the avalanche), or decelerate a flowing avalanche (Anderson and McClung, 2012) (location in runout or propagation zone). For instance, our results do not reflect the possibility that the avalanche release zone II has potentially become inactive in 2017. This is highly likely considering that, in 2017, the density of forest pixels between 1800 and 1900 m a.s.l. and above is very high compared with considered previous dates (Fig. 2e). More broadly, forest protection efficiency also depends on forest structural characteristics, the magnitude of the avalanche (small, medium, large) and more specifically its velocity, density and type (powder or dense flow avalanche). For example, large avalanches initiating way above the timberline can destroy forests with negligible deceleration (Brang and others, 2006) but noticeable lateral spread (Christen and others, 2010). This statement remains under discussion in light of studies showing that the runout distance of large avalanches can be decreased despite the destruction of trees (Takeuchi and others, 2018). Our approach does not grasp these complex modulations of forest protective effect as a function of avalanche and forest characteristics, notably how protection efficiency may have changed over the study period as function of changes in tree species and density. Eventually, a basic assumption underlying our modeling is that an increase in forest cover within the path always reduces (or at least does not increase) the risk. This sounds generally true, except when more forest leads to a higher incorporation of trunks within avalanche flows that can generate high impact pressures, and, hence, potentially increase the risk to buildings and people inside. These different shortcuts could potentially be relaxed in future developments, notably by including more advanced numerical modeling techniques within the workflow to represent the complex interactions between avalanche flows, forest stands and buildings (Védrine and others, 2021). However, a computational effort much higher than the one required to implement our approach would then be necessary.

By contrast, the novelty of the proposed holistic quantitative approach for the snow avalanche field is clear. It lays in (1) consideration of the full variability of avalanche hazard, conditional on the temporal evolution of forest fraction in the avalanche path, through comprehensive hazard statistical–dynamical modeling, (2) taking into account different building types and failure states in the evaluation of risk changes. This step forward with regards to the state-of-the-art (1) successfully highlights the local decrease in avalanche risk for buildings according to observed changes in forest cover, (2) features the important protective capacity of forests against snow avalanches and (3) points toward the potential of different combinations between forest cover and building types to manage and reduce the risk to the desired level (taking into consideration that the protective effect of forests may be sharply reduced by different disturbances, see above). Note eventually that all computations remain feasible on a standard computer.

Conclusion and further outlook

This study proposed an innovative holistic risk analysis methodology to evaluate diachronic risk estimates that account for land cover changes within the avalanche path in a comprehensive way. First, a Bayesian statistical–dynamical model was expanded to account for multiple release areas, and then locally calibrated for the Ravin de Côte-Belle avalanche path in Abriés (Queyras massif, French Alps). Changes in the distribution of avalanche hazard were then evaluated according to the observed changes in the aerial percentage of the terrain covered by forests. Results were eventually combined with RC fragility curves for different types of buildings and failure states. This allowed quantifying to which extent avalanche hazard and risk to buildings decreases with increasing forest cover. Notably, between 1825 and 2017, avalanche risk in the Ravin de Côte-Belle decreased by 20–60% depending on how forest fraction is accounted for in avalanche statistical–dynamical modeling. In addition, we showed that the decrease in risk depends not only on the forest fraction but also on the construction technology used for the building in the avalanche path. This should, in the future, allow objective cost/benefit analyses of the economic viability of protective forests versus structural mitigation measures (Moos and others, 2018). It could also facilitate the setup of successful disaster risk reduction strategies that consider both the building technology used and the type of avalanche protection (i.e. forest, structural measures or combination of the latter) to ensure the reduction of risk to acceptable level, in other words to intelligently combine Eco-DRR and technical protection measures. Such an integrated approach could be an important step forward toward an improved understanding of ecosystem services in mountainous terrain fully accounting for process uncertainties (Stritih and others, 2019) that may help answering stakeholders needs.

Eventually, several studies posit drastic ongoing changes in risks due to snow avalanches and other natural hazards in the alps and wider mountain areas (Eckert and others, 2013; Ballesteros-Cánovas and others, 2018; Hock and others, 2019; Giacona and others, 2021; Schlögl and others, 2021). Since evolution of avalanche risk depends on the interactions between hazard, exposure and vulnerability driven by socio-economic, land cover and climatic drivers (Zgheib and others, 2020, 2022), a dynamic quantitative risk assessment should become the basis of future disaster risk reduction, mitigation plans and policies. Although being arguably a first step, our approach is still far from this ambitious objective. We indeed considered only the impact of the temporal evolution of forest on risk explicitly, and through partially ad-hoc simulations. The latter limitation could be addressed, as suggested above, by calibrating the link between the forest fraction and the Voellmy friction coefficients. However, to take into account the combined effect of a changing climate and land cover on avalanche hazard, the current model should now be made fully non-stationary to derive frequency-magnitude relationships changing over time in a comprehensive and realistic way. Combined with elements at risk evolving in terms of vulnerability and/or exposition as a function of socio-economic transitions, this would provide crucial support to the great challenge of climate change adaptation.