# Extreme avalanche cycles: Return levels and probability distributions depending on snow and meteorological conditions

Guillaume Evin, Pascal Dkengne Sielenou, Nicolas Eckert, Philippe Naveau, Pascal Hagenmuller, Samuel Morin

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

## ► To cite this version:

Guillaume Evin, Pascal Dkengne Sielenou, Nicolas Eckert, Philippe Naveau, Pascal Hagenmuller, et al.. Extreme avalanche cycles: Return levels and probability distributions depending on snow and meteorological conditions. Weather and Climate Extremes, 2021, 33, 10.1016/j.wace.2021.100344. hal-03401778

HAL Id: hal-03401778

https://hal.science/hal-03401778

Submitted on 28 Oct 2021

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Extreme avalanche cycles: Return levels and probability distributions depending on snow and meteorological conditions

Guillaume Evin<sup>a,*</sup>, Pascal Dkengne Sielenou<sup>a</sup>, Nicolas Eckert<sup>a</sup>, Philippe Naveau<sup>b</sup>, Pascal Hagenmuller<sup>c</sup>, Samuel Morin<sup>c</sup>

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

## ARTICLE INFO

### Keywords:

- Avalanche cycles

- Extreme value theory

- Discrete distributions

- French alps

## ABSTRACT

Remarkable episodes of avalanche events, so-called snow avalanche cycles, are recurring threats to people and infrastructures in mountainous areas. This study focuses on the hazard assessment of snow avalanche cycles defined by daily occurrence numbers exceeding the 2-year return level. To this aim, extreme value distributions are tailored to account for discrete observations and potential covariates. A comprehensive statistical framework is provided including model fitting, model selection and evaluation, and derivation of quantities of interest such as return levels. In each of the 23 massifs of the French Alps, two discrete generalized Pareto (dGP) models are applied to extreme avalanche cycles extracted from 60 years of daily avalanche activity observations from 1958 to 2018, an unconditional version and a conditional version incorporating snow and meteorological covariates. In the conditional dGP model, the scale parameter is allowed to depend on snow and meteorological conditions from a local reanalysis, leading the corresponding distributions to outperform their unconditional counterparts in about half of the French Alps massifs. Unconditional dGP models provide valuable estimates of high return levels of avalanche numbers. In particular, it is shown that the number of avalanches per path which can be expected on average every 100 and 300 years for the French Alps is approximately equal to 0.25 (roughly one avalanche for four paths) and 0.32 (one avalanche for three paths). As exemplified with the January 2018 Eleanor winter storm, conditional dGP models refine the statistical description of the largest avalanche cycles by providing the information conditional to specific meteorological and snow conditions, with potential applications to avalanche forecasting and climate change impact studies. The same framework could be put to work in other mountain areas and for analyzing extreme counts of various other damaging phenomena.

## 1. Introduction

In the Alpine environment, snow avalanches are an everyday threat resulting in casualties as well as direct and indirect economic losses. As exemplified in December 2008 in southeast France (Eckert et al., 2010b), recurring severe inconvenient consequences of snow avalanches are the inaccessibility to tourism destinations and infrastructures, and public facilities such as schools and hospitals. As shown by the example of February 1999, during most extreme winters, snow avalanches can hold 100,000 tourists hostage in resorts cut off from the outside world. Estimates of the material damage incurred due to snow avalanches during that period are close to 400 million Euros in Switzerland, extrapolated to about 1 billion Euros for all of the Alps (Ammann and Bebi, 2000). No evacuation is possible after a snow avalanche

has been released because of the extreme rapidity of the flow, making anticipation of future avalanche activity crucial both in short- (Morin et al., 2020) and long-term (Eckert et al., 2010c) management. Notably, a correct assessment of the statistical characteristics of the most extreme events is critical due to their potential catastrophic consequences (Schweizer et al., 2009). To this aim, the concept of avalanche cycle is often used to highlight a remarkable cluster of avalanche events at a given spatial scale (the mountain range, the district, etc.) and over a short period of time (typically a few days), and being able to characterize the severity of avalanche cycles in terms of probabilities is valuable.

<footnote><sup>a</sup> Univ. Grenoble Alpes, INRAE, UR ETGR, Grenoble, France</footnote>

<footnote><sup>b</sup> Laboratoire des Sciences du Climat et de l'Environnement, UMR 8212 CEA, CNRS, UVSQ, IPSL and Université Paris-Saclay, CE Orme des Merisiers, Gif-sur-Yvette, France</footnote>

<footnote><sup>c</sup> Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, CEN, 38000 Grenoble, France</footnote>

<footnote>* Corresponding author.
E-mail address: guillaume.evin@inrae.fr (G. Evin).</footnote>
Many geophysical phenomena can be, at the temporal scale of interest, considered as continuous (e.g. temperature). In the case of extreme values of such quantities, a well-established statistical strategy is to take advantage of Extreme Value Theory (EVT) (Embrechts et al., 1997; Coles, 2001; Beirlant et al., 2004; De Haan and Ferreira, 2007). This approach is well tailored for continuous variables and the exceedances (i.e. peaks over threshold, POT) values are classically modeled by a Generalized Pareto (GP) distribution (see, e.g. Davison, 1984), often considering that observed extremes are independently and identically distributed (see, e.g. Serinaldi and Kilsby, 2014, for an application to rainfall extremes). A variant of the POT analysis can be applied using covariates in order to improve the statistical representation of extremes and integrate physical relationships with the quantity of interest (Davison and Smith, 1990). In snow avalanche hazard assessments, the statistical approaches to high return period avalanche evaluation are often based on the EVT framework, even if it is not always explicit. For instance, the runout ratio model (McClung and Lied, 1987) applies Gumbel or a Generalized Extreme Value (GEV) distribution to normalized runouts of extremes avalanches collected over a sample of paths (Keylock, 2005; McClung and Schaerer, 2006) and may be seen as a specific application of the block-maxima GEV approach. Ancey (2012) has discussed the behavior of extreme avalanches with regards to outliers' theory. Eventually, how EVT can be included within a risk-based framework was investigated in Favier et al. (2016).

While it is relatively straightforward to precisely quantify the magnitude/frequency relationship on given snow avalanche paths, quantifying the return period at larger spatial scales corresponding to an avalanche cycle is more difficult. First, a quantitative definition of an avalanche cycle needs to be specified. The existing definitions are sometimes based on the sole number of recorded events (e.g. Eckert et al., 2010b), but more generally on cumulative indices taking into account one or several magnitude variables such as the total mass of avalanche deposits (Birkeland et al., 2001; Laternser and Schneebeli, 2002). This use of cumulative indices has the advantage of lowering the bias induced by the non-observation of small avalanche events, but generally involves estimating the snow mass from the size classes of the avalanches. Second, another specific difficulty is that discrete random variables such as avalanche counts cannot be fitted by classical continuous generalized extreme value distributions as easily as in the continuous case (Anderson, 1980; Leadbetter et al., 1983). Indeed, the maximum term of an integer-valued random variable can be approximated by a continuous extreme value distribution under certain conditions only (Anderson et al., 1997; Nadarajah and Mitov, 2002; Dkengne Sielenou et al., 2016). A dedicated framework based on discretized extreme value distributions should therefore be used, as already employed in other fields such as reliability (Nakagawa and Osaki, 1975; Nakagawa, 1978), dentistry (Krishna and Singh Pundir, 2009), hydrometeorology (Chakraborty and Chakravarty, 2014) or accidentology (Prieto et al., 2014).

Avalanche activity is controlled by permanent and variable factors (International Commission of Snow and Ice, 1981). Permanent factors are related to the mountain landscape (elevation, slope, aspect, configuration of terrain, roughness of the ground, etc.) and variable factors are related to meteorological conditions (snowfall, rainfall, wind, temperature, etc.). Notably, avalanche cycles are generally caused by a severe storm bringing large snowfall accompanied by substantial drifting snow, but strong temperature variations causing snowmelt and/or fluctuations of the freezing level can also be involved. Studying avalanche cycles therefore mainly aims at understanding their relationships with meteorological factors, such as precipitation, temperature, wind effects (Birkeland and Mock, 2001; Höller, 2009). Even if the complex non-linear relationship between avalanche release and these factors makes an exact deterministic prediction of avalanche release out of reach (Schweizer et al., 2003, 2009; Vernay et al., 2015), this control by meteorological and snow conditions of avalanche activity is strong enough for weather forecasting to be the basis of operational

avalanche forecasting, together with meteorological and snow cover observations (Morin et al., 2020). For instance, an increasingly wide and complex range of statistical avalanche forecasting models has been developed to relate snow and meteorological data to avalanche activity (e.g. Obled and Good, 1980; Gassner and Brabec, 2002; Hendriks et al., 2005; Baggi and Schweizer, 2009; Dreier et al., 2016). Operationally, these statistical models are valuable tools supplementing physically-based simulations and expert knowledge for making decisions such as closing ski resorts and/or evacuate the threatened mountain communities when a critical level is reached. However, they generally provide a deterministic avalanche/non-avalanche day classification, which does not come up with a probabilistic information about the severity of an avalanche cycle.

