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

Taline Zgheib, Florie Giacona, Samuel Morin, Anne-Marie Granet-Abisset, Philomène Favier, Nicolas Eckert

► To cite this version:

Taline Zgheib, Florie Giacona, Samuel Morin, Anne-Marie Granet-Abisset, Philomène Favier, et al. Diachronic quantitative snow avalanche risk assessment as a function of forest cover changes. Journal of Glaciology, 2022, pp.1-19. 10.1017/jog.2022.103. hal-03924185

HAL Id: hal-03924185

https://hal.science/hal-03924185

Submitted on 16 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

**Cite this article:** Zgheib T, Giacona F, Morin S, Granet-Abisset A-M, Favier P, Eckert N (2022). Diachronic quantitative snow avalanche risk assessment as a function of forest cover changes. *Journal of Glaciology* 1–19. https://doi.org/10.1017/jog.2022.103

Received: 29 November 2021

Revised: 27 August 2022

Accepted: 26 September 2022

**Keywords:**

Avalanches; glaciological natural hazards;
snow

**Author for correspondence:**

Taline Zgheib,

E-mail: taline.zgheib@inrae.fr

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

Taline Zgheib¹, Florie Giacona¹, Samuel Morin² ID, Anne-Marie Granet-Abisset³,
Philoméne Favier⁴ and Nicolas Eckert¹ ID

¹INRAE, ETNA, University Grenoble Alpes, 38000 Grenoble, France; ²CNRS, CNRM, Centre d'Etudes de la Neige, University Grenoble Alpes, Université de Toulouse, Météo-France, Grenoble, France; ³UMR CNRS 5190 Laboratoire de Recherche Historique Rhône-Alpes (LARHRA), University Grenoble Alpes, Saint-Martin-d'Hères, France and ⁴INRAE, LESSEM, University Grenoble Alpes, 38000 Grenoble, France

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 $\mu$ and a turbulent term depending on the squared velocity and on the inverse of a coefficient $\xi$ (Voellmy, 1955). The static friction coefficient $\mu$ 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 $\xi$ 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

© The Author(s), 2022. Published by
Cambridge University Press on behalf of The
International Glaciological Society. This is an
Open Access article, distributed under the
terms of the Creative Commons Attribution
licence (http://creativecommons.org/licenses/
by/4.0/), which permits unrestricted re-use,
distribution and reproduction, provided the
original article is properly cited.

<figure><img src="image_2.png" /></figure>
land cover properties (Salm and others, 1990; Ancey and others, 2003). This explains why, in frictional approaches, more often than not, $\xi$ is lowered to values in the order of $400 \text{ 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, 2001). However, some research also considered a concomitant slight increase in $\mu$, with a $\Delta\mu$ 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 stakeholders. 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 $\mu$ and $\xi$ 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 \text{ 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.

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
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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
<figure><img src="image_4.png" /><figcaption>**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².</figcaption></figure>

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:

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

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

where $v$ is the flow velocity, $h$ is the flow depth, $\phi$ 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):

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

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

The full stochastic model (Fig. 3) representing the variability of avalanche events at the study site is noted $p(y, a|\theta_M, \lambda)$. 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):

$$ p(y, a|\theta_M, \lambda) = p(y|\theta_M)p(a|\lambda). \qquad (4) $$

Avalanche frequency is modeled as a Poisson-distributed process, with a scalar parameter $\theta_F = \lambda$ representing the mean annual avalanche number i.e. $a|\lambda \sim P(\lambda)$, necessary for the computation

of return periods. The magnitude model specified below is more complex, with 13 parameters $\theta_M = (\alpha_1, \alpha_2, \beta_1, \beta_2, p, b_1, b_2, \sigma_h, c, d, e, \sigma, \xi)$. Also, the magnitude model evaluates, for each avalanche, the latent friction $\mu$ 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):

$$ \begin{aligned} x_{\text{start}_i} = & P(x_{\text{start}_i} \in [x_{\text{min}_1}, x_{\text{max}_1}])x_{\text{start}_1} + \\ & (1 - P(x_{\text{start}_i} \in [x_{\text{min}_2}, x_{\text{max}_2}]))x_{\text{start}_2}, \\ & P(x_{\text{start}_i} \in [x_{\text{min}_1}, x_{\text{max}_1}]) \sim B(p) \end{aligned} \qquad (5) $$

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

$$ \begin{aligned} x_{\text{start}_{\text{norm}_1}} &= \frac{x_{\text{start}_1} - x_{\text{min}_1}}{x_{\text{max}_1} - x_{\text{min}_1}} 1(x_{\text{start}_1} \in [x_{\text{min}_1}, x_{\text{max}_1}]) \sim \text{Beta}(\alpha_1, \beta_1), \\ x_{\text{start}_{\text{norm}_2}} &= \frac{x_{\text{start}_2} - x_{\text{min}_2}}{x_{\text{max}_2} - x_{\text{min}_2}} 1(x_{\text{start}_2} \in [x_{\text{min}_2}, x_{\text{max}_2}]) \sim \text{Beta}(\alpha_2, \beta_2), \end{aligned} \qquad (6) $$

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 $\alpha_1, \alpha_2, \beta_1$ and $\beta_2$ the parameters of the corresponding Beta distributions, and 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 $\sigma_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:

$$ h_{\text{start}}|b_1, b_2, \sigma_h, x_{\text{start}} \sim \text{Gamma}\left(\frac{1}{\sigma_h^2}(b_1 + b_2x_{\text{start}_{\text{norm}}})^2, \frac{1}{\sigma_h^2}(b_1 + b_2x_{\text{start}_{\text{norm}}})\right), \qquad (7) $$

