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, Reference Delage2003; Moos and others, Reference Moos2018). 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, Reference Kumar2020; Poratelli and others, Reference Poratelli2020a). 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, Reference Brang, Schönenberger, Ott and Gardner2001, Reference Brang2006). The protective effects of forest stands against, e.g. rockfalls (Dupire and others, Reference Dupire2016) and snow avalanches (Bebi and others, Reference Bebi, Kulakowski and Rixen2009), is of a great importance in alpine regions, preventing human deaths and destruction of buildings and infrastructures (Getzner and others, Reference Getzner, Gutheil-Knopp-Kirchwald, Kreimer, Kirchmeir and Huber2017). 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, Reference De Quervain1978; Salm, Reference Salm1978; Gubler and Rychetnik, Reference Gubler and Rychetnik1991; Viglietti and others, Reference Viglietti, Letey, Motta, Maggioni and Freppaz2010; Teich and others, Reference Teich, Bartelt, Grêt-Regamey and Bebi2012a). They also have the capacity to decelerate flowing avalanches (Anderson and McClung, Reference Anderson and McClung2012; Takeuchi and others, Reference Takeuchi, Nishimura and Patra2018). Their protective potential is directly linked to forest structural parameters, e.g. the stem density for small-medium avalanches (Teich and others, Reference Teich, Bartelt, Grêt-Regamey and Bebi2012a, Reference Teich2014), and the distance traveled before penetrating into the forest for large avalanches triggered above the timberline (McClung, Reference McClung2003; Takeuchi and others, Reference Takeuchi, Torita, Nishimura and Hirashima2011; Teich and others, Reference Teich, Bartelt, Grêt-Regamey and Bebi2012a).
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, Reference Salm1978; Gubler and Rychetnik, Reference Gubler and Rychetnik1991; Bartelt and Stöckli, Reference Bartelt and Stöckli2001; Takeuchi and others, Reference Takeuchi, Torita, Nishimura and Hirashima2011). 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, Reference Voellmy1955). The static friction coefficient μ generally varies within the 0.1–0.7 range, and is often thought to summarize snow properties (Salm and others, Reference Salm, Burkard and Gubler1990). The turbulent friction coefficient ξ generally varies within the 1000–10 000 m s−2 range for forest-free terrain, and aims at representing the roughness of the path potentially related to land cover properties (Salm and others, Reference Salm, Burkard and Gubler1990; Ancey and others, Reference Ancey, Meunier and Richard2003). This explains why, in frictional approaches, more often than not, ξ is lowered to values in the order of 400 m s−2 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, Reference Bartelt and Stöckli2001). 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, Reference Gruber and Bartelt2007; Christen and others, Reference Christen, Bartelt and Kowalski2010). More recently, as an alternative to frictional approaches, Feistl and others (Reference Feistl2014) 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, Reference Védrine, Li and Gaume2021).
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, Reference Harbitz, Issler and Keylock1998). 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, Reference Bründl, Romang, Bischof and Rheinberger2009). First implementations relied on simple scenarios representing extreme avalanches (Fuchs and others, Reference Fuchs, Bründl and Stötter2004), 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, Reference Keylock, McClung and Magnússon1999; Barbolini and others, Reference Barbolini, Cappabianca and Savi2004; Cappabianca and others, Reference Cappabianca, Barbolini and Natale2008; Eckert and others, Reference Eckert, Parent, Faug and Naaim2008, Reference Eckert, Parent, Faug and Naaim2009; Favier and others, Reference Favier, Eckert, Bertrand and Naaim2014b, Reference Favier, Eckert, Faug, Bertrand and Naaim2016). This was clearly a step toward more integrated avalanche risk mitigation solutions potentially accounting, e.g.for acceptability thresholds (Arnalds and others, Reference Arnalds, Jónasson and Sigurosson2004) 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, Reference Grêt-Regamey and Straub2006) or scenarios (Teich and Bebi, Reference Teich and Bebi2009) 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, Reference Moos2018; Farvacque and others, Reference Farvacque2019) and it is a rather common strategy (although still difficult to properly achieve) for flood risk (Rogger and others, Reference Rogger2017; Bathurst and others, Reference Bathurst, Fahey, Iroumé and Jones2020).
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, Reference Mather, Fairbairn and Needle1999; Bebi and others, Reference Bebi2017). This time frame roughly corresponds to the classical 100–300 year reference periods that define legal thresholds in land-use planning (Eckert and others, Reference Eckert2018). 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 (Reference Giacona2018) 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 (Reference García-Hernández2017) 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 (Reference Mainieri2020) and Zgheib and others (Reference Zgheib2022) 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 (Reference Eckert, Naaim and Parent2010b) 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, Reference Zgheib, Giacona, Granet-Abisset, Morin and Eckert2020) 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, Reference Favier, Eckert, Bertrand and Naaim2014b), 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 m2.
Winter climate in the Queyras massif is cold, with relatively low precipitation in comparison with the rest of the French Alps (Durand and others, Reference Durand2009a, Reference Durand2009b). However, local avalanche activity is significant (Corona and others, Reference Corona2013; de Bouchard d'Aubeterre and others, Reference de Bouchard d'Aubeterre2019), 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, Reference Le Roux, Evin, Eckert, Blanchet and Morin2021), leading to marked avalanche cycles (Eckert and others, Reference Eckert2010). 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, Reference Bourova2016)), 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, Reference Vernay2019).
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 (Reference Zgheib, Giacona, Granet-Abisset, Morin and Eckert2020). 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, Reference Zgheib2022).
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 (Reference Eckert, Naaim and Parent2010b). According to Naaim and others (Reference Naaim, Naaim-Bouvet, Faug and Bouchet2004), 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, Reference Naaim1998). 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:
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, Reference Voellmy1955):
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, Reference Voellmy1955).
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, Reference Eckert, Naaim and Parent2010b):
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 = (α 1, α 2, β 1, β 2, p, b 1, b 2, σ 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 = (α 1, α 2, β 1, β 2, p, b 1, b 2, σ 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, Reference Eckert, Naaim and Parent2010b), 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, Reference Lavigne, Bel, Parent and Eckert2012):
Since normalization is a requirement when dealing with the Beta distribution, normalized release abscissas are calculated and considered in the model as follows:
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 α 1, α 2, β 1 and β 2 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 1 and b 2 reflecting its dependency on release abscissa and a constant dispersion around the mean σ h (Eckert and others, Reference Eckert, Naaim and Parent2010b). 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:
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, Reference Eckert, Naaim and Parent2010b):
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 2, 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:
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, α 1, α 2, β 1 and β 2, 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:
where θ = (b 1, b 2, σ 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 (Reference Eckert, Naaim and Parent2010b). 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, Reference Eckert, Parent and Richard2007, Reference Eckert, Naaim and Parent2010b; Fischer and others, Reference Fischer, Kofler, Fellin, Granig and Kleemayr2015, Reference Fischer2020) 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:
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:
The inverse problem is then solved, to evaluate the runout distance quantile corresponding to the return period T as:
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:
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:
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 (Reference Sovilla, Schaer, Kern and Bartelt2008) and Naaim and others (Reference Naaim2008) 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−3 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, Reference Dent and Lang1980; Borstad and McClung, Reference Borstad and McClung2009), which has been recently confirmed by more systematic sensitivity analyses by Heredia and others (Reference Heredia, Prieur and Eckert2021, Reference Heredia, Prieur and Eckert2022). Similarly, Teich and others (Reference Teich2012b) 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):
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, Reference Salm, Burkard and Gubler1990; Ancey and others, Reference Ancey, Gervasoni and Meunier2004; Naaim and others, Reference Naaim, Durand, Eckert and Chambon2013) 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:
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−2) (Feistl, Reference Feistl2015). 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, Reference Wohlgemuth, Schwitter, Bebi, Sutter and Brang2017) 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 (Reference Eckert2012), the specific (dimensionless) avalanche risk r z for an element at risk z is:
where V z(y) is the vulnerability of the element z to the avalanche magnitude y. Favier and others (Reference Favier, Eckert, Bertrand and Naaim2014b) further show that this generic expression can be rewritten at abscissa x b as:
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, Reference Favier, Bertrand, Eckert and Naaim2014a). Within a reliability framework, Favier and others (Reference Favier, Bertrand, Eckert and Naaim2014a) 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, Reference Bertrand, Naaim and Brun2010). 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, Reference Favier, Bertrand, Eckert and Naaim2014a), 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 (Reference Favier, Bertrand, Eckert and Naaim2014a). 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(Pressurehouse⩾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(Pressurehouse ⩾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 (Reference Teich and Bebi2009) 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 (Reference Mainieri2020) and Zgheib and others (Reference Zgheib2022) 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, Reference Granet-Abisset1991) and in the European Alps in general (MacDonald and others, Reference MacDonald2000) 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, Reference Zgheib2022). 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, Reference Eckert2018), 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, Reference Teich2019; Stritih and others, Reference Stritih, Bebi, Rossi and Grêt-Regamey2021; Caduff and others, Reference Caduff, Brožová, Kupferschmid, Krumm and Bebi2022). In turn, forest disturbances are sensitive to climate changes (Maroschek and others, Reference Maroschek, Rammer and Lexer2015; Seidl and others, Reference Seidl2017) that alters their frequency, intensity and duration (Dale and others, Reference Dale2001). 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, Reference Faug, Gauer, Lied and Naaim2008; Eckert and others, Reference Eckert2012; Favier and others, Reference Favier2022). 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, Reference Teich and Bebi2009). 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, Reference Poratelli, Accastello, Freppaz and Brun2020) 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, Reference Barbolini, Gruber, Keylock, Naaim and Savi2000; Heredia and others, Reference Heredia, Prieur and Eckert2021, Reference Heredia, Prieur and Eckert2022). 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, Reference Salm1978; Viglietti and others, Reference Viglietti, Letey, Motta, Maggioni and Freppaz2010) (location in the starting zone of the avalanche), or decelerate a flowing avalanche (Anderson and McClung, Reference Anderson and McClung2012) (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, Reference Brang2006) but noticeable lateral spread (Christen and others, Reference Christen, Bartelt and Kowalski2010). 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, Reference Takeuchi, Nishimura and Patra2018). 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, Reference Védrine, Li and Gaume2021). 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, Reference Moos2018). 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, Reference Stritih, Bebi and Grêt-Regamey2019) 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, Reference Eckert, Keylock, Castebrunet, Lavigne and Naaim2013; Ballesteros-Cánovas and others, Reference Ballesteros-Cánovas, Trappmann, Madrigal-González, Eckert and Stoffel2018; Hock and others, Reference Hock2019; Giacona and others, Reference Giacona2021; Schlögl and others, Reference Schlögl, Fuchs, Scheidl and Heiser2021). 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, Reference Zgheib, Giacona, Granet-Abisset, Morin and Eckert2020, Reference Zgheib2022), 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.
Acknowledgements
This work is supported by the French National Research Agency (i) Within the CDP-Trajectories frameworkin the framework of the ‘Investissements d'avenir’ program (ANR-15-IDEX-02), (ii) through the statistical modeling for the assessment and mitigation of mountain risks in a changing environment – SMARTEN program (ANR-20-Tremplin-ERC8-0001). The authors are grateful to the numerous people from ONF-RTM that contributed to the EPA data collection and to the referees and editor whose comments helped produce a better paper. INRAE/ETNA and CNRM/CEN are members of Labex OSUG.
APPENDIX A. A Data collected for mapping forest cover evolution
Table 2 introduces the set of aerial photographs and old maps used for the analysis of the evolution of the forest fraction in the Ravin de Cête-Belle avalanche path.
Table 2. Source data for mapping forest cover evolution at Ravin de Côte-Belle avalanche path from 1825 to 2017
APPENDIX B. Frequentist inference of the mixture model for the release position
According to the topography of the Ravin de Côte-Belle avalanche path, x start is modeled as a binomial mixture of two Beta distributions (Eqn (5)). The five parameters of the mixture, α 1, α 2, β 1, β 2 and p, are estimated using the method of moments as follows:
where N is the total number of avalanches considered for the calibration of the magnitude model (i.e. 17 avalanches at the study site), and, among these, n 1 is the number of avalanches released from zone I (i.e. 13 avalanches at the study site).
Inversion of the standard formula for the mean and variance of a Beta distribution as a function of its parameters leads:
where $\bar {X_{1}}$
, $\bar {X_{2}}$
, v 1, v 2 are, respectively, the empirical mean and variance of normalized release positions in zones I and II calculated as follows:
where n j is the number of avalanches released in the specific release zone j (n 1 = 13 and n 2 = 4 at the study site).
APPENDIX C. Fragility curves for RC buildings exposed to snow avalanches
Favier and others (Reference Favier, Bertrand, Eckert and Naaim2014a) used a reliability analysis to assess fragility curves for RC buildings exposed to snow avalanches, considering four failure limit states of the RC wall facing the avalanche and ten different building configurations (Table 3). The four limit states (q ELS, q ULS, q ALS, q YLT) are defined in the stress–displacement graph in Fig. 11. The loading pressure is the sole avalanche magnitude variable considered. The 40 resulting fragility curves are those of Fig. 12.
Fig. 11. Generic stress–displacement representation of an RC wall subject to quasi-static avalanche loading, with loading pressure as the sole stress variable considered. The diagram highlights the four limit states considered: ELS (elastic limit state), ULS (ultimate limit state), ALS (accidental limit state), YLT (yield line theory, collapse of the building). Each of them leads a specific fragility curve for a given building type (Favier and others, Reference Favier, Bertrand, Eckert and Naaim2014a).
Fig. 12. Forty vulnerability relationships for RC buildings subject to quasi-static avalanche loading considered in this study. Each of them corresponds to one of the four limit states and to one of the ten boundary conditions and was obtained using a reliability analysis (Favier and others, Reference Favier, Bertrand, Eckert and Naaim2014a).
Table 3. Ten RC building types considered, defined by the boundary conditions applying to the wall facing the avalanche
APPENDIX D. Influence of the truncation on the distribution of μ
According to Eqn 16, μ is sampled in a truncated Gaussian distribution. Table 4 provides the mean and SD of the truncated distribution as well as the rejection rate (percentage of negative values in the non-truncated distribution) as a function of the forest fraction.
Table 4. Characteristics of the samples of μ as a function of the forest fraction
APPENDIX E. Avalanche dataset used for the calibration of the statistical–dynamical model
Table 5 provides the details of the 17 avalanches recorded in the EPA database (Bourova and others, Reference Bourova2016) for the Ravin de Cête-Belle avalanche path that were used for the calibration of the magnitude component of the statistical–dynamical model. Elevations were converted into positions on the 1-D topography of the path.
Table 5. Avalanche dataset used for the calibration of the statistical dynamical model
Open access