On this basis, the goal of this work is to construct predictive probabilistic models for extreme avalanche cycles providing (i) unconditional return periods, and (ii) statistical relationships between meteorological and snow conditions and extreme daily avalanche numbers. To this aim, we take advantage of the comprehensive and high-quality data available all over the French Alps: 60 years of daily reanalyzed snow and meteorological conditions and avalanche activity observations in the 23 massifs of the French Alps from 1958 to 2018, obtained from the so-called SAFRAN - SURFEX/ISBa - crocus - MEPRA (S2M) modeling chain (Durand et al., 2009a,b; Vernay et al., 2019). Recorded snow avalanche occurrence numbers are available over the same 1958-2018 period thanks to the *Enquête Permanente sur les Avalanches* (EPA, see Bourova et al., 2016). For a given massif, avalanche cycles are defined as daily occurrence numbers exceeding those observed in average once every two years. In each of the 23 massifs, we first fit a discrete version of the generalized Pareto (dGP) distribution (Asadi et al., 2001) to avalanche cycles. We then select the best conditional dGP model having a log-linear varying scale parameter with respect to snow and meteorological covariates provided by S2M reanalyses. The predictive performances of the unconditional and conditional models are evaluated using a dedicated scoring rule. In Section 2, we present the processed datasets, as well as the study area. Section 3 contains our proposed strategy for the statistical modeling of avalanche cycles. Section 4 contains the main results from the application of our approaches to the 23 massifs of the French Alps. Section 5 provides a final discussion and concludes.

## 2. Data

This study focuses on avalanche activity observed on specified paths over the period 1958-2018 in 23 French mountain regions known as massifs (Durand et al., 1999), which is also the spatial scale of the avalanche bulletin (see map in Fig. 1).

### 2.1. Avalanche activity in the French Alps

In the French Alps, daily observed avalanche data are provided by the *Enquête Permanente sur les Avalanches (EPA)* which monitors avalanche events on approximately 3900 paths since the beginning of the 20th century (see Mougin, 1922). Quantitative (runout elevations, deposit volumes, etc.) and qualitative (flow regime, snow quality, etc.) information is recorded for each event. The report quality varies in time and space as a function of local observers (mostly forestry rangers). Another source of uncertainty with regards to natural avalanche activity is that records mainly concern paths that are visible from the valleys, presumably underestimating activity at high elevation. However, the EPA clearly stands among the worldwide longest and most comprehensive records of natural avalanche activity (Bourova et al., 2016). Studies by Jomelli et al. (2007), Eckert et al. (2010a), Eckert et al. (2010b), Lavigne et al. (2012), Eckert et al. (2013) and Lavigne et al. (2015) include further discussions of EPA's record strengths and weaknesses.
<figure><img src="image_5.png" /><figcaption>Fig. 1. The 23 SAFRAN massifs of the French Alps. Numbers between brackets indicate the number of monitored snow avalanche paths in each massif.</figcaption></figure>

In this study, we describe avalanche activity by the daily number of avalanche events recorded in the EPA. An avalanche event may be recorded a few days after the avalanche release. In this case, the exact day of the event estimated by the observer may be approximated. On the one hand, the selection of the events for which the day is known is too restrictive (too many events are lost). On the other hand, in the case where the avalanche could have occurred many days before the observation, the inclusion of the observation could introduce a bias in our analysis (inaccurate number of avalanches for that day, difficulty to relate this event to snow and meteorological conditions). As a compromise, we consider only avalanche events which happened less than three days before their observation. We aggregate the daily numbers of avalanche events at the massif scale, for the 60 winter seasons covering the period 1958 to 2018. This corresponds to 63,994 avalanche events, with 502 to 5,997 avalanches per massif during the entire period. Note that the number of paths covered by the EPA varies substantially from one massif to another, ranging from 35 for the Vercors massif to 297 for the Oisans massif (see Fig. 1).

## 2.2. Snow and meteorological conditions

In this study, snow and meteorological conditions are provided by the S2M modeling chain (Durand et al., 2009a,b; Vernay et al., 2019). The surface area of each massif ranges between 500 km² and 1500 km² and the key assumption regarding snow and meteorological numerical simulations is their spatial homogeneity, i.e. within each massif, meteorological and snowpack properties are assumed to depend only on elevation, slope and aspect (see Durand et al., 1999). The meteorological data are provided by the SAFRAN system (Durand et al., 2009b), a meteorological downscaling and surface analysis system performing an objective analysis of meteorological data available from various observation networks. Snow conditions are simulated, based on meteorological data from SAFRAN, by the Crocus snow cover model, part of the land surface model ISBA (Durand et al., 2009a; Vionnet et al., 2012), for different elevation bands and aspects (North, East, South, and West).

In this study, for each massif, different snow and meteorological data are considered as potential explanatory variables for the avalanche
activity (see Table 1). Daily values spanning the same period as avalanche records (1958–2018) of rainfall and snow precipitation, minimum, mean and maximum daily temperature and wind speed are included. For both solid and liquid precipitation and temperature, 3-day aggregated values are also included. Concerning snow conditions, snow thickness and the thickness of wet snow for the different elevation bands and aspects (North, East, South, and West) are considered. For both meteorological and snow conditions, variables are made available at the elevations of 1800 m, 2400 m and 3000 m when possible (the maximum elevation being less than 2400 m for several massifs, e.g. the Vercors massif).

# 3. Methods

## 3.1. Identification of avalanche cycles

The concepts of return period and return level are commonly used to convey information about the likelihood of extreme events such as floods, earthquakes, volcanic eruptions, avalanche disasters, etc. These concepts are used here to define avalanche cycles for each massif of the French Alps. Avalanche cycles are denoted by $e_t$, the $t$th daily occurrence number strictly greater than a threshold $h$, where $h$ is the 2-year return level of daily avalanche occurrence numbers, computed for each massif. Here, the threshold value $h$ above which an avalanche occurrence number of a given day is considered to be remarkable corresponds to the sample quantile of order $\alpha = 1 - (2 \times 365.25)^{-1} \approx 99.86\%$ of all avalanche occurrence numbers for this massif. This high quantile ensures that $e_t$ corresponds to rare events. For each massif, we thus obtain a sample $e_1, \dots, e_N$ of remarkable daily avalanche occurrence numbers where $N$ is the number of avalanche cycles in the massif. Obviously, we have by construction $N \approx 30$ for all massifs since the studied data spread over 60 years (see Table 3 below) and the return periods of $e_t$ are greater than 2 years.

This definition of avalanche cycles using a percentile has the advantage to be independent from the chosen spatial scale. Also, it is advantageous for operational purposes, because in areas where high avalanche activity is usual, operational services are generally well prepared, whereas a few events can cause considerable problems in areas where activity is usually low. In addition, our definition facilitates the comparison of the main drivers of high activity from one massif to another. Indeed, even in two massifs where the avalanche activity is rather different (say a massif that is usually affected by lots of avalanches and another by very few), we thus compare snow and meteorological drivers that induce a similarly rare avalanche activity. Note eventually that the two-year threshold was chosen as a compromise: sufficiently high to isolate days with a really high avalanche activity and not too high to get “enough” days as remarkable. Yet, our approach can be implemented with other threshold choices.

Fig. 2 represents resulting avalanche cycles $e_t$ for the different massifs. The massif Oisans clearly shows the highest avalanche activity, with up to 80 avalanches which have been observed in one day. Massifs at lower average and maximum elevations (e.g. Chartreuse, Vercors) do not reach $e_t$ higher than 20 avalanches in one day. These numbers can be put in perspective with the total number of paths monitored in each massif, which varies between 35 (Vercors) and 297 (Oisans), as shown in Fig. 1.

## 3.2. Pre-selection and transformation of the datasets of meteorological and snow covariates

For each massif, we first perform a pre-selection and a transformation of the meteorological and snow covariates. Considering the different elevations and aspects, the total number of covariates can be up to 72 (12 × 3 = 36 meteorological variables and 3 × 4 × 3 = 36 snow variables). Most of these meteorological and snow variables are very strongly correlated to each other (see Fig. A.9 in Appendix for an

illustration). In order to avoid too much redundancy in the covariates and numerical instabilities while fitting the statistical models, we apply a first selection based on the cross-correlations between the covariates. Some covariates are thus discarded using a pair-wise absolute correlation cutoff of 0.95. As a result, about half of the covariates are retained in each massif.

In addition, we apply a principal component analysis (PCA) in order to summarize the reduced set of covariates with a smaller number of representative and uncorrelated variables that collectively explain most of the variability in the original covariates. These transformed covariates (i.e. principal components), expressed as a linear transform of the reduced set of covariates, are statistically uncorrelated from each other. They are denoted by $F_j$ for the principal component $j$. Different numbers of $F_j$ are tested thereafter for the modeling of remarkable avalanche numbers. Illustration of resulting transformed covariates are provided in Appendix for the massif Haute-Maurienne (Fig. A.10). For this massif, the interpretation of the first two covariates is not straightforward. The third covariate is related to the thickness of wet snow (TWS) and the fourth covariate is clearly related to 3-day aggregated values of rain and snowfall.

## 3.3. Discrete peak over threshold modeling

### 3.3.1. Unconditional and conditional dGP models