where $x_{start_{norm}} = (x_{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_{start_{norm}}$ and the mean release depth $h_{start}$, the friction coefficient $\mu$ 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):

$$ \mu|c, d, e, \sigma, x_{\text{start}}, h_{\text{start}} \sim N(c + dx_{\text{start}_{\text{norm}}} + eh_{\text{start}}, \sigma). \qquad (8) $$

The parameters $c, d$ and $e$ represent the dependency of $\mu$ on the release abscissa and mean release depth, with a constant dispersion $\sigma$ 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 $\mu$. Therefore the distribution of the $\mu$ and $h_{start}$ varies according to the release altitude. Alternatively, the velocity-dependent friction coefficient $\xi$ is modeled as a parameter in the strict statistical sense of the term. Both $\xi$ and $\mu_p, i \in [1, n]$ where $n$ is the data sample size of fully documented events are estimated from the data.
<figure><img src="image_5.png" /><figcaption>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: $\theta_M = (\alpha_1, \alpha_2, \beta_1, \beta_2, p, b_1, b_2, \sigma_h, c, d, e, \sigma, \xi)$ and $\theta_F = \lambda$. (2) Simulation of the forest integration. Here three cases are considered: case I representing the dependency on $\mu$, case II on $\xi$ and case III on both $\mu$ and $\xi$. (3) Risk assessment and (4) risk calculation as a function of forest cover and building technology.</figcaption></figure>

Small Gaussian differences between the observed runout distances $x_{stop, data}$ and computed runout distances $x_{stop}$ are postulated:

$$x_{stop, data} | \sigma_{num}, x_{stop} \sim N(x_{stop}, \sigma_{num}). \quad (9)$$

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 $\sigma_{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_\lambda, b_\lambda$) prior for the parameter $\lambda$, the posterior distribution of $\lambda$ remains Gamma-distributed: Gamma($a_\lambda + T_{obs}, b_\lambda + N_{obs}$), where the pair ($a_\lambda, b_\lambda$) 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}, \alpha_1, \alpha_2, \beta_1$ and $\beta_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:

$$\begin{aligned} p(\theta, \mu, x_{stop} | \text{data}, \sigma_{num}) &\propto \pi(\theta) \\ &\times p(h_{start}, x_{stop, data} | \theta, \mu, x_{start}, x_{stop}, \sigma_{num}) \\ &\times p(\mu, x_{stop} | \theta, h_{start}, x_{start}, x_{stop, data}, \sigma_{num}) \end{aligned} \quad (10)$$

where $\theta = (b_1, b_2, \sigma_h, c, d, e, \sigma, \xi)$, $\pi(\theta)$ 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).

## 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
<table><caption>**Table 1.** Parameter estimates of the statistical dynamical model</caption><thead><tr><th></th><th>Prior</th><th>Posterior mean/point estimate</th><th>SD</th><th>2.5%</th><th>97.5%</th></tr></thead><tbody><tr><td>b<sub>1</sub></td><td>b<sub>1</sub> ~ N(3, 1)</td><td>1.56</td><td>0.22</td><td>1.25</td><td>2.16</td></tr><tr><td>b<sub>2</sub></td><td>b<sub>2</sub> ~ N(0, 1)</td><td>-0.54</td><td>0.27</td><td>-1.01</td><td>-0.04</td></tr><tr><td>σ<sub>h</sub></td><td>σ<sub>h</sub> ~ Gamma(5, 10)</td><td>0.59</td><td>0.23</td><td>0.31</td><td>1.2</td></tr><tr><td>c</td><td>c ~ N(0.5, 0.2)</td><td>0.34</td><td>0.05</td><td>0.24</td><td>0.44</td></tr><tr><td>d</td><td>d ~ N(0, 0.25)</td><td>-0.02</td><td>0.05</td><td>-0.12</td><td>0.086</td></tr><tr><td>e</td><td>e ~ N(0, 0.125)</td><td>-0.02</td><td>0.015</td><td>-0.06</td><td>-0.001</td></tr><tr><td>σ</td><td>σ ~ Gamma(1, 0.03)</td><td>0.16</td><td>0.04</td><td>0.09</td><td>0.25</td></tr><tr><td>ξ</td><td>ξ ~ N(1300, 100)</td><td>2297</td><td>13.5</td><td>2275</td><td>2323</td></tr><tr><td>α<sub>1</sub></td><td></td><td>α<sub>1</sub> = 0.28</td><td></td><td></td><td></td></tr><tr><td>α<sub>2</sub></td><td></td><td>α<sub>2</sub> = 0.02</td><td></td><td></td><td></td></tr><tr><td>β<sub>1</sub></td><td></td><td>β<sub>1</sub> = 0.67</td><td></td><td></td><td></td></tr><tr><td>β<sub>2</sub></td><td></td><td>β<sub>2</sub> = 0.01</td><td></td><td></td><td></td></tr><tr><td>p</td><td></td><td>p = 13/17</td><td></td><td></td><td></td></tr><tr><td>λ</td><td>λ ~ Gamma(0.01, 0.001)</td><td>0.31</td><td>0.07</td><td>0.17</td><td>0.45</td></tr></tbody></table>

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).

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_{start_{norm}} = (x_{start} - x_{min1})/(x_{max2} - x_{min1})$, where $x_{start}$ is the simulated release abscissa depending on its location on the path, $x_{max2}$, the maximal abscissa of zone II and $x_{min1}$, 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 $\mu$ and writes as follows:

$$ 
\begin{aligned}
p(y|\hat{\theta}_M) &= \\ 
& \int p(x_{start}|\hat{\alpha}_1, \hat{\alpha}_2, \hat{\beta}_1, \hat{\beta}_2, \hat{p}) \\ 
& p(h_{start}|\hat{b}_1, \hat{b}_2, \hat{\sigma}_h, x_{start})p(x_{stop}|x_{start}, h_{start}, \mu, \xi)d\mu,
\end{aligned} 
\quad (11) 
$$

where “.” classically denotes a statistical estimate. This simulation strategy leads, for instance, the joint distribution of the input variables ($x_{start}$, $h_{start}$, $\mu|\hat{\theta}_M$) and of the output variables $p(x_{stop}, v_{xt}, h_{xt}|\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).

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

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

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

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

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:

$$ CI_{\alpha} = \hat{F}(x_{stop}) \pm q_{N_{\alpha_c}} \sqrt{\frac{\hat{F}(x_{stop})(1 - \hat{F}(x_{stop}))}{n}}, \quad (14) $$

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 $\alpha_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:

$$ P = C_x \frac{1}{2} \rho v^2, \quad (15) $$

where $C_x$ is the dimensionless coefficient of resistance i.e. drag coefficient, $\rho$ 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 $\rho = 300 \text{ 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 $\xi$ significantly, and increasing the value of $\mu$ 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 $\mu$ (e.g. Dent and Lang, 1980; Borstad and McClung, 2009), which has been recently confirmed by more systematic
<figure><img src="image_6.png" /><figcaption>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 $\mu$, (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.</figcaption></figure>

sensitivity analyses by Heredia and others (2021, 2022). Similarly, Teich and others (2012b) showed that focusing on the turbulent friction $\xi$ only when investigating how forest cover impacts avalanche dynamics is not enough since the effects on runout distances of sole slight decrease of $\xi$ remain too limited. Hence, we consider that changes in land cover can affect both $\xi$ and $\mu$ 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 $\mu$ 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 $\mu$ based on the value of the forest fraction $f_k$ relative to $\bar{f}$. This uses the fact that $\mu$ is a latent variable in our statistical-dynamical model. A linear dependency between the mean of the distribution of $\mu$ and $f_k$ is chosen, similar to the one between the mean of $\mu$, release abscissa and mean release depth. To avoid the small number of unphysical negative values of $\mu$ that may result from this formulation with low values of $f_k$, we simply add a positivity constraint and sample the static friction coefficient $\mu$ 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):

$$ \mu \sim N(c + dx_{\text{start,norm}} + eh + g(f_k - \bar{f}), \sigma)1(\mu > 0), \quad (16) $$

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_{\text{start}}$ and $h$ are neglected, the minimal value of the mean of $\mu$, 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 $\mu$ are ~0.15. Within the equation $\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 $\xi$ decays as a power law with an increasing forest fraction:

$$ \xi = \hat{\xi}b(f_k - \bar{f}) \quad (17) $$

where $\hat{\xi}$ is the posterior estimate resulting from the calibration step. This relationship uses the fact that $\xi$ is a parameter in the strict statistical sense of the term. It is considered that, moving from deforested to forested terrain, $\xi$ 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 $\mu$ and $\xi$ 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:

$$r_z = \lambda \int p(y)V_z(y)dy, \quad (18)$$

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:

$$r_z(x_b) = \lambda \int p(P|x_b \le x_{\text{stop}})p(x_b \le x_{\text{stop}}) \times V_z(P)dP. \quad (19)$$

Here, the avalanche magnitude distribution corresponds to the joint distribution of runout distances $x_{\text{stop}}$ and pressures $P$ i.e. $p(P, x_{\text{stop}}) = p(P|x_b \le x_{\text{stop}})p(x_b \le x_{\text{stop}})$ where $p(P|x_b \le x_{\text{stop}})$ is the pressure distribution at abscissa $x_b$ knowing that $x_b$ has been reached by an avalanche and $p(x_b \le x_{\text{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_{\text{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 $\mu$ and the turbulent friction $\xi$), followed very closely by case I (forest fraction acts on $\mu$), whereas modification is much lower in case II (the forest fraction is linked to the velocity-dependent friction $\xi$). 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_{\text{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 1840 m than for abscissas >1840 m (Fig. 5).

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_{\text{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_{\text{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_{\text{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 $\ge 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.

### 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 $\ge 30$ kPa at $x_{\text{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(\text{Pressure}_{\text{house}}\ge 30$ kPa) at $x_{\text{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(\text{Pressure}_{\text{house}}\ge 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 $\ge 30$ kPa at $x_{\text{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
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Annual exceedance probability of runout distances for: (a) case I ($f_k$ acts only on the static friction coefficient $\mu$), (b) case II ($f_k$ acts only on the turbulent friction coefficient $\xi$) and (c) case III ($f_k$ acts on both the static friction coefficient $\mu$ and the turbulent friction coefficient $\xi$). 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).</figcaption></figure>

path, high impact pressures are more likely to occur in the propaga-
tion 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 pres-
sures exceeding 250 kPa was 4% at $x = 1000$, 2% at $x = 1840$ m and
1% only at $x = 400$ m (Fig. 7).

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 = ſ = 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.

From the results, it is immediately clear that, for a given build-
ing type and limit state, the risk decreases with an increasing for-
est 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 tem-
poral 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*<sub>*k*</sub> = *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*<sub>*k*</sub> = 0 and *f*<sub>*k*</sub> = 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).

Discussion

Avalanche risk changes with changes in forest cover

This study aimed at assessing the impact of the changes in the for-
est fraction on avalanche hazard and risk to a hypothetical
<figure><img src="image_8.png" /><figcaption>**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 $\mu$), (b) case II ($f_k$ acts only on the turbulent friction coefficient $\xi$) and (c) case III ($f_k$ acts on the static and turbulent friction coefficient $\mu$ and $\xi$ 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).</figcaption></figure>

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} \ge 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 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
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Annual exceedance probability of impact pressure modeled for the three cases I ($f_k$ acts only on the static friction coefficient $\mu$), II ($f_k$ acts only on the turbulent friction coefficient $\zeta$) 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.</figcaption></figure>

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 $\le 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
<figure><img src="image_10.png" /><figcaption>**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 $\mu$), II ($f_k$ acts only on the turbulent friction coefficient $\xi$) and III ($f_k$ acts on both the static and turbulent friction coefficient $\mu$ and $\xi$ 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.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

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

<figure><img src="image_12.png" /><figcaption>**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 $\mu$) and (b) case II ($f_k$ acts only on the turbulent friction coefficient $\xi$). g and b are the parameters representing the dependency of $\mu$ and $\xi$ 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).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**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 $\mu$) and II ($f_k$ acts only on the turbulent friction coefficient $\xi$). $g$ and $b$ are the parameters representing the dependency of $\mu$ and $\xi$ 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$).</figcaption></figure>

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 $\mu$ 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 $\xi$ (case II). This relates to the high sensitivity to changes in $\mu$ 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 $\mu$ and modification of the parameter $\xi$ with an exponential dependency on $f_k$, as well as with physical and empirical knowledge regarding current values of $\mu$ and $\xi$

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 $\mu$, 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 $\mu$ 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 $\mu$ 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 $\mu$ and $\xi$ 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 (Stritth 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.

**Acknowledgements.** This work is supported by the French National Research Agency (i) Within the CDP-Trajectories framework in 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.

## References

* Ancey C, Gervasoni C and Meunier M (2004) Computing extreme avalanches. *Cold Regions Science and Technology* 39(2-3), 161-180.

* Ancey C, Meunier M and Richard D (2003) Inverse problem in avalanche dynamics models. *Water Resources Research* 39(4), 1-13. doi: 10.1029/2002WR001749.

* Anderson G and McClung D (2012) Snow avalanche penetration into mature forest from timber-harvested terrain. *Canadian Geotechnical Journal* 49(4), 477-484. doi:10.1139/t2012-018

* Arnalds P, Jónasson K and Sigurosson S (2004) Avalanche hazard zoning in Iceland based on individual risk. *Annals of Glaciology* 38, 285-290. doi:10.3189/172756404781814816
**Ballesteros-Cánovas JA, Trappmann D, Madrigal-González J, Eckert N and Stoffel M** (2018) Climate warming enhances snow avalanche risk in the Western Himalayas. *Proceedings of the National Academy of Sciences* **115**(13), 3410–3415.

**Barbolini M, Cappabianca F and Savi F** (2004) Risk assessment in avalanche-prone areas. *Annals of Glaciology* **38**, 115–122. doi:10.3189/172756404781815103

**Barbolini M, Gruber U, Keylock CJ, Naaim M and Savi F** (2000) Application of statistical and hydraulic-continuum dense-snow avalanche models to five real European sites. *Cold Regions Science and Technology* **31**(2), 133–149. doi:10.1016/S0165-232X(00)00008-2

**Bartelt P and Stöckli V** (2001) The influence of tree and branch fracture, overturning and debris entrainment on snow avalanche flow. *Annals of Glaciology* **32**, 209–216. doi:10.3189/172756401781819544

**Bathurst JC, Fahey B, Iroumé A and Jones J** (2020) Forests and floods: using field evidence to reconcile analysis methods. *Hydrological Processes* **34**, 3295–3310. doi:10.1002/hyp.13802

**Bebi PSEP** and 9 others (2017) Changes of forest cover and disturbance regimes in the mountain forests of the alps. *Forest Ecology and Management* **388**, 43–56. Ecology of Mountain Forest Ecosystems in Europe. doi:10.1016/j.foreco.2016.10.028

**Bebi P, Kulakowski D and Rixen C** (2009) Snow avalanche disturbances in forest ecosystems – state of research and implications for management. *Forest ecology and Management* **257**(9), 1883–1892. doi:10.1016/j.foreco.2009.01.050.

**Bertrand D, Naaim M and Brun M** (2010) Physical vulnerability of reinforced concrete buildings impacted by snow avalanches. *Natural Hazards and Earth System Sciences* **10**(7), 1531. doi:10.5194/nhess-10-1531-2010

**Borstad CP and McClung DM** (2009) Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events. *Canadian Geotechnical Journal* **46**(9), 1024–1033. doi:10.1139/T09-042

**Bourova E and 6 others** (2016) A new web-based system to improve the monitoring of snow avalanche hazard in France. *Natural Hazards and Earth System Sciences* **16**(5), 1205–1216. doi:10.5194/nhess-16-1205-2016

**Brang P and 5 others** (2006) Management of protection forests in the European Alps: an overview. *Forest Snow and Landscape Research* **80**(1), 23–44.

**Brang P, Schönenberger W, Ott E and Gardner B** (2001) Forests as protection from natural hazards. In: The Forests Handbook, Vol. 2, pp. 53–81. doi:10.1002/9780470757079.ch3

**Bründl M, Romang HE, Bischof N and Rheinberger CM** (2009) The risk concept and its application in natural hazard risk management in Switzerland. *Natural Hazards and Earth System Sciences* **9**(3), 801.

**Caduff ME, Brožová N, Kupferschmid AD, Krumm F and Bebi P** (2022) How large-scale bark beetle infestations influence the protective effects of forest stands against avalanches: a case study in the Swiss Alps. *Forest Ecology and Management* **514**, 120201. doi:10.1016/j.foreco.2022.120201

**Cappabianca F, Barbolini M and Natale L** (2008) Snow avalanche risk assessment and mapping: a new method based on a combination of statistical analysis, avalanche dynamics simulation and empirically-based vulnerability relations integrated in a GIS platform. *Cold Regions Science and Technology* **54**(3), 193–205. doi:10.1016/j.coldregions.2008.06.005

**Christen M, Bartelt P and Kowalski J** (2010) Back calculation of the In den Arenen avalanche with RAMMS: interpretation of model results. *Annals of Glaciology* **51**(54), 161–168. doi:10.3189/17275640791386553

**Corona C and 5 others** (2013) Seven centuries of avalanche activity at Echalp (Queyas massif, Southern French Alps) as inferred from tree rings. *The Holocene* **23**(2), 292–304.

**Dale VH and 12 others** (2001) Climate change and forest disturbances: climate change can affect forests by altering the frequency, intensity, duration, and timing of fire, drought, introduced species, insect and pathogen outbreaks, hurricanes, windstorms, ice storms, or landslides. *Bioscience* **51**(9), 723–734. doi:10.1641/0006-3568(2001)051[0723:CCAFD]2.0.CO;2

**de Bouchard d'Aubeterre G**, 2019 Tree-ring reconstruction of snow avalanche activity: does avalanche path selection matter?. *Science of the Total Environment* **684**, 496–508.

**Delage P** (2003) Risk in civil engineering: from natural to man-made hazards. In *France–Stanford Conference on 'Risk issues in contemporary science and engineering', Stanford*.

**Dent JD and Lang TE** (1980) Modeling of snow flow. *Journal of Glaciology* **26**(94), 131–140. doi:10.3189/S0022143000010674

**De Quervain M** (1978) Wald und lawinen. In *Proceedings of the IUFRO Seminar Mountain Forests and Avalanches, Davos, Switzerland*, pp. 219–231.

**Dupire S and 6 others** (2016) The protective effect of forests against rockfalls across the French Alps: influence of forest diversity. *Forest Ecology and Management* **382**, 269–279. doi:10.1016/j.foreco.2016.10.020

**Durand Y and 5 others** (2009a) Reanalysis of 47 years of climate in the French Alps (1958–2005): climatology and trends for snow cover. *Journal of Applied Meteorology and Climatology* **48**(12), 2487–2512.

**Durand Y and 5 others** (2009b) Reanalysis of 44 yr of climate in the French Alps (1958–2002): methodology, model validation, climatology, and trends for air temperature and precipitation. *Journal of Applied Meteorology and Climatology* **48**(3), 429–449.

**Eckert N and others** (2010a) Cross-comparison of meteorological and avalanche data for characterising avalanche cycles: the example of December 2008 in the eastern part of the French Alps. *Cold Regions Science and Technology* **64**(2), 119–136.

**Eckert N and 6 others** (2012) Quantitative risk and optimal design approaches in the snow avalanche field: review and extensions. *Cold Regions Science and Technology* **79-80**, 1–19. doi:10.1016/j.coldregions.2012.03.003

**Eckert N and 7 others** (2018) Reflexions on the basis of legal zoning for recurrent mountain risks [Repenser les fondements du zonage réglementaire des risques en montagne]. *Houille Blanche* **2**, 38–67. doi:10.1051/lhb/2018019

**Eckert N, Keylock C, Castebronet H, Lavigne A and Naaim M** (2013) Temporal trends in avalanche activity in the French Alps and subregions: from occurrences and runout altitudes to unsteady return periods. *Journal of Glaciology* **59**(213), 93–114.

**Eckert N, Naaim M and Parent E** (2010b) Long-term avalanche hazard assessment with a Bayesian depth-averaged propagation model. *Journal of Glaciology* **56**(198), 563–586. doi:10.3189/002214310793146331

**Eckert N, Parent É, Faug T and Naaim M** (2008) Optimal design under uncertainty of a passive defense structure against snow avalanches: from a general Bayesian framework to a simple analytical model. *Natural Hazards and Earth System Sciences* **8**(5), 1067–1081.

**Eckert N, Parent E, Faug T and Naaim M** (2009) Bayesian optimal design of an avalanche dam using a multivariate numerical avalanche model. *Stochastic Environmental Research and Risk Assessment* **23**(8), 1123–1141. doi:10.1007/s00477-008-0287-6

**Eckert N, Parent E and Richard D** (2007) Revisiting statistical-topographical methods for avalanche predetermination: Bayesian modelling for runout distance predictive distribution. *Cold Regions Science and Technology* **49**(1), 88–107.

**Farvacque M and 5 others** (2019) How is rockfall risk impacted by land-use and land-cover changes? Insights from the French Alps. *Global and Planetary Change* **174**, 138–152. doi:10.1016/j.gloplacha.2019.01.009

**Faug T, Gauer P, Lied K and Naaim M** (2008) Overrun length of avalanches overtopping catching dams: cross-comparison of small-scale laboratory experiments and observations from full-scale avalanches. *Journal of Geophysical Research: Earth Surface* **113**(F3009), 1–17. doi: 10.1029/2007JF000854.

**Favier P and 6 others** (2022) A framework to account for structural damage, functional efficiency and reparation costs within the optimal design of countermeasures: application to snow avalanche risk mitigation. *Cold Regions Science and Technology* **199**, 103559.

**Favier P, Bertrand D, Eckert N and Naaim M** (2014a) A reliability assessment of physical vulnerability of reinforced concrete walls loaded by snow avalanches. *Natural Hazards and Earth System Sciences* **14**(3), 689–704. doi:10.5194/nhess-14-689-2014

**Favier P, Eckert N, Bertrand D and Naaim M** (2014b) Sensitivity of avalanche risk to vulnerability relations. *Cold Regions Science and Technology* **108**, 163–177. doi:10.1016/j.coldregions.2014.08.009

**Favier P, Eckert N, Faug T, Bertrand D and Naaim M** (2016) Avalanche risk evaluation and protective dam optimal design using extreme value statistics. *Journal of Glaciology* **62**(234), 725–749. doi:10.1017/jog.2016.64

**Feistl T and 6 others** (2014) Observations and modeling of the braking effect of forests on small and medium avalanches. *Journal of Glaciology* **60**(219), 124–138. doi:10.3189/2014JoG12J05

**Feistl T** (2015) Vegetation Effects on Avalanche Dynamics. Ph.D. thesis, Technische Universität München.

**Fischer JT and 5 others** (2020) Bayesian inference in snow avalanche simulation with r.avaflow. *Geosciences* **10**(5), 191. doi:10.3390/geosciences10050191

**Fischer JT, Kofler A, Fellin W, Granig M and Kleemayr K** (2015) Multivariate parameter optimization for computational snow avalanche simulation. *Journal of Glaciology* **61**(229), 875–888.

**Fuchs S, Bründl M and Stötter J** (2004) Development of avalanche risk between 1950 and 2000 in the municipality of Davos, Switzerland.
García-Hernández C and 5 others (2017) Reforestation and land use change as drivers for a decrease of avalanche damage in mid-latitude mountains (NW Spain). *Global and Planetary Change* **153**, 35–50. doi:10.1016/j.glopla-cha.2017.05.001.

Getzner M, Gutheil-Knopp-Kirchwald G, Kreimer E, Kirchmeir H and Huber M (2017) Gravitational natural hazards: valuing the protective function of alpine forests. *Forest Policy and Economics* **80**, 150–159. doi:10.1016/j.forpol.2017.03.015

Giacona F and 8 others (2018) Avalanche activity and socio-environmental changes leave strong footprints in forested landscapes: a case study in the Vosges medium-high mountain range. *Annals of Glaciology* **59**(77), 111–133.

Giacona F and 7 others (2021) Upslope migration of snow avalanches in a warming climate. *Proceedings of the National Academy of Sciences* **118** (44), e2107306118. doi: 10.1073/pnas.2107306118.

Granet-Abisset AM (1991) *La route réinventée. Quelques éléments sur les migrations des Queyrassins au siècle dernier*, Vol. 19. Persée-Portail des revues scientifiques en SHS.

Grét-Regamey A and Straub D (2006) Spatially explicit avalanche risk assessment linking Bayesian networks to a GIS. *NHESS* **6**(6), 911–926. doi:10.5194/nhess-6-911-2006

Gruber U and Bartelt P (2007) Snow avalanche hazard modelling of large areas using shallow water numerical methods and GIS. *Environmental Modelling & Software* **22**(10), 1472–1481. doi:10.1016/j.envsoft.2007.01.001

Gubler H and Rychetnik J (1991) Effects of forests near the timberline on avalanche formation. In *Symposium at Vienna*, pp. 19–37 (doi: 10.1.1.841.5352).

Harbitz CB, Issler D and Keylock CJ (1998) Conclusions from a recent survey of avalanche computational models. In *Proceedings of the anniversary conference*, Vol. 25, pp. 128–135.

Heredia MB, Prieur C and Eckert N (2021) Nonparametric estimation of aggregated sobol'indices: application to a depth averaged snow avalanche model. *Reliability Engineering & System Safety* **212**, 107422.

Heredia MB, Prieur C and Eckert N (2022) Global sensitivity analysis with aggregated Shapley effects, application to avalanche hazard assessment. *Reliability Engineering & System Safety* **222**, 108420.

Hock R and 9 others (2019) High mountain areas: In: IPCC special report on the ocean and cryosphere in a changing climate.

Keylock CJ, McClung DM and Magnússon MM (1999) Avalanche risk mapping by simulation. *Journal of Glaciology* **45**(150), 303–314. doi:10.3189/s0022143000001805

Kumar P and 30 others (2020) Towards an operationalisation of nature-based solutions for natural hazards. *Science of the Total Environment* **731**, 138855. doi:10.1016/j.scitotenv.2020.138855

Lavigne A, Bel L, Parent E and Eckert N (2012) A model for spatio-temporal clustering using multinomial probit regression: application to avalanche counts. *Environmetrics* **23**(6), 522–534. doi:10.1002/env.2167

Le Roux E, Evin G, Eckert N, Blanchet J and Morin S (2021) Elevation-dependent trends in extreme snowfall in the French Alps from 1959 to 2019. *The Cryosphere* **15**(9), 4335–4356.

MacDonald D and 7 others (2000) Agricultural abandonment in mountain areas of Europe: environmental consequences and policy response. *Journal of Environmental Management* **59**(1), 47–69. doi:10.1006/jema.1999.0335

Mainieri R and 8 others (2020) Impacts of land-cover changes on snow avalanche activity in the French Alps. *Anthropocene* **30**, 100244. doi:10.1016/j.ancene.2020.100244

Maroschek M, Rammer W and LEXER MJ (2015) Using a novel assessment framework to evaluate protective functions and timber production in Austrian mountain forests under climate change. *Regional Environmental Change* **15**(8), 1543–1555. doi:10.1007/s10113-014-0691-z

Mather AS, Fairbairn J and Needle CL (1999) The course and drivers of the forest transition: the case of France. *Journal of Rural Studies* **15**(1), 65–90. doi:10.1016/S0743-0167(98)00023-0

McClung DM (2003) Magnitude and frequency of avalanches in relation to terrain and forest cover. *Arctic, Antarctic, and Alpine Research* **35**(1), 82–90. doi:10.1657/1523-0430(2003)035[0082:MAFOAI]2.0.CO;2

Moos C and 5 others (2018) Ecosystem-based disaster risk reduction in mountains. *Earth-Science Reviews* **177**, 497–513. doi:10.1016/j.earscirev.2017.12.011

Naaim M (1998) *Contribution to Snow Drift and Avalanches Flows Modelling*. Habilitation thesis, University Joseph Fourier, Grenoble, France.

Naaim M and 6 others (2008) Snow avalanche pressure on obstacles. In *Proceedings Whistler 2008 International Snow Science Workshop* September 21–27, 2008, 740.

Naaim M, Durand Y, Eckert N and Chambon G (2013) Dense avalanche friction coefficients: influence of physical properties of snow. *Journal of Glaciology* **59**(216), 771–782.

Naaim M, Naaim-Bouvet F, Faug T and Bouchet A (2004) Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects. *Cold Regions Science and Technology* **39**(2–3), 193–204.

Poratelli F and 5 others (2020b) State-of-the-art on ecosystem-based solutions for disaster risk reduction: the case of gravity-driven natural hazards in the alpine region. *International Journal of Disaster Risk Reduction* **51**, 101929. doi:10.1016/j.ijdrr.2020.101929

Poratelli F, Accastello C, Freppaz M and Brun F (2020a) Integrated grey-green management of avalanche risk: economic and ecologic evidences from the Western Italian Alps. *International Journal of Disaster Risk Reduction* **46**, 101502. doi:10.1016/j.ijdrr.2020.101502

Rogger M and 9 others (2017) Land use change impacts on floods at the catchment scale: challenges and opportunities for future research. *Water Resources Research* **53**(7), 5209–5219.

Salm B (1978) Snow forces on forest plants. In *Proc. Seminar on Mountain Forests and Avalanches*, Davos, Switzerland, pp. 157–181.

Salm B, Burkard A and Gubler HU (1990) Calcul des avalanches: une méthode pour le praticien avec des exemples (transl. Christophe Ancey). Technical report, l'Institut fédéral pour l'étude de la neige et des avalanches. (SLF report No. 47).

Schlögl M, Fuchs S, Scheidl C and Heiser M (2021) Trends in torrential flooding in the Austrian Alps: A combination of climate change, exposure dynamics, and mitigation measures. *Climate Risk Management* **32**, 100294.

Seidl R and 9 others (2017) Forest disturbances under climate change. *Nature Climate Change* **7**(6), 395–402. doi:10.1038/nclimate3303

Sovilla B, Schaer M, Kern M and Bartelt P (2008) Impact pressures and flow regimes in dense snow avalanches observed at the Vallée de la Sionne test site. *Journal of Geophysical Research: Earth Surface* **113**(F1), F01–010. doi: 10.3189/002214310793146287.

Stritih A, Bebi P and Grét-Regamey A (2019) Quantifying uncertainties in earth observation-based ecosystem service assessments. *Environmental Modelling & Software* **111**, 300–310. doi:j.envsoft.2018.09.005

Stritih A, Bebi P, Rossi C and Grét-Regamey A (2021) Addressing disturbance risk to mountain forest ecosystem services. *Journal of Environmental Management* **296**, 113188. doi:10.1016/j.jenvman.2021.113188

Takeuchi Y, Nishimura K and Patra A (2018) Observations and numerical simulations of the braking effect of forests on large-scale avalanches. *Annals of Glaciology* **59**, 1–9. doi:10.1017/aog.2018.22

Takeuchi Y, Torita H, Nishimura K and Hirashima H (2011) Study of a large-scale dry slab avalanche and the extent of damage to a cedar forest in the Makunosawa Valley, Myoko, Japan. *Annals of Glaciology* **52**(58), 119–128. doi:10.3189/172756411797252059

Teich M and 5 others (2012b) Avalanche simulations in forested terrain: A framework towards a Bayesian probabilistic model calibration. pp. 628–632.

Teich M and 5 others (2014) Computational snow avalanche simulation in forested terrain. *NHESS* **14**(8), 2233–2248. doi:10.5194/nhess-14-2233-2014

Teich M and 5 others (2019) Effects of bark beetle attacks on forest snow-pack and avalanche formation – implications for protection forest management. *Forest Ecology and Management* **438**, 186–203. doi:10.1016/j.foreco.2019.01.052

Teich M, Bartelt P, Grét-Regamey A and Bebi P (2012a) Snow avalanches in forested terrain: influence of forest parameters, topography, and avalanche characteristics on runout distance. *Arctic, Antarctic, and Alpine Research* **44** (4), 509–519. doi:10.1657/1938-4246-44.4.509

Teich M and Bebi P (2009) Evaluating the benefit of avalanche protection forest with GIS-based risk analyses – a case study in Switzerland. *Forest Ecology and Management* **257**(9), 1910–1919. doi:10.1016/j.foreco.2009.01.046

Védrine L, Li X and Gaume J (2021) Detrainment and braking of snow avalanches interacting with forests. *Natural Hazards and Earth System Sciences Discussions* **22**, 1015–1028.

Vernay M and 5 others (2019) The S2M meteorological and snow cover reanalysis in the French mountainous areas (1958–present), [data set], AERIS.

Viglietti D, Letey S, Motta R, Maggioni M and Freppaz M (2010) Snow avalanche release in forest ecosystems: a case study in the Aosta Valley Region (NW-Italy). *Cold Regions Science and Technology* **64**(2), 167–173. International Snow Science Workshop 2009, Davos. doi:10.1016/j.coldregions.2010.08.007

Voellmy A (1955) Ueber die Zerstöerungskraft von Lawinen Schweizerische Bauzeitung. English version 'On the destructive force of avalanches' translated by Tate RE (1964), ed. US Department of Agriculture Forest Service.
Wohlgemuth T, Schwitter R, Bebi P, Sutter F and Brang P (2017) Post-windthrow management in protection forests of the Swiss Alps. *European Journal of Forest Research* **136**(5), 1029–1040. doi:10.1007/s10342-017-1031-x

Zgheib T and 5 others (2022) Spatio-temporal variability of avalanche risk in the French Alps. *Regional Environmental Change* **22**(1), 1–18.

Zgheib T, Giacona F, Granet-Abisset AM, Morin S and Eckert N (2020) One and a half century of avalanche risk to settlements in the upper Maurienne valley inferred from land cover and socio-environmental changes. *Global Environmental Change* **65**, 102–149. doi:10.1016/j.gloenvcha.2020.102149

## 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><caption>**Table 2.** Source data for mapping forest cover evolution at Ravin de Côte-Belle avalanche path from 1825 to 2017</caption><thead><tr><th>Year</th><th>Data type</th><th>Scale/resolution</th><th>Source</th></tr></thead><tbody><tr><td>1825</td><td>Cadastral map</td><td>1/1250</td><td>Cadastral service of the Hautes-Alpes department</td></tr><tr><td>1948</td><td>Black-and-white aerial photographs</td><td>0.5 m</td><td>IGN</td></tr><tr><td>1980</td><td>Black-and-white aerial photographs</td><td>0.5 m</td><td>IGN</td></tr><tr><td>2017</td><td>Color aerial photographs</td><td>1.5 m</td><td>IGN</td></tr></tbody></table>

IGN is the French National Geographical Institute.

## 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, $\alpha_1, \alpha_2, \beta_1, \beta_2$ and $p$, are estimated using the method of moments as follows:

$$ \hat{p} = P(x_{start_i} \in [x_{min_i}, x_{max_i}]) = \frac{n_1}{N}, \qquad (B.1) $$

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:

$$ \begin{aligned} \hat{\alpha}_1 &= \bar{X}_1 \left( \frac{\bar{X}_1(1-\bar{X}_1)}{v_1} - 1 \right), \\ \hat{\alpha}_2 &= \bar{X}_2 \left( \frac{\bar{X}_2(1-\bar{X}_2)}{v_2} - 1 \right), \\ \hat{\beta}_1 &= (1 - \bar{X}_1) \left( \frac{\bar{X}_1(1-\bar{X}_1)}{v_1} - 1 \right), \\ \hat{\beta}_2 &= (1 - \bar{X}_2) \left( \frac{\bar{X}_2(1-\bar{X}_2)}{v_2} - 1 \right), \end{aligned} \qquad (B.2) $$

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:

$$ \begin{aligned} \bar{X}_j &= \frac{1}{n_j} \sum_{i=1}^{n_j} X_{norm_i}, \\ v_j &= \frac{1}{n_j} \sum_{i=1}^{n_j} (X_{norm_i} - \bar{X}_j)^2, \end{aligned} \qquad (B.3) $$

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 (2014a) 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

<table><caption>**Table 3.** Ten RC building types considered, defined by the boundary conditions applying to the wall facing the avalanche</caption><thead><tr><th>Building type</th><th>Boundary condition</th></tr></thead><tbody><tr><td>I</td><td>Four simply supported edges</td></tr><tr><td>II</td><td>Simply supported on the two large edges, clamped on the two small edges</td></tr><tr><td>III</td><td>Simply supported on one large edge, clamped on the three other edges</td></tr><tr><td>IV</td><td>One free large edge, simply supported on the three other edges</td></tr><tr><td>V</td><td>One free large edge, clamped on the three other edges</td></tr><tr><td>VI</td><td>Clamped on the small edge, simply supported on the three other edges</td></tr><tr><td>VII</td><td>Simply supported side by side, clamped on the two other edges</td></tr><tr><td>VIII</td><td>Four clamped edges</td></tr><tr><td>IX</td><td>One free large edge, one clamped large edge, simply supported on the two small edges</td></tr><tr><td>X</td><td>One free large edge, one simply supported large edge clamped on the two small edges</td></tr></tbody></table>

Each of these conditions, once a limit state is chosen, leads a specific fragility curve usable as input for risk assessment (Favier and others, 2014a).

<figure><img src="image_14.png" /><figcaption>**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, 2014a).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**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, 2014a).</figcaption></figure>
avalanche magnitude variable considered. The 40 resulting fragility curves are those of Fig. 12.

## APPENDIX D. Influence of the truncation on the distribution of $\mu$

According to Eqn 16, $\mu$ 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><caption>**Table 4.** Characteristics of the samples of $\mu$ as a function of the forest fraction</caption><thead><tr><th>f<sub>k</sub></th><th>Number of simulations</th><th>Rejection rate %</th><th>Mean</th><th>SD</th></tr></thead><tbody><tr><td>1</td><td>10 000</td><td>0</td><td>0.6</td><td>0.15</td></tr><tr><td>0.46</td><td>10 000</td><td>0</td><td>0.38</td><td>0.15</td></tr><tr><td>0.35</td><td>10 000</td><td>0</td><td>0.31</td><td>0.14</td></tr><tr><td>0.24</td><td>10 000</td><td>0.6</td><td>0.26</td><td>0.13</td></tr><tr><td>0.16</td><td>10 000</td><td>11</td><td>0.23</td><td>0.12</td></tr><tr><td>0</td><td>10 000</td><td>26</td><td>0.17</td><td>0.11</td></tr></tbody></table>

Samples are simulated with $g = 0.6$. Mean and SD are those of the truncated distribution.

## 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, 2016) 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><caption>**Table 5.** Avalanche dataset used for the calibration of the statistical dynamical model</caption><thead><tr><th>Occurrence date</th><th>Release elevation m a.s.l.</th><th>Runout elevation m a.s.l.</th><th>Estimated snow deposit, length × width × depth m<sup>3</sup></th></tr></thead><tbody><tr><td>2/7/1960</td><td>2400</td><td>1600</td><td>189 × 143 × 2</td></tr><tr><td>4/18/1963</td><td>2400</td><td>1590</td><td>140 × 143 × 1</td></tr><tr><td>3/13/1969</td><td>2200</td><td>1590</td><td>108 × 671 × 1</td></tr><tr><td>1/12/1970</td><td>2200</td><td>1590</td><td>90 × 671 × 1</td></tr><tr><td>3/21/1971</td><td>2450</td><td>1600</td><td>184 × 67 × 2</td></tr><tr><td>2/19/1972</td><td>2450</td><td>1600</td><td>134 × 67 × 1</td></tr><tr><td>1/15/1978</td><td>1900</td><td>1590</td><td>145 × 1211 × 1</td></tr><tr><td>2/19/1978</td><td>1850</td><td>1590</td><td>144 × 1377 × 1</td></tr><tr><td>3/20/1978</td><td>2300</td><td>1590</td><td>174 × 408 × 2</td></tr><tr><td>1/15/1980</td><td>1800</td><td>1590</td><td>77 × 1470 × 1</td></tr><tr><td>1/24/1980</td><td>1800</td><td>1590</td><td>91 × 1470 × 1</td></tr><tr><td>3/4/1993</td><td>2450</td><td>1590</td><td>130 × 67 × 1</td></tr><tr><td>1/13/2004</td><td>2450</td><td>1600</td><td>119 × 67 × 1</td></tr><tr><td>12/16/2008</td><td>2450</td><td>1600</td><td>137 × 67 × 1</td></tr><tr><td>2/27/2010</td><td>2400</td><td>1650</td><td>154 × 143 × 1</td></tr><tr><td>1/8/2018</td><td>2450</td><td>1650</td><td>143 × 67 × 1</td></tr><tr><td>4/12/2018</td><td>2450</td><td>1600</td><td>181 × 67 × 2</td></tr></tbody></table>