In what follows, we make the assumption that the discrete versions of the classical continuous generalized extreme value distributions provide adequate models for the extreme values of integer-valued random variables such as our massif-scale avalanche counts. Let $E_t$ denote the random variable associated to remarkable daily avalanche occurrence numbers $e_t$. For each massif, because of the discrete nature of $E_t$, the discrete version of the classical continuous generalized Pareto distribution (Asadi et al., 2001; Prieto et al., 2014) is considered. More explicitly, we consider the probability mass function of $E_t$ to be of the form

$$ \mathrm{Pr} \{ E_t = k; \sigma(t), \gamma \} = \left[ 1 + \gamma \left( \frac{k - \lfloor h \rfloor - 1}{\sigma(t)} \right) \right]^{-1/\gamma} - \left[ 1 + \gamma \left( \frac{k - \lfloor h \rfloor}{\sigma(t)} \right) \right]^{-1/\gamma}, $$

where $t = 1, \dots, N$ and $k > \lfloor h \rfloor$. Here, the number of avalanches $k$ is the possible value taken by the random variable $E_t$ for this day $t$ and can only exceed the threshold $h$ rounded to its floor value $\lfloor h \rfloor$. The shape parameter $\gamma \neq 0$ and the scale parameter function $\sigma(t) > 0$ are unknown. As indicated above, two versions of the discrete generalized Pareto (dGP) distribution are considered in this study, and only differ in the definition of the scale parameter function $\sigma(t)$:

1. **Unconditional dGP model:** In this case, the scale parameter functions $\sigma(t)$ corresponds to a unique parameter $\sigma_1$, which needs to be estimated.

2. **Conditional dGP model:** In this second version, the scale parameter varies as a log-linear function of the first principal components $F_j$, that is

$$ \log \sigma(t) = \sigma_1 + \sum_{j=1}^{J} \sigma_{j+1} F_j(t), \quad (1) $$

where the coefficients $\sigma_j$ have to be estimated. Different numbers $J = 1, \dots, 4$ of principal components $F_j(t)$ are considered thereafter (more than 4 covariates not leading to better performances, results not shown). Here, the real number $F_j(t)$ is the coordinate on the $j$th principal component $F_j$ of the vector of original covariates associated to the observed number $e_t$. In other words, $F_j(t)$ is the value taken by the $j$th new covariate at the date where $e_t$ is observed.
<caption>**Table 1**</caption>

SAFRAN meteorological variables and Crocus snow variables used for the study. Snow variables are available for the four aspects (N, S, W, E), all variables are available for the elevations 1800 m, 2400 m and 3000 m.

<table><thead><tr><th>Meteorological variables</th><th>Aspect</th><th>Elevation</th></tr></thead><tbody><tr><td>RAIN: Daily rain precipitation [kg/m²]</td><td></td><td></td></tr><tr><td>RAIN3D: 3-day moving sum of rain precipitation [kg/m²]</td><td></td><td></td></tr><tr><td>SNOW: Daily snow precipitation [kg/m²]</td><td></td><td></td></tr><tr><td>SNOW3D: 3-day moving sum of snow precipitation [kg/m²]</td><td></td><td></td></tr><tr><td>TMIN: Daily minimum temperature [°C]</td><td></td><td>1800 m</td></tr><tr><td>TMIN3D: 3-day moving average of minimum temperature [°C]</td><td></td><td>2400 m</td></tr><tr><td>TMEAN: Daily mean temperature [°C]</td><td>/</td><td>3000 m</td></tr><tr><td>TMEAN3D: 3-day moving average of mean temperature [°C]</td><td></td><td></td></tr><tr><td>TMAX: Daily maximum temperature [°C]</td><td></td><td></td></tr><tr><td>TMAX3D: 3-day moving average of maximum temperature [°C]</td><td></td><td></td></tr><tr><td>WIND: Daily average of wind speed [m/s]</td><td></td><td></td></tr><tr><td>WIND3D: 3-day moving average of wind speed [m/s]</td><td></td><td></td></tr><tr><td colspan="3">Snow variables</td></tr><tr><td>SD: Total snow depth [m]</td><td>N, S, W, E</td><td>1800 m 2400 m</td></tr><tr><td>TWS: Thickness of wet snow at the top of the snowpack [m]</td><td>N, S, W, E</td><td>3000 m</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Box plots of remarkable avalanche numbers in the 23 massifs of the French Alps. The boxes span the inter-quartile range from the 1st to 3rd quartile with the bold vertical line showing the median. The black whiskers show the range of observed values that fall within 1.5 times the interquartile range and the black circles are outliers above or below it. Stars indicate the total number of snow avalanche paths monitored in each massif.</figcaption></figure>

The conditional dGP model aims at refining the estimation of the probability of having *k* avalanches for a remarkable day *t* given the snow and meteorological conditions of this day *t*. For example, the probability of having very large numbers of avalanches is likely to increase if large amounts of snow have been recorded on this day *t* and/or just before. Note that the terminology “conditional” indicates here the indirect relationship with the snow and meteorological variables and does not refer to a conditional probability distribution (i.e. Pr{X|Y}). This type of model is often referred to as a “non-stationary” model in hydro-meteorological applications (see, e.g. Coles, 2001, Section 6), where the parameters of the distribution vary as functions of time or climate indices (El Adlouni et al., 2007). This model (1) is also a special case of vector generalized linear and additive models (VGLMs and VGAMs) (Yee and Stephenson, 2007).

### 3.3.2. Parameter estimation

In order to estimate the shape parameter $\gamma$ and the scale parameter $\sigma_1$ (for the unconditional dGP model) or $(\sigma_1, \dots, \sigma_{J+1})$ (for the conditional dGP model), we apply the generalized maximum likelihood

estimation (GMLE) method (Martins and Stedinger, 2000; El Adlouni et al., 2007). The shape parameter is difficult to estimate for small samples and setting a pertinent prior information on it helps to reduce this issue (Martins and Stedinger, 2000). In our case, a prior distribution is used to restrict $\gamma$ values to a statistically/physically reasonable range, namely a beta distribution supported on the interval (-0.5, +0.5) and having the probability density function

$$ \pi(\gamma) = \frac{(0.5 + \gamma)^{p-1}(0.5 - \gamma)^{q-1}}{B(p, q)}, $$

where $p = 9$, $q = 6$ and $B(p, q) = \Gamma(p)\Gamma(q)/\Gamma(p+q)$ in which $\Gamma(\cdot)$ is the gamma function. The generalized likelihood function is thus computed as the product between the classical likelihood function and that prior distribution. The generalized maximum likelihood estimator of the parameters $\hat{\gamma}$ and $\hat{\sigma}$ or $(\hat{\sigma}_1, \hat{\sigma}_2, \hat{\sigma}_3, \hat{\sigma}_4, \hat{\sigma}_5)$ can be obtained by maximizing the obtained generalized log-likelihood function

$$ \ell(\sigma(t), \gamma) = \log \pi(\gamma) + \sum_{t=1}^{N} \log \Pr \{E_t = e_t; \sigma(t), \gamma\}. $$
<figure><img src="image_7.png" /><figcaption>Fig. 3. Fitted unconditional dGP distribution versus observed frequencies of remarkable avalanche numbers for the Thabor massif (01-08-1958–31-07-2018).</figcaption></figure>

Here, estimates of the parameters were computed numerically, using the function maxLik() from the maxLik package (Henningsen and Toomet, 2011) implemented in the R software (R, 2017) using the robust Simulated Annealing algorithm (SANN) (Bélisle, 1992) as the main numerical optimization method. The standard errors (SE) of these parameters can be obtained in the usual way from standard likelihood theory (Coles, 2001). These standard errors can be used to assess the significance of the parameter estimates with the statistical T-test implemented in the maxLik package.

### 3.3.3. Model evaluation

To test the goodness-of-fit (GOF) of the sample exceedances with the dGP model, we apply the Kolmogorov–Smirnov GOF test using the ks.testk() function of the dgof package in R, which provides an implementation adapted for discrete distributions (see Arnold and Emerson, 2011, for further details). Here, we consider two-sided tests, and the null hypothesis $H_0$ “the observed remarkable avalanche occurrence numbers follow a dGP distribution” can be rejected with the level of significance $\alpha = 0.05$. Exact corresponding p-values can be obtained, and the null hypothesis is rejected if the p-value < 0.05.

### 3.3.4. Model selection

We evaluate the predictive performance of the dGP models by comparing the fitted probabilistic distribution to observed remarkable avalanche cycles $e_t$, using the log-score (Gneiting and Raftery, 2007), defined as:

$$ \mathrm{LOG\_SCORE}(G, y) = - \log g(y), $$

where y is the observation, G is the predictive cumulative distribution function and g is its corresponding probability density function. In other words, LOG SCORE(G, y) is the negative logarithm of the predictive density, g, evaluated at the observation y. This score is strictly proper and negatively oriented.

Because of the discrete nature of the dGP distribution used in this study, we translate the log-score into:

$$ \mathrm{LOG\_SCORE}(G_{E_t}, e_t) = - \log \mathrm{Pr} \{ E_t = e_t; \sigma(t), \gamma \}, $$

for $t = 1, ..., N$, where $G_{E_t}$ is the predictive cumulative distribution function for the observation $e_t$ which can be either an unconditional or a conditional dGP model.

The log-scores are computed using a leave-one-out cross-validation algorithm. For all observed $e_t$, $t = 1, ..., N$, parameter estimation of the dGP model is performed without $e_t$ (including related covariates $F_j(t)$

for the conditional dGP model). The resulting fitted dGP distribution is denoted by $\hat{G}_{-t}$. LOG SCORE($\hat{G}_{-t}$, $e_t$) is computed for this observed $e_t$. The final score for this dGP model is

$$ \overline{\mathrm{LOG\_SCORE}} = \frac{1}{N} \sum_{t=1}^{N} \mathrm{LOG\_SCORE}(\hat{G}_{-t}, e_t), \qquad (2) $$

i.e., the average value of the log-scores computed at each observed $e_t$.

For each massif, we select the dGP model (conditional or not) generating the smallest average log-score $\overline{\mathrm{LOG\_SCORE}}$.

### 3.3.5. Return level estimation

The T-year ($T > 2$) return level $e(T, t)$ (the number of daily avalanche occurrences which is expected to be exceeded in average once every T years, and which considers the meteorological and snow conditions of the day t for the conditional dGP model) is defined by

$$ e(T, t) = [h] + 1 + \left\lceil \frac{\sigma(t)}{\gamma} \left[ (T n_y)^{\gamma} - 1 \right] \right\rceil, $$

where $n_y = N/60$ is the average number of remarkable avalanche cycles per year. Our data span a period of 60 years from 1958 to 2018. Hence, in all massifs, $n_y \approx 0.5$ since we consider daily occurrence numbers exceeded once every 2 years in average.

## 4. Results

In this Section we present the application of unconditional and conditional dGP models to $e_t$, $t = 1, ..., N$ in the 23 French Alps massifs. For each massif, four different conditional dGP models are tested, using a total number of $J = 1, ..., 4$ transformed covariates $F_j(t)$ (i.e. $J = 1, ..., 4$ first principal components obtained using a PCA analysis).

Fig. 3 provides an illustration of a fitted unconditional dGP distribution for the massif Thabor. Obviously, we see that the fitted dGP distribution provides smoother decaying frequencies than observed frequencies, which are necessarily restricted to few numbers (related to the possible round numbers of avalanches), but overall, model fit to data appears as fair as confirmed by the Kolmogorov–Smirnov GOF test ($p$-value=0.97).

Fig. 4 shows the unconditional and conditional (with three covariates, see selection in Table 2 below) dGP distributions for the Haute-Maurienne massif. The conditional one corresponds to the snow and meteorological conditions for the 4 January 2018. Due to the difference of $\sigma(t)$ values between the two models, the conditional dGP distribution of that day ($\sigma(t) = 44.18$) is much more spread than the
<figure><img src="image_8.png" /><figcaption>Fig. 4. (a) Unconditional (in red) and conditional (in cyan) dGP distributions for the Haute-Maurienne massif for the 4 January 2018. We also provide the predictive intervals (PI) associated to intervals of probabilities [0.25, 0.75] and [0.17, 0.83] (i.e. 50% and 66% predictive intervals respectively) for both models. The dotted vertical line indicates that 60 avalanches were observed that day. (b–g) Number of avalanches and snow and meteorological conditions for the period between the 23 December 2017 and the 12 January 2018 in the Haute-Maurienne massif. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

unconditional version ($\sigma(t) = 9.02$). As a consequence, the probability for a remarkable avalanche cycle to exceed 20 avalanches in one day is 0.58 for the unconditional version, and decreases to 0.04 for the probability of exceeding 50 avalanches. For the conditional version, these probabilities increase to 0.89 and 0.46, respectively, which illustrates that given the meteorological conditions that prevailed at that time, a very high avalanche activity is much more likely. We remind here that these predictive distributions are obtained for the subset of days where a remarkable avalanche number has occurred, which are rare events that happen with a small probability (i.e. the probability of exceeding the threshold $h$).

The 50% and 66% predictive intervals provided for both models also highlight an important discrepancy between the two distributions. The very large number of avalanches recorded that date is related to the storm “Eleanor”, a major European storm that affected western Europe on the 2–3 January 2018. As shown in Fig. 4c–g, this storm caused very high precipitation amounts (snow and rain), and strong winds. Furthermore, mild temperatures on the 2–3 January 2018 led to a rise of the rain/snow limit above 2000 m, and to the formation of a thick layer of wet snow favorable to avalanche triggering. Note also that this event was shortly followed by another avalanche cycle on the 8 January 2018, which was due to a very large amount of snow (see Fig. 4c).

### 4.1. Model performance and parameter estimates

Table 2 provides the average values of the log-score of each model. As indicated above, since the log-score is negatively oriented, the smallest values indicate the best predictive performances. In 12 out of the 23 massifs, the unconditional dGP model is preferred to the conditional alternative versions. For the 11 remaining massifs, the conditional dGP models exhibit the best performances. Among these 11 massifs, five conditional dGP models employ one transformed covariate $F_j(t)$, one model employs two transformed covariates, four models employ three transformed covariates, and one incorporates four transformed covariates. The inclusion of meteorological and snow covariates thus improves the predictive performances of $E_i$ for about half of the massifs. As discussed in Section 5, possible avenues could be considered in order to increase the predictive performance of the conditional dGP models.

Table 3 provides, in each massif, the parameter estimates $\hat{\sigma}_1, \hat{\sigma}_2, \hat{\sigma}_3, \hat{\sigma}_4, \hat{\sigma}_5$ and $\hat{\gamma}$ for the best dGP model, i.e. the unconditional or conditional dGP model exhibiting the best predictive performance according to the log-score LOG SCORE. As indicated above, the T-test is used to assess if the parameter estimates differ from 0 with a significance level of 0.05. The parameter estimates $\hat{\sigma}_1$ are always significant which is not
<caption>**Table 2**</caption>

Average log-score LOG SCORE for each tested dGP model, i.e., in each massif, the unconditional dGP model (Unc.) and conditional dGP models Cond-1–Cond-4, corresponding to 1 to 4 transformed covariates $F_t(t)$ for the scale parameter of the dGP distribution. Bold numbers highlight the best model in each massif.

<table><caption>**Table 3**</caption><thead><tr><th>Massif</th><th>Unc.</th><th>Cond-1</th><th>Cond-2</th><th>Cond-3</th><th>Cond-4</th></tr></thead><tbody><tr><td>Chablais</td><td>2.85</td><td>2.89</td><td><strong>2.79</strong></td><td>3.06</td><td>3.08</td></tr><tr><td>Aravis</td><td><strong>2.92</strong></td><td>2.95</td><td>2.96</td><td>3.00</td><td>3.01</td></tr><tr><td>Mont-Blanc</td><td>3.19</td><td><strong>3.18</strong></td><td>3.27</td><td>3.22</td><td>3.40</td></tr><tr><td>Bauges</td><td><strong>2.25</strong></td><td>2.30</td><td>2.30</td><td>2.38</td><td>2.46</td></tr><tr><td>Beaufortain</td><td><strong>2.32</strong></td><td>2.35</td><td>2.37</td><td>2.46</td><td>3.41</td></tr><tr><td>Haute-Tarentaise</td><td><strong>3.02</strong></td><td>3.03</td><td>3.06</td><td>3.11</td><td>3.23</td></tr><tr><td>Chartreuse</td><td><strong>2.44</strong></td><td>2.46</td><td>2.51</td><td>2.53</td><td>2.61</td></tr><tr><td>Belledonne</td><td><strong>2.25</strong></td><td>2.31</td><td>2.27</td><td>2.31</td><td>2.32</td></tr><tr><td>Maurienne</td><td>3.23</td><td>3.45</td><td>3.13</td><td><strong>3.04</strong></td><td>3.12</td></tr><tr><td>Vanoise</td><td>3.03</td><td><strong>3.02</strong></td><td>3.03</td><td>3.05</td><td>3.12</td></tr><tr><td>Haute-Maurienne</td><td>3.37</td><td>3.45</td><td>3.60</td><td><strong>3.25</strong></td><td>3.31</td></tr><tr><td>Grandes Rousses</td><td>2.57</td><td>2.54</td><td>2.56</td><td><strong>2.54</strong></td><td>2.54</td></tr><tr><td>Thabor</td><td><strong>2.53</strong></td><td>2.60</td><td>2.63</td><td>2.62</td><td>2.67</td></tr><tr><td>Vercors</td><td>1.90</td><td><strong>1.64</strong></td><td>1.65</td><td>1.67</td><td>1.66</td></tr><tr><td>Oisans</td><td><strong>3.60</strong></td><td>3.63</td><td>3.65</td><td>3.62</td><td>3.68</td></tr><tr><td>Pelvoux</td><td>2.67</td><td><strong>2.62</strong></td><td>2.65</td><td>2.80</td><td>3.01</td></tr><tr><td>Queyras</td><td>2.84</td><td><strong>2.72</strong></td><td>2.84</td><td>2.83</td><td>6.51</td></tr><tr><td>Devoluy</td><td><strong>1.84</strong></td><td>1.99</td><td>1.98</td><td>2.02</td><td>2.09</td></tr><tr><td>Champsaur</td><td>2.80</td><td>2.85</td><td>2.88</td><td>2.91</td><td><strong>2.75</strong></td></tr><tr><td>Embrunais Parpaillon</td><td><strong>1.77</strong></td><td>1.86</td><td>1.95</td><td>1.96</td><td>2.19</td></tr><tr><td>Ubaye</td><td><strong>2.78</strong></td><td>2.80</td><td>2.85</td><td>2.86</td><td>3.01</td></tr><tr><td>Haut-Var Haut-Verdon</td><td><strong>1.84</strong></td><td>1.97</td><td>1.98</td><td>1.94</td><td>2.03</td></tr><tr><td>Mercantour</td><td>2.27</td><td>2.30</td><td>2.21</td><td><strong>2.16</strong></td><td>2.19</td></tr></tbody></table>

Parameter estimates $\hat{\sigma}_1, \hat{\sigma}_2, \hat{\sigma}_3, \hat{\sigma}_4, \hat{\sigma}_5$ and $\hat{\gamma}$ for the best fitted dGP distribution fitted to remarkable daily avalanche occurrence numbers in each of the 23 massifs of the French Alps. Values between brackets indicate the corresponding standard errors (SE). N and $[h]$ indicate the number of these remarkable avalanche cycles over the 60 year of processed data, and the threshold applied to select them, respectively. Bold numbers highlight significantly non-zero estimates, according to the result of the T-test ($p$-value lower than 0.05).

<table><thead><tr><th></th><th>N</th><th>[h]</th><th>σ̂<sub>1</sub></th><th>σ̂<sub>2</sub></th><th>σ̂<sub>3</sub></th><th>σ̂<sub>4</sub></th><th>σ̂<sub>5</sub></th><th>γ̂</th></tr></thead><tbody><tr><td>Chablais</td><td>26</td><td>15</td><td>1.67 (0.22)</td><td>0.03 (0.10)</td><td>-0.26 (0.11)</td><td></td><td></td><td>0.03 (0.14)</td></tr><tr><td>Aravis</td><td>26</td><td>13</td><td>6.10 (1.64)</td><td></td><td></td><td></td><td></td><td>0.07 (0.12)</td></tr><tr><td>Mont-Blanc</td><td>29</td><td>12</td><td>2.12 (0.21)</td><td>-0.08 (0.07)</td><td></td><td></td><td></td><td>0.01 (0.14)</td></tr><tr><td>Bauges</td><td>28</td><td>10</td><td>3.30 (0.70)</td><td></td><td></td><td></td><td></td><td>0.03 (0.13)</td></tr><tr><td>Beaufortain</td><td>25</td><td>6</td><td>3.36 (0.78)</td><td></td><td></td><td></td><td></td><td>0.07 (0.13)</td></tr><tr><td>Haute-Tarentaise</td><td>27</td><td>13</td><td>6.99 (1.30)</td><td></td><td></td><td></td><td></td><td>0.04 (0.12)</td></tr><tr><td>Chartreuse</td><td>16</td><td>4</td><td>3.42 (1.01)</td><td></td><td></td><td></td><td></td><td>0.11 (0.12)</td></tr><tr><td>Belledonne</td><td>24</td><td>8</td><td>2.86 (0.75)</td><td></td><td></td><td></td><td></td><td>0.13 (0.12)</td></tr><tr><td>Maurienne</td><td>27</td><td>14</td><td>1.91 (0.22)</td><td>-0.10 (0.07)</td><td>0.34 (0.10)</td><td>0.21 (0.14)</td><td></td><td>0.03 (0.14)</td></tr><tr><td>Vanoise</td><td>28</td><td>15</td><td>1.92 (0.21)</td><td>0.09 (0.07)</td><td></td><td></td><td></td><td>0.03 (0.13)</td></tr><tr><td>Haute-Maurienne</td><td>27</td><td>15</td><td>2.12 (0.22)</td><td>0.09 (0.08)</td><td>0.02 (0.08)</td><td>-0.31 (0.11)</td><td></td><td>0.07 (0.14)</td></tr><tr><td>Grandes Rousses</td><td>29</td><td>10</td><td>1.37 (0.22)</td><td>0.17 (0.09)</td><td>0.13 (0.11)</td><td>0.13 (0.10)</td><td></td><td>0.03 (0.13)</td></tr><tr><td>Thabor</td><td>28</td><td>7</td><td>4.25 (0.91)</td><td></td><td></td><td></td><td></td><td>0.05 (0.13)</td></tr><tr><td>Vercors</td><td>25</td><td>4</td><td>0.57 (0.23)</td><td>-0.29 (0.08)</td><td></td><td></td><td></td><td>-0.00 (0.14)</td></tr><tr><td>Oisans</td><td>30</td><td>23</td><td>11.77 (2.11)</td><td></td><td></td><td></td><td></td><td>0.10 (0.12)</td></tr><tr><td>Pelvoux</td><td>30</td><td>11</td><td>1.57 (0.20)</td><td>0.12 (0.06)</td><td></td><td></td><td></td><td>-0.01 (0.13)</td></tr><tr><td>Queyras</td><td>27</td><td>9</td><td>1.58 (0.23)</td><td>-0.19 (0.07)</td><td></td><td></td><td></td><td>0.08 (0.13)</td></tr><tr><td>Devoluy</td><td>25</td><td>6</td><td>1.99 (0.49)</td><td></td><td></td><td></td><td></td><td>0.10 (0.12)</td></tr><tr><td>Champsaur</td><td>28</td><td>14</td><td>1.55 (0.22)</td><td>0.03 (0.08)</td><td>-0.17 (0.09)</td><td>-0.09 (0.11)</td><td>0.38 (0.13)</td><td>0.07 (0.12)</td></tr><tr><td>Embrunais Parpaillon</td><td>26</td><td>4</td><td>1.80 (0.45)</td><td></td><td></td><td></td><td></td><td>0.11 (0.12)</td></tr><tr><td>Ubaye</td><td>28</td><td>12</td><td>5.20 (1.09)</td><td></td><td></td><td></td><td></td><td>0.08 (0.12)</td></tr><tr><td>Haut-Var Haut-Verdon</td><td>17</td><td>6</td><td>1.93 (0.54)</td><td></td><td></td><td></td><td></td><td>0.09 (0.12)</td></tr><tr><td>Mercantour</td><td>26</td><td>14</td><td>0.99 (0.22)</td><td>-0.05 (0.07)</td><td>0.26 (0.10)</td><td>0.19 (0.11)</td><td></td><td>0.02 (0.14)</td></tr></tbody></table>

surprising considering that they correspond to the constant term of the scale parameter. For the selected conditional models, one $\hat{\sigma}_j, j > 1$ estimate is often significant, indicating that the importance of some covariates are clearly identified. In 21 cases out of the 23 massifs, the shape parameter $\hat{\gamma}$ is greater than zero. This indicates that the distribution of $E_t$ in most massifs of the French Alps is unbounded. However, when $\hat{\gamma}$ is positive, its value is often close to zero (i.e. $\hat{\gamma} < 0.1$) and non-significantly different from zero, indicating a light tail property for the probability distribution of extreme avalanche cycles in these French Alps massifs. However, it must be noticed that the uncertainties for $\hat{\gamma}$ are rather important (standard errors between 0.12 and 0.14), which limit their interpretation.

## 4.2. Return levels of extreme avalanche cycles from unconditional dGP models

Fig. 5 shows the return level curves derived from the unconditional dGP models, for all the massifs of the French Alps, grouped in two sub-regions according to their geographic location. Return level curves clearly highlight the differences between the magnitudes of avalanche cycles for the different massifs of the French Alps. Massifs having lower elevations (e.g. Vercors) show smaller return levels, the corresponding daily avalanche occurrences being less than 20 even for a 300-year return period. By contrast, Oisans and Haute-Maurienne massifs clearly show the highest avalanche numbers, with a 100-year return level of 80 and 59 avalanches in one day, respectively. Maps in Fig. 6 sum-up the centennial and tri-centennial return levels in terms of daily number of avalanches and daily number of avalanche per path. They
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Return level plots of extreme avalanche numbers (Daily number of avalanches versus return period) for the 23 massifs of the French Alps evaluated using the fitted unconditional dGP models.</figcaption></figure>

confirm, e.g., that daily extreme occurrence numbers are higher in Oisans than in all other massifs of the French Alps, whereas weakest extreme avalanche cycles are located in Embrunais-Parpaillon. Maps of the return levels of the number of avalanches by path in Fig. 6b–d shows a slightly different picture, some massifs (e.g. Vercors) being possibly among the most active when the return levels are related to the number of monitored paths.

Table 4 further provides a few statistics (mean, standard deviation and coefficient of variation) summarizing the inter-massif variability of these return levels. While the coefficient of variations of the daily number of avalanches is 0.49 and 0.50 for return periods of 100 and 300 years, respectively, it decreases to 0.30 and 0.31 for the daily number of avalanches per path corresponding to the same return periods. Hence, interestingly, inter-massif dispersion is much higher for the daily number of avalanches than for the daily number of avalanche per path corresponding to a given return period. This illustrates the obvious scaling effect of the total number of paths by massif, which can be used to summarize extreme avalanche activity in the French Alps in a meaningful and potentially useful way. For example, 0.25 avalanches per path (roughly one avalanche for four paths) and 0.32 avalanches per path (≈ one avalanche for three paths) seem rather robust approximations of the activity to be expected during extreme avalanche cycles corresponding to return periods of 100 and 300 years, respectively, all over the French Alps.

<caption>**Table 4**</caption>

Empirical inter massif mean, standard deviation (SD) and coefficient of variation (CV) of the daily number of avalanches (raw and per path) corresponding to return periods of 100 and 300 years.

<table><thead><tr><th rowspan="2">T (years)</th><th colspan="2">Nb of aval.</th><th colspan="2">Nb of aval. / path</th></tr><tr><th>100</th><th>300</th><th>100</th><th>300</th></tr></thead><tbody><tr><td>Mean</td><td>33.7</td><td>42.0</td><td>0.25</td><td>0.32</td></tr><tr><td>SD</td><td>16.5</td><td>20.9</td><td>0.08</td><td>0.10</td></tr><tr><td>CV</td><td>0.49</td><td>0.50</td><td>0.30</td><td>0.31</td></tr></tbody></table>

### 4.3. Predictive intervals of extreme avalanche cycles from conditional dGP models

Fig. 7 compares the predictive intervals of remarkable avalanche cycles obtained with the unconditional and with the best conditional dGP model, for 4 massifs where the conditional dGP model is preferred. While the predictive intervals obtained with unconditional and conditional dGP models are often similar, differences are huge for some days. For these specific days, the scale parameter obtained with the fitted conditional dGP model (see Eq. (1)) is really much larger than the estimated scale parameter of the unconditional dGP model. The reason is that unusual meteorological and snow conditions for these days lead to high absolute values for the transformed covariates $F_j(t)$. As a result, a stretched conditional dGP distribution is obtained, leading
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Daily number of avalanches (a–c) and Daily number of avalanches per path (b–d) corresponding to return periods of 100 and 300 years for the 23 massifs of the French Alps. In each massif, return levels are evaluated using the fitted unconditional dGP models.</figcaption></figure>

to very large return levels (see illustration in Fig. 4) generally well corroborated by observations (very large number of avalanches occurrences). As an example, on 2018-01-04 in the massif Haute-Maurienne, using the selected conditional dGP model, the observed number of avalanches corresponds roughly to the median of the distribution and the predictive intervals cover well this very high number of avalanches. For the unconditional counterpart, even the 95% predictive interval does not reach the observed value. In other words, given the snow and meteorological corresponding to this specific day, the occurrence of a large number of avalanches becomes more likely than using the unconditional model. By contrast, Fig. 7 shows that, for a few snow avalanche episodes, (i) the predictive intervals are actually smaller for the conditional dGP model than for the unconditional dGP model, and (ii) the observed value exceeds the 95% predictive interval of the conditional dGP model. One of these examples is 1980-05-05 for the massif Champsaur. Also, in some cases, it seems that the snow and meteorological conditions were favorable to a very high number of avalanches according to the predictive intervals obtained with the conditional dGP (see, e.g. on 1990-02-14 in Grandes Rousses or on 2004-01-13 in Maurienne), but that this was not corroborated by the observation of a high avalanche activity. However, all in all, Fig. 7 shows that most of the time predictive intervals provided by the conditional dGP model leads to a better coverage of actual avalanche observations than with the unconditional version. This overall improvement of the predictive performances is quantified by the log-scores provided in Table 2, which lead to the selection of the conditional dGP models for these massifs.

Fig. 8 shows the shape parameters $\gamma$ of the unconditional and conditional dGP distributions for the massifs where a conditional model has been selected. For all these massifs, the shape parameters have decreased, indicating that the scale function $\sigma(t)$ of the conditional model, using covariates, compensates the need of having a high shape parameter in order to fit the tail of the distributions. Since the number of possible avalanches is physically limited, these distributions should, in theory, be bounded. The shape parameter shifting systematically towards bounded dGP distributions (i.e. decreasing) with the conditional model provides an additional evidence that incorporating snow and meteorological covariates improves the probabilist representation of these avalanche cycles.

# 5. Discussion, conclusion and outlooks

## 5.1. Methodology

High-magnitude events are by definition rare, which pleads for using robust methods to extrapolate beyond observational records. In this study, we introduce a rigorous statistical framework based on extreme value theory which aims at representing the distribution of avalanche cycles, defined by daily occurrence numbers at the massif scale exceeding the 2-year return level. Due to the discrete nature of the data, we apply a discrete version of the generalized Pareto (dGP) distribution, this study being the first application of the dGP distribution to avalanche data. Two dGP models are proposed, the first one being
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** 50%, 66%, 90% and 95% predictive intervals of remarkable daily number of avalanches corresponding to intervals of probabilities [0.25, 0.75], [0.17, 0.83], [0.05, 0.95] and [0.275, 0.975], respectively, from the fitted unconditional (left) and conditional (right) dGP models, for 4 massifs of the French Alps. For each massif, dates with remarkable avalanche cycles are ordered according to the number of observed avalanches. For these 4 massifs, the conditional dGP model is the best predictive model according to the log-score.</figcaption></figure>

a direct application of the dGP distribution to avalanche cycles, and the second being a conditional dGP model where the scale parameter depends on meteorological and snow covariates. The methodology presented in this study also includes a comprehensive framework for model fitting, selection and evaluation, based on cross-validation and the log-score.

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Shape parameters $\gamma$ of the dGP distribution for the selected conditional models versus their unconditional counterparts. The red box highlights cases for which the unconditional dGP distribution is unbounded and the conditional version is bounded. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 5.2. Main outcomes of the unconditional dGP model

For the 23 massifs of the French Alps, the number of avalanches corresponding to high return periods (e.g. 100 years, 300 years) can be directly obtained from the unconditional dGP model fitted to 60 years of data. These estimates are a first important outcome for public safety. In addition, when expressed as mean numbers of avalanches per path, the unconditional dGP model also provide robust estimates, with rather small variations from one massif to another. These values/numbers will be very useful for technical services to anticipate the most critical situations they may likely have to face under climate conditions characterized by the last 60 years.

It must be noticed that this approach goes far beyond the state of the art of the snow avalanche field, where return periods for extreme avalanche cycles where so far based on empirical distributions (Eckert et al., 2010b) or with continuous approximations (Dkengne Sielenou et al., 2016). Furthermore, this approach could be easily put to work for analyzing extreme counts of other potentially damageable phenomena.

## 5.3. Benefits of the conditional dGP model

The second dGP model evaluated in this study incorporates snow and meteorological covariates in order to improve the probabilistic representation of avalanche cycles. For each day with a very large avalanche activity, this conditional dGP model is more or less stretched according to the snow and meteorological conditions of this day. The conditional dGP model leads to improved predictive performances for about half of the massifs. Clearly, as illustrated in different massifs, the observed number of avalanches is considered as improbable according to the unconditional dGP model. By bringing a physical relationship between meteorological and snow factors and the number of avalanches, the conditional dGP model is able to adapt its statistical representation of avalanche cycles according to the conditions of each day.

By contrast, adding snow and meteorological information did not improve the predictive capacities of dGP models in the other massifs. This recalls that the link between snow conditions and avalanche activity is far from trivial and that the development of more elaborated avalanche activity indices and of statistical tools that can represent
<figure><img src="image_13.png" /><figcaption>Fig. A.9. Correlation matrix of S2M covariates for the Haute-Maurienne massif. (a) All covariates. (b) Retained covariates with variable selection based on a 0.95 absolute correlation threshold. See Table 1 for the description of each covariate, with available elevations (e.g. 1800 m) and aspects (e.g. 180°).</figcaption></figure>

complex relationships between avalanche activity and their drivers is still required. In particular, in this study, as a first attempt, the scale parameter only varies from one day to another according to the values of the covariates.

## 5.4. Potential improvements of conditional models

In the future, different alternative versions of a conditional dGP model could be tested:

* Here, transformed covariates were obtained as a linear transformation of the original covariates. A refined study, for each massif, could be conducted in order to identify the most relevant

meteorological drivers and snow variables for these remarkable avalanche numbers. In addition, many different covariates could be integrated from the outputs of the detailed Crocus snowpack model such as the mechanical stability criteria (Morin et al., 2020);

* In this study, for the conditional dGP model, the logarithm of the scale parameter was expressed as a linear sum of the transformed covariates. Different non-linear combinations could be considered, either using alternative mathematical functions, or more advanced regression models (e.g. random forests);

* The 23 massifs were considered as independent entities. An appealing – yet challenging – avenue would be to account for potential dependencies in avalanche activity and its drivers among
<figure><img src="image_14.png" /><figcaption>**Fig. A.10.** Principal Component Analysis of the retained covariates for the Haute-Maurienne massif. Values indicate the unscaled projections of all covariates, while blue and red colors highlight positive and negative values for each axis, respectively. See Table 1 for the description of each covariate, with available elevations (e.g. 1800 m) and aspects (e.g. 180°).</figcaption></figure>

massifs, which may allow refining the evaluation of the highest return levels by information transfer.

## 5.5. Outlooks for avalanche forecasting and climate change impact studies

A last important practical outcome of this work is the potential application of the conditional dGP model in the context of avalanche forecasting. Combined to efficient statistical avalanche forecasting models which can be used to identify the days during which extreme avalanche activity occurs (Dreier et al., 2016; Choubin et al., 2019; Dkengne Sielenou et al., 2021), and given the daily snow and meteorological conditions, this framework would provide the first-ever probabilistic assessment of avalanche activity likely to occur during the most extreme avalanche cycles. Such conditional probabilistic forecasts are now rather routinely provided and assessed, e.g., in the flood forecasting community (Hamill and Scheuerer, 2018), but have not been developed so far in the snow avalanche community with the exception of few studies and operational implementations (see, e.g. Vernay et al., 2015; Morin et al., 2020). In addition, in the latter, the probabilistic prediction is linked to ensemble snow and meteorological forecasts (i.e. different snow and meteorological conditions provided by a probabilistic forecasting system). Our approach introduces for the first time a different kind of uncertainty in the probabilistic prediction, the one related to the link between avalanche occurrences summed-up at a massif scale and corresponding to prevailing snow and meteorological conditions. Hence, our approach may well deliver new probabilistic insights with high potential to evaluate avalanche risk level during the most severe avalanche cycles as expected in an avalanche bulletin. To this aim, a limitation of the present work is that our snow and meteorological data resulted from reanalyses, by definition not available for real-time forecasting. The impact of forecasted snow and meteorological information on the method accuracy should now be tested, and ultimately combined with ensemble forecast techniques such as in Vernay et al. (2015). This may produce probabilistic avalanche activity forecasts for extreme avalanche cycles at various spatio-temporal scales and horizons that would integrate:

* the forecast uncertainty,

* the epistemic uncertainty related to our imperfect knowledge of the massif scale drivers of snow avalanche activity, as investigated in this paper.

Ultimately, these conditional dGP models could be applied conditionally to available future projections of snow and meteorological conditions (e.g. Verfaillie et al., 2018), which could be used to assess future behaviors of remarkable avalanche numbers under expected future climate conditions. This would compliment the (rare) existing results of that kind that target longer seasonal to annual time scales (Castebrunet et al., 2014).

### CRediT authorship contribution statement

**Guillaume Evin:** Preparation of the data, Development of the final methodology, Codes and finalisation of the manuscript. **Pascal Dkengne Sielenou:** First design of the methodology, Original draft preparation with the first main results. **Nicolas Eckert:** Conceptualization and organisation of the study, Writing - reviewing and editing. **Philippe Naveau:** Reviewing and editing. **Pascal Hagenmuller:** Reviewing and editing. **Samuel Morin:** Reviewing and editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This study was supported by the French Ministry of the Environment, Risk Division (DGPR) through the ECANA project and the “wet snow avalanche” action. INRAE and CNRM/CEN are members of LabEx OSUG.

## Appendix. Selection of covariates: illustration for the Haute Maurienne massif

See Figs. A.9 and A.10.

### References

Ammann, W., Bebi, P., 2000. WSL Institute for Snow and Avalanche Research SLF, Der Lawinenwinter 1999, Ereignisanalyse. Technical Report, SLF Davos, p. 588.

Ancey, C., 2012. Are there “dragon-kings” events (i.e. genuine outliers) among extreme avalanches? Eur. Phys. J. Spec. Top. 205 (1), 117–129. http://dx.doi.org/10.1140/epjst/e2012-01565-7.

Anderson, C.W., 1980. Local limit theorems for the maxima of discrete random variables. Math. Proc. Camb. Phil. Soc. 88 (1), 161–165. http://dx.doi.org/10.1017/S0305004100057443.

Anderson, C.W., Coles, S.G., Hüsler, J., 1997. Maxima of Poisson-like variables and related triangular arrays. Ann. Appl. Probab. 7 (4), 953–971.

Arnold, T.B., Emerson, J.W., 2011. Nonparametric goodness-of-fit tests for discrete null distributions. R J. 3 (2), 34–39.

Asadi, M., Rao, C.R., Shanbhag, D.N., 2001. Some unified characterization results on generalized Pareto distributions. J. Statist. Plann. Inference 93 (1), 29–50. http://dx.doi.org/10.1016/S0378-3758(00)00171-3.

Baggi, S., Schweizer, J., 2009. Characteristics of wet-snow avalanche activity: 20 years of observations from a high alpine valley (Dischma, Switzerland). Nat. Hazards 50 (1), 97–108. http://dx.doi.org/10.1007/s11069-008-9322-7.

Beirlant, J., Goegebeur, Y., Segers, J., Teugels, J.L., 2004. Statistics of Extremes: Theory and Applications. In: Wiley Series in Probability and Statistics, John Wiley & Sons, Google-Books-ID: GtIYLAITcKEC.

Bélisle, C.J.P., 1992. Convergence theorems for a class of simulated annealing algorithms on $R^d$. J. Appl. Probab. 29 (4), 885–895. http://dx.doi.org/10.2307/3214721.

Birkeland, K.W., Mock, C.J., 2001. The major snow avalanche cycle of february 1986 in the Western United States. Nat. Hazards 24 (1), 75–95. http://dx.doi.org/10.1023/A:1011192619039.

Birkeland, K.W., Mock, C.J., Shinker, J.J., 2001. Avalanche extremes and atmospheric circulation patterns. Ann. Glaciol. 32, 135–140. http://dx.doi.org/10.3189/172756401781819030.
Bourova, E., Maldonado, E., Leroy, J., Alouani, R., Eckert, N., Bonnefoy-Demongeot, M., Deschatres, M., 2016. A new web-based system to improve the monitoring of snow avalanche hazard in France. Nat. Hazards Earth Syst. Sci. 16 (5), 1205–1216.
http://dx.doi.org/10.5194/nhess-16-1205-2016.

Castebrunet, H., Eckert, N., Giraud, G., Durand, Y., Morin, S., 2014. Projected changes of snow conditions and avalanche activity in a warming climate: the French Alps over the 2020–2050 and 2070–2100 periods. Cryosphere 8 (5), 1673–1697.
http://dx.doi.org/10.5194/tc-8-1673-2014.

Chakraborty, S., Chakravarty, D., 2014. A discrete gumbel distribution. arXiv:1410.7568 [Math, Stat].

Choubin, B., Borji, M., Mosavi, A., Sajedi-Hosseini, F., Singh, V.P., Shamshirband, S., 2019. Snow avalanche hazard prediction using machine learning methods. J. Hydrol. 577, 123929. http://dx.doi.org/10.1016/j.jhydrol.2019.123929.

Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. In: Springer Series in Statistics, Springer-Verlag London.

Davison, A.C., 1984. Modelling excesses over high thresholds, with an application. In: de Oliveira, J.T. (Ed.), Statistical Extremes and Applications. NATO ASI Series, Springer Netherlands, Dordrecht, pp. 461–482. http://dx.doi.org/10.1007/978-94-017-3069-3-34.

Davison, A.C., Smith, R.L., 1990. Models for exceedances over high thresholds. J. R. Stat. Soc. Ser. B Stat. Methodol. 52 (3), 393–425. http://dx.doi.org/10.1111/j. 2517-6161.1990.tb01796.x

De Haan, L., Ferreira, A., 2007. Extreme Value Theory: An Introduction. Springer Series in Operations Research and Financial Engineering, Google-Books-ID: t6tfXnykazEC.

Dkengne Sielenou, P., Eckert, N., Naveau, P., 2016. A limiting distribution for maxima of discrete stationary triangular arrays with an application to risk due to avalanches. Extremes 19 (1), 25–40. http://dx.doi.org/10.1007/s10687-015-0234-0.

Dkengne Sielenou, P., Viallon-Galinier, L., Hagenmuller, P., Naveau, P., Morin, S., Dumont, M., Verfaillie, D., Eckert, N., 2021. Combining random forests and class-balancing to discriminate between three classes of avalanche activity in the french alps. Cold Reg. Sci. Technol. 187, 103276. http://dx.doi.org/10.1016/j.coldregions. 2021.103276.

Dreier, L., Harvey, S., van Herwijnen, A., Mitterer, C., 2016. Relating meteorological parameters to glide-snow avalanche activity. Cold Reg. Sci. Technol. 128, 57–68. http://dx.doi.org/10.1016/j.coldregions.2016.05.003.

Durand, Y., Giraud, G., Brun, E., Mérindol, L., Martin, E., 1999. A computer-based system simulating snowpack structures as a tool for regional avalanche forecasting. J. Glaciol. 45 (151), 469–484. http://dx.doi.org/10.3189/S0022143000001337.

Durand, Y., Giraud, G., Laternser, M., Etchevers, P., Mérindol, L., Lesaffre, B., 2009a. Reanalysis of 47 years of climate in the french alps (1958–2005): Climatology and trends for snow cover. J. Appl. Meteorol. Climatol. 48 (12), 2487–2512. http://dx.doi.org/10.1175/2009JAMC1810.1.

Durand, Y., Laternser, M., Giraud, G., Etchevers, P., Lesaffre, B., Mérindol, L., 2009b. Reanalysis of 44 yr of climate in the french alps (1958–2002): Methodology, model validation, climatology, and trends for air temperature and precipitation. J. Appl. Meteorol. Climatol. 48 (3), 429–449. http://dx.doi.org/10.1175/2008JAMC1808.1.

Eckert, N., Baya, H., Deschatres, M., 2010a. Assessing the response of snow avalanche runout altitudes to climate fluctuations using hierarchical modeling. Application to 61 winters of data in France. J. Clim. 23 (12), 3157–3180. http://dx.doi.org/10. 1175/2010JCLI3312.1

Eckert, N., Coleou, C., Castebrunet, H., Deschatres, M., Giraud, G., Gaume, J., 2010b. 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 Reg. Sci. Technol. 64 (2), 119–136. http://dx.doi.org/10.1016/j. coldregions.2010.08.009.

Eckert, N., Keylock, C.J., Castebrunet, H., Lavigne, A., Naaim, M., 2013. Temporal trends in avalanche activity in the French Alps and subregions: from occurrences and runout altitudes to unsteady return periods. J. Glaciol. 59 (213), 93–114. http://dx.doi.org/10.3189/2013JoG12J091.

Eckert, N., Naaim, M., Parent, E., 2010c. Long-term avalanche hazard assessment with a bayesian depth-averaged propagation model. J. Glaciol. 56 (198), 563–586. http://dx.doi.org/10.3189/002214310793146331.

El Adlouni, S., Ouarda, T.B.M.J., Zhang, X., Roy, R., Bobée, B., 2007. Generalized maximum likelihood estimators for the nonstationary generalized extreme value model. Water Resour. Res. 43 (3), W03410. http://dx.doi.org/10.1029/2005WR004545.

Embrechts, P., Klüppelberg, C., Mikosch, T., 1997. Modelling Extremal Events: For Insurance and Finance. In: Stochastic Modelling and Applied Probability, Springer-Verlag, Berlin Heidelberg.

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

Gassner, M., Brabec, B., 2002. Nearest neighbour models for local and regional avalanche forecasting. Nat. Hazards Earth Syst. Sci. 2 (3/4), 247–253. http://dx. doi.org/10.5194/nhess-2-247-2002.

Gneiting, T., Raftery, A.E., 2007. Strictly proper scoring rules, prediction, and estimation. J. Amer. Statist. Assoc. 102 (477), 359–378. http://dx.doi.org/10.1198/ 016214506000001437.

Hamill, T.M., Scheuerer, M., 2018. Probabilistic precipitation forecast postprocessing using quantile mapping and rank-weighted best-member dressing. Mon. Weather Rev. 146 (12), 4079–4098. http://dx.doi.org/10.1175/MWR-D-18-0147.1.

Hendrikx, J., Owens, I., Carran, W., Carran, A., 2005. Avalanche activity in an extreme maritime climate: The application of classification trees for forecasting. Cold Reg. Sci. Technol. 43 (1), 104–116. http://dx.doi.org/10.1016/j.coldregions.2005.05. 006.

Henningsen, A., Toomet, O., 2011. maxLik: A package for maximum likelihood estimation in R. Comput. Statist. 26 (3), 443–458. http://dx.doi.org/10.1007/s00180- 010-0217-1.

Höller, P., 2009. Avalanche cycles in Austria: an analysis of the major events in the last 50 years. Nat. Hazards 48 (3), 399–424. http://dx.doi.org/10.1007/s11069- 008-9271-1.

International Commission of Snow and Ice, 1981. Avalanche Atlas: Illustrated International Avalanche Classification. UNESCO, Paris, OCLC: 858384427.

Jomelli, V., Delval, C., Grancher, D., Escande, S., Brunstein, D., Hetu, B., Filion, L., Pech, P., 2007. Probabilistic analysis of recent snow avalanche activity and weather in the French Alps. Cold Reg. Sci. Technol. 47 (1–2), 180–192. http://dx.doi.org/ 10.1016/j.coldregions.2006.08.003.

Keylock, C.J., 2005. An alternative form for the statistical distribution of extreme avalanche runout distances. Cold Reg. Sci. Technol. 42 (3), 185–193. http://dx. doi.org/10.1016/j.coldregions.2005.01.004.

Krishna, H., Singh Pundir, P., 2009. Discrete Burr and discrete Pareto distributions. Stat. Methodol. 6 (2), 177–188. http://dx.doi.org/10.1016/j.stamet.2008.07.001.

Laternser, M., Schneebeli, M., 2002. Temporal trend and spatial distribution of avalanche activity during the last 50 years in Switzerland. Nat. Hazards 27 (3), 201–230. http://dx.doi.org/10.1023/A:1020327312719.

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

Lavigne, A., Eckert, N., Bel, L., Parent, E., 2015. Adding expert contributions to the spatiotemporal modelling of avalanche activity under different climatic influences. J. R. Stat. Soc. Ser. C. Appl. Stat. 64 (4), 651–671. http://dx.doi.org/10.1111/rssc. 12095.

Leadbetter, M.R., Lindgren, G., Rootzen, H., 1983. Extremes and Related Properties of Random Sequences and Processes. In: Springer Series in Statistics, Springer Verlag, New York.

Martins, E.S., Stedinger, J.R., 2000. Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data. Water Resour. Res. 36 (3), 737–744. http://dx.doi.org/10.1029/1999WR900330.

McClung, D.M., Lied, K., 1987. Statistical and geometrical definition of snow avalanche runout. Cold Reg. Sci. Technol. 13 (2), 107–119. http://dx.doi.org/10.1016/0165- 232X(87)90049-8.

McClung, D., Schaerer, P.A., 2006. The Avalanche Handbook, The Mountaineers Books.

Morin, S., Horton, S., Techel, F., Bavay, M., Coléou, C., Fierz, C., Gobiet, A., Hagenmuller, P., Lafaysse, M., Li' zar, M.z., Mitterer, C., Monti, F., Müller, K., Olefs, M., Snook, J.S., van Herwijnen, A., Vionnet, V., 2020. Application of physical snowpack models in support of operational avalanche hazard forecasting: A status report on current implementations and prospects for the future. Cold Reg. Sci. Technol. 170, 102910. http://dx.doi.org/10.1016/j.coldregions.2019.102910.

Mougin, P., 1922. Les avalanches en Savoie. Technical Report, Ministère de l'Agriculture, Direction Générale des Eaux et Forêts, Service des Grandes Forces Hydrauliques, Paris.

Nadarajah, S., Mitov, K., 2002. Asymptotics of maxima of discrete random variables. Extremes 5 (3), 287–294. http://dx.doi.org/10.1023/A:1024081112501.

Nakagawa, T., 1978. Discrete extreme distributions. IEEE Trans. Reliab. R-27 (5), 367–368. http://dx.doi.org/10.1109/TR.1978.5220425.

Nakagawa, T., Osaki, S., 1975. The discrete Weibull distribution. In: Conference Name: IEEE Transactions on Reliability. IEEE Trans. Reliab. R-24 (5), 300–301. http://dx.doi.org/10.1109/TR.1975.5214915.

Obled, C., Good, W., 1980. Recent developments of avalanche forecasting by discriminant analysis techniques: A methodological review and some applications to the Parsenn area (Davos, Switzerland). J. Glaciol. 25 (92), 315–346. http://dx.doi.org/10.3189/S0022143000010522.

Prieto, F., Gómez-Déniz, E., Sarabia, J.M., 2014. Modelling road accident blackspots data with the discrete generalized Pareto distribution. Accid. Anal. Prev. 71, 38–49. http://dx.doi.org/10.1016/j.aap.2014.05.005.

R, 2017. R: A language and environment for statistical computing. ISBN: 3-900051-07-0, https://www.r-project.org/.

Schweizer, J., Jamieson, J.B., Schneebeli, M., 2003. Snow avalanche formation. Rev. Geophys. 41 (4), http://dx.doi.org/10.1029/2002RG000123.

Schweizer, J., Mitterer, C., Stoffel, L., 2009. On forecasting large and infrequent snow avalanches. Cold Reg. Sci. Technol. 59 (2), 234–241. http://dx.doi.org/10.1016/j. coldregions.2009.01.006.
Serinaldi, F., Kilsby, C.G., 2014. Rainfall extremes: Toward reconciliation after the battle of distributions. Water Resour. Res. 50 (1), 336–352. http://dx.doi.org/10.1002/2013WR014211.

Verfaillie, D., Lafaysse, M., Déqué, M., Eckert, N., Lejeune, Y., Morin, S., 2018. Multi-component ensembles of future meteorological and natural snow conditions for 1500 m altitude in the Chartreuse mountain range, Northern French Alps. Cryosphere 12 (4), 1249–1271. http://dx.doi.org/10.5194/tc-12-1249-2018.

Vernay, M., Lafaysse, M., Hagenmuller, P., Nheili, P., Verfaillie, D., Morin, S., 2019. The S2M meteorological and snow cover reanalysis in the French mountainous areas (1958 - present). AERIS, https://doi.org/10.25326/37.

Vernay, M., Lafaysse, M., Mérindol, L., Giraud, G., Morin, S., 2015. Ensemble forecasting of snowpack conditions and avalanche hazard. Cold Reg. Sci. Technol. 120, 251–262. http://dx.doi.org/10.1016/j.coldregions.2015.04.010.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., Willemet, J.-M., 2012. The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2. Geosci. Model Dev. 5 (3), 773–791. http://dx.doi.org/10.5194/gmd-5-773-2012.

Yee, T.W., Stephenson, A.G., 2007. Vector generalized linear and additive extreme value models. Extremes 10 (1), 1–19. http://dx.doi.org/10.1007/s10687-007-0032-4.