# Estimating changes in extreme snow load in Europe as a function of global warming levels

G. Evin, E. Le Roux, E. Kamir, S. Morin

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

## To cite this version:

G. Evin, E. Le Roux, E. Kamir, S. Morin. Estimating changes in extreme snow load in Europe as a function of global warming levels. Cold Regions Science and Technology, 2025, 231, pp.104424. 10.1016/j.coldregions.2025.104424. hal-04910810

HAL Id: hal-04910810

https://hal.science/hal-04910810v1

Submitted on 12 Sep 2025

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
# Estimating changes in extreme snow load in Europe as a function of global warming levels

G. Evin<sup>a,*</sup>, E. Le Roux<sup>b</sup>, E. Kamir<sup>c</sup>, S. Morin<sup>c,d</sup>

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

## ARTICLE INFO

### Keywords:

Snow load

Discrete-continuous distribution

Nonstationary model

Global warming level

## ABSTRACT

Most European regions regularly experience large snow loads due to large snow accumulations and/or intense snowfalls. In the future, the magnitude of extreme snow loads is generally expected to decrease due to global warming. However, this decrease depends strongly on the latitude and the elevation. This study aims to provide a flexible statistical framework to estimate changes in extreme snow load (50-year return levels) in Europe as a function of global warming levels, using a multi-model ensemble of snow cover projections. Different sets of nonstationary models are introduced to accommodate the nonzero probability of experiencing years without snow. In particular, parsimonious distributions (exponential, Gamma, Inverse-Gamma) are considered when zero SWE maxima are present. This approach is illustrated by an application to European regions. In a + 3 °C world, extreme snow loads are projected to decrease strongly in all of Europe compared to a warming level of +1 °C, which corresponds roughly to the current climate.

## 1. Introduction

Snow load corresponds to the pressure (in N m⁻² or Pa) applied by the snow cover on underlying surfaces (ground, buildings, infrastructures). Following Eurocodes specifications, extreme snow loads denote return levels with a return period of 50 years (a.k.a. 50-year return levels), see Croce et al. (2021) and references therein. Such extreme snow loads are critical for designing infrastructures (Sadovský et al., 2012; Croce et al., 2018a), as they can lead to roof collapses (O'Rourke, 1997; Strasser, 2008; Ousset et al., 2023). They may also cause large-scale damages to ecosystems (forests) or infrastructures (transportation and electricity networks). While these different impacts are often caused by intense snow storms, they may be the result of snow accumulation at the multi-day scale, possibly aggravated by short and intense snowfalls.

Snow load equals the gravitational acceleration constant g = 9.8ms⁻² times the snow water equivalent (SWE), i.e. the mass of snow per unit surface area (in kg m⁻², also reported in mm w.e.). SWE corresponds to the height of liquid water (in mm) if all the considered snow in a 1 m² area would instantaneously melt. Assessing extreme snow load or extreme SWE values is thus equivalent. Extreme SWE values at the

seasonal scale are of keen interest for water resource management (Marty et al., 2017). Because SWE is a variable that is more commonly observed and dealt with than snow load, this study focuses on SWE values; corresponding snow load information can directly be inferred by applying the factor g to the results expressed in terms of SWE.

In a warming climate, changes in snow load are expected to depend strongly on the latitude and the elevation. For the European Alps, the special report from the Intergovernmental Panel on Climate Change (IPCC) on the Ocean and Cryosphere in a Changing Climate, and in particular its Chapter 2 on High mountain areas (IPCC, 2022) assessed a decrease of the mean winter SWE at the end of the century, at all elevations (see Fig. 2.3 in this IPCC report). Concerning extreme snow loads, Le Roux et al. (2022) shows that the severity of projected decreases depends on the elevation and the emission scenario, or, equivalently, to the warming levels. In low-elevation areas in Europe, Croce et al. (2021) estimated that extreme snow loads are generally projected to decrease for the period 2011–2050 with respect to 1951–1990 under RCP4.5 and RCP8.5 scenarios, except for some cold regions in Ireland and Scandinavian countries.

Studies focusing on projected changes in snow load usually rely on high-resolution (typically, 10 km horizontal resolution) climate

<footnote><sup>a</sup> Univ. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, Grenoble 38000, France</footnote>

<footnote><sup>b</sup> Lab-STICC, IMT Atlantique, Brest, France</footnote>

<footnote><sup>c</sup> Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble 38000, France</footnote>

<footnote><sup>d</sup> CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France</footnote>

<footnote>* Corresponding author.
E-mail address: guillaume.evin@inrae.fr (G. Evin).</footnote>
projections to account for the high spatial variability of snow characteristics (especially in mountainous areas). Generally, projections from general circulation models (GCMs) are dynamically downscaled using regional climate models (RCMs). Corrections of systematic biases are then handled by statistically adjusting RCM outputs, and future projections of snow characteristics are finally obtained by using the adjusted climate projections as inputs for a snowpack model (Verfaillie et al., 2017, 2018). Some studies directly use snow cover information computed within the regional climate models (Croce et al., 2021), which are internally consistent with atmospheric fields from the same model runs but can exhibit biases and discrepancies with respect to observations in some cases (Monteiro et al., 2022; Kotlarski et al., 2023).

Projected changes in extreme SWE can be obtained by applying statistical models to the annual maxima of snow load, where the maxima are taken over years centered on the winter season (i.e. August to July in the Northern Hemisphere). In regions where snowfall is frequent, nonstationary extreme value models can be applied to estimate extreme snow loads (Le Roux et al., 2020). In regions where snowfall is rare (e.g. Mediterranean area or at low elevations in Western European countries), the same models can be inadequate (O'Donnell et al., 2020). Indeed, in these regions, annual SWE maxima are often small or equal to zero, whereas extreme value theory relies on the assumption that maxima are sampled from the upper tail of the distributions (Coles, 2001). If the distribution of daily SWE values is composed of a large probability of having a zero and a heavy-tail distribution for non-zero SWE values, zero or small annual SWE maxima are not representative draws of the extreme tail of the distribution. Possible remedies could be to consider larger blocks to extract the maxima (e.g. 10-year maxima) or apply Peak-Over-Threshold approaches with larger thresholds. However, this would lead to small samples and very uncertain estimates. Alternatively, it is possible to apply standard distributions to these "unextreme" maxima and explicitly treat the probability of having a zero maxima.

This study aims to provide a flexible statistical framework that can be used to estimate return levels of projected SWE maxima in all European regions, including those where snowfall is rare. Different statistical distributions are fitted to annual SWE maxima and a statistical criterion is used to select the best model. Different sets of distributions are considered according to the presence of zero values in annual SWE maxima. The proposed approach is illustrated using data derived from the Mountain Tourism Meteorological and Snow Indicators (MTMSI) dataset (Morin et al., 2021), obtained from statistically adjusted EURO-CORDEX climate projections used as input to the snow cover model Crocus. Results are generated for NUTS-3 (Nomenclature des Unités Territoriales Statistiques) geographical areas spanning all of Europe. The choice of using NUTS-3 regions is motivated by their relevance for providing information to local and national authorities in Europe.

Our contributions are twofold. First, we introduce a flexible statistical framework to estimate return levels of snow loads, including the fitting of several nonstationary models of varying complexity which can accommodate the presence of annual SWE maxima equal to zero, the selection of the best model, a goodness-of-fit test for this best model, and the uncertainty of the return level estimates. When zero SWE maxima are present, we consider parsimonious distributions (exponential, Gamma, Inverse-Gamma). Then, our second contribution is to illustrate this approach at a pan-European scale, for many different climates and snow conditions. Following a commonly used approach (IPCC, 2018; IPCC, 2021) and many national and international commitments to reduce emissions such as the Paris Agreement (UNFCCC, 2015), we provide the return levels for specific global warming levels.

This paper is organized as follows. Section 2 presents the MTMSI dataset and global warming levels for each GCM. Section 3 describes the statistical models used to estimate 50-year return levels of snow loads, which consider different distributions with parameters varying as functions of global warming levels. Section 4 illustrates our approach by presenting detailed results obtained for two NUTS-3 regions and presents the main results at a European scale. In Section 5, we discuss

various methodological aspects and some limitations of this work. Section 6 concludes.

## 2. Data

### 2.1. MTMSI dataset

This study relies on data underlying the MTMSI dataset (Morin et al., 2021), which provides a set of snow-related indicators based on the UERRA reanalysis (Soci et al., 2016) and CMIP5 EURO-CORDEX projections (Kotlarski et al., 2014) at NUTS-3 level, by steps of 100 m elevation. In mountainous areas, the data is provided for several elevation steps, while for non-mountainous areas the data is provided at the mean elevation of the NUTS-3 region. The ADAMONT method (Verfaillie et al., 2017) was used to adjust the EURO-CORDEX GCM/ RCM pairs using the UERRA 5.5 km reanalysis as an observation reference. In the MTMSI dataset, 20 future simulation chains (RCP/GCM/ RCM) are available for the 21st century (9 GCM/RCM pairs for RCP4.5 and RCP8.5, including 2 for RCP2.6). Here, we only consider the 9 simulation chains obtained with the emission scenario RCP8.5 (see Table 1) because return levels are estimated for different warming levels which are not necessarily reached by all climate simulations for low or medium emission scenarios. Finally, annual SWE maxima taken over years centered on the winter season (i.e. from August to July) are considered for the period 1951–2099 for each of the 9 simulation chains, under the historical scenario for the period 1951–2005 and the scenario RCP8.5 for the period 2006–2099. Here we focus on the annual maximum of daily SWE values from this dataset.

In addition, we consider only the average elevation for each of the 1524 NUTS-3 regions, i.e., the elevation available in the MTMSI dataset that is closest to the mean elevation of the NUTS-3 region. Fig. 1 shows the distribution of elevation values considered for the 1524 NUTS-3 regions. In this study, we also illustrate some results for two NUTS-3 regions:

* Loire-Atlantique, France is a NUTS-3 region located at a low mean elevation (0 m) and has a suboceanic climate with quite mild and rainy winters,

* Osttirol, Austria is a NUTS-3 region located in the Alps at a high mean elevation (2000 m) and is typical of an Alpine climate with cold and snowy winters.

### 2.2. Global warming levels for each GCM

For each simulation chain RCP/GCM/RCM, the climate simulation obtained from the corresponding GCM provides the average temperature at the planetary level (Evin et al., 2024). The global mean surface temperatures are averaged annually over the period 1860–2100 (historical runs for 1860–2005, and GCM runs forced by the emission scenario RCP8.5 for 2006–2100). To reduce the effect of internal variability, the mean annual global temperatures are then smoothed

<caption>**Table 1**</caption>

CMIP5 EURO-CORDEX GCM/RCM pairs used in the MTMSI dataset. The part of model names indicated in blue refers to short names used in the rest of the paper.

<table><thead><tr><th>GCM</th><th>RCM</th></tr></thead><tbody><tr><td>MPI-M-MPI-ESM-LR</td><td>CCLM4-8-17</td></tr><tr><td>MPI-M-MPI-ESM-LR</td><td>MPI-CSC-REMO2009</td></tr><tr><td>MPI-M-MPI-ESM-LR</td><td>SMHI-RCA4</td></tr><tr><td>CNRM-CERFACS-CNRM-CM5</td><td>CNRM-ALADIN53</td></tr><tr><td>CNRM-CERFACS-CNRM-CM5</td><td>SMHI-RCA4</td></tr><tr><td>IPSL-CM5A-MR</td><td>IPSL-INERIS-WRF331F</td></tr><tr><td>IPSL-CM5A-MR</td><td>SMHI-RCA4</td></tr><tr><td>ICHEC-EC-EARTH</td><td>SMHI-RCA4</td></tr><tr><td>MOHC-HadGEM2-ES</td><td>SMHI-RCA4</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Mean elevation of the 1524 NUTS-3 regions. The two NUTS-3 regions in cyan (Loire-Atlantique, France, and Osttirol, Austria) are used to illustrate some results. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

using cubic splines with a degree of freedom equal to 6 (Rigal et al., 2019; Ribes et al., 2022; Evin et al., 2024). Anomalies of the global mean temperatures relative to the pre-industrial period 1860–1900 are considered and referred to as global warming levels hereafter.

Fig. 2 shows the global warming levels for the different GCMs. These smoothed global warming levels simulated by each GCM can be compared to observed global warming levels from HadCRUT5 (Morice et al., 2021) which provides a gridded dataset of surface temperatures. For future periods, the warming level reached by the different climate projections depends on the climate sensitivity of each GCM. At the end of the century, GCMs IPSL-CM5A-MR and HadGEM2-ES project global warming levels close to +6 °C compared to the pre-industrial period, while they are slightly lower than +5 °C for the other GCMs.

SWE maxima and global warming levels are available each year for the different simulation chains, as shown in Fig. 2. Fig. 3 shows the SWE maxima obtained in the two illustrative NUTS-3 regions, for the different simulation chains, as a function of the global warming levels. In Loire-Atlantique, France (Fig. 3a), SWE maxima are small at +1 °C

(current situation) and often equal to zero. For higher warming levels (e. g. +3 °C, +4 °C), the proportion of zeros increases but large maxima are still possible. This example shows that all the statistical properties of the SWE maxima are impacted by the warming level (proportion of zeros, mean, dispersion) and must potentially be considered nonstationary. In Osttirol, Austria (Fig. 3b), SWE maxima are larger and far from zero at +1 °C, and decrease as a function of the warming level.

### 3. Nonstationary models

This section presents the statistical framework used to estimate 50-year return levels of SWE for each of the 1524 NUTS-3 regions as a function of the global warming level $w_g(t)$, where g denotes the GCM and t the year. In the remainder of this section, the global warming level is simply denoted by w to lighten the notations. Several statistical models are applied to the annual maxima of daily SWE values, denoted by $X_i(t)$, where i refers to the simulation chain. As discussed in Section 1, a major difficulty is the adequate treatment of SWE maxima for areas with rare snowfalls. In these regions, SWE maxima are often close to or equal to 0. Different nonstationary models are thus considered according to the proportion of zero maxima in a simulation chain, where SWE maxima smaller than 1 mm w.e. are considered zeros.

Table 2 provides the cumulative distribution functions (CDF) of the candidate nonstationary models. These models consider that parameters of the distribution evolve linearly with the global warming level w after some transformations. For each simulation chain, the set of candidate models differs depending on the number of SWE maxima equal to zero:

* **>90 % of zeros:** If the percentage of nonzero SWE maxima is less than 10%, we consider that there is insufficient data to provide any estimates.

* **With zeros:** Otherwise, if at least one zero SWE maximum is present, we consider a discrete-continuous distribution with a probability mass $p(w)$ in zero and a continuous distribution for non-zero values, with parameters $\mu(w)$ and $\sigma(w)$. The continuous distribution can be either exponential, Gamma, or Inverse-Gamma, and the corresponding nonstationary models are referred to as pEXP, pGA and pIG, respectively. Considering that few nonzero values can be available, parsimonious distributions are thus preferred for these cases and we do not apply extreme value distribution.

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Global warming levels from the GCMs and HadCRUT5 in degrees Celsius (one colour per GCM), i.e. anomalies for global annual mean surface temperature compared to the pre-industrial period 1860–1900. For the period 2006–2100, GCMs are forced by the scenario RCP8.5. Raw and smooth global warming levels are shown with dash and plain lines, respectively, and HadCRUT5 values are shown in black for the sake of comparison.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. SWE maxima as a function of the global warming level, for the different simulation chains and the two illustrative NUTS-3 regions. (a) Loire-Atlantique, France (b) Osttirol, Austria.</figcaption></figure>

* Without zeros: Finally, if there are no zero SWE maxima, different continuous distributions are considered: Gamma, Generalized Gamma, Gumbel, or the generalized extreme value distribution (models GA, GG, GUM and GEV, respectively). They are parameterized by μ(w) and σ(w), and an extra shape parameter ξ for the models GG and GEV. These different continuous distributions are generally suited to skewed data, like daily precipitation (Papalexiou et al., 2013). The Gumbel distribution is a special case of the GEV distribution and has been applied to determine the European standards for structural design, i.e. Eurocodes (for more details, see Croce et al., 2019; Le Roux et al., 2020). Other distributions (e.g. Log-normal, Inverse-Gamma, Pareto) were considered during the first attempts but were discarded due to their poor performances.

The parameters p(w) and σ(w) are linear functions of w after logit and log-transformation for p and σ, respectively:

$$ 
\begin{aligned}
\log\left(\frac{p(w)}{1-p(w)}\right) &= \alpha_p + \beta_p \times w, \\
\mathcal{T}(\mu(w)) &= \alpha_{\mu} + \beta_{\mu} \times w, \\
\log(\sigma(w)) &= \alpha_{\sigma} + \beta_{\sigma} \times w.
\end{aligned}
$$

The transformation T for μ(w) is logarithmic or not depending on the distribution (Table 2), and follows the default implementation of the packages gamlss and gamlssx. For the GG and GEV distributions, an additional shape parameter ξ must be estimated. This parameter does not evolve as a function of the global warming level. For the GEV distribution, this choice has been discussed in many contributions (see, e.g. Coles, 2001; Cooley, 2013). This parameter is difficult to estimate, and a changing ξ implies that critical characteristics of the distribution would evolve with the warming level, such as its support (possible minimum and maximum values for the SWE maxima) and the behavior of the tail. The nonstationary models are fitted using maximum likelihood methods implemented in the R packages gamlss (Stasinopoulos and Rigby,
<caption>**Table 2**</caption>

Name and cumulative distribution functions of the nonstationary models considered, where $0 < p < 1$, $\mu(w) > 0$ ($\mu(w) \in \mathbb{R}$ for GUM and GEV), $\sigma(w) > 0$, $\xi \neq 0$, $\gamma(a, x) = \int_{0}^{x} t^{a-1}e^{-t}dt$ is the incomplete gamma function, $\Gamma(a, x) = \int_{x}^{\infty} t^{a-1}e^{-t}dt$ is the complement of the incomplete gamma function, $\Gamma(a) = \int_{0}^{\infty} t^{a-1}e^{-t}dt$ is the ordinary gamma function, $\Phi(x)$ is the standard Gaussian cumulative distribution.

<table><thead><tr><th>Model name</th><th>CDF(x)</th><th>Log-link for F(μ(w))</th></tr></thead><tbody><tr><td>pEXP</td><td>= p(w) if x = 0<br>= p(w) + (1 - p(w))[1 - exp(-x/μ(w))] if x > 0</td><td>✓</td></tr><tr><td>pGA</td><td>= p(w) if x = 0<br>= p(w) + (1 - p(w))<sup>γ(σ(w)<sup>-2</sup>, xμ(w)<sup>-1</sup>σ(w)<sup>-2</sup>)</sup> / Γ(σ(w)<sup>-2</sup>) if x > 0</td><td>✓</td></tr><tr><td>pIG</td><td>= p(w) if x = 0<br>= p(w) + (1 - p(w))φ(x<sup>s</sup>) + exp{2(μσ<sup>2</sup>)<sup>-1</sup>}φ(-x<sup>s</sup>) if x > 0</td><td>✓</td></tr><tr><td>GA</td><td>with x<sup>s</sup> = (σ<sup>2</sup>x)<sup>-1/2</sup>(x/μ - 1)<br>γ(σ(w)<sup>-2</sup>, xμ(w)<sup>-1</sup>σ(w)<sup>-2</sup>) / Γ(σ(w)<sup>-2</sup>)</td><td>✓</td></tr><tr><td>GG</td><td>= γ(ξ<sup>s</sup>, ξ<sup>s</sup>(x/μ)<sup>s</sup>) / Γ(ξ<sup>s</sup>) if ξ > 0<br>= Γ(ξ<sup>s</sup>, ξ<sup>s</sup>(x/μ)<sup>s</sup>) / Γ(ξ<sup>s</sup>) if ξ < 0</td><td>✓</td></tr><tr><td>GUM</td><td>with ξ<sup>s</sup> = 1/(σ<sup>2</sup>ξ<sup>s</sup>)<br>exp{−exp{−[((x − μ(w))/σ(w))]}}</td><td></td></tr><tr><td>GEV</td><td>exp{−[1 + ξ((x − μ(w))/σ(w))]<sup>-1/s</sup>}</td><td></td></tr></tbody></table>

2008) and *gamlssx* (Northrop, 2024) for the model GEV.

### 3.1. Model selection

The best nonstationary model is selected according to the AIC criterion (Burnham and Anderson, 2004) which is shown to be adequate for the selection of nonstationary models (Kim et al., 2017; Xavier et al., 2019):

$$AIC = 2k - 2\log(L), \quad (1)$$

where *k* is the number of free parameters and *L* is the likelihood. For each NUTS-3 region and each simulation *i*, the model with the minimum AIC value is selected as the best one and models with small AIC differences (i.e. less than 1) are considered equally reliable.

### 3.2. Goodness-of-fit test

We check if the normalized quantile residuals of the best nonstationary model are likely to follow a Gaussian distribution using the Shapiro-Wilk normality test (Royston, 1995). Under the null hypothesis, normalized quantile residuals are normally distributed and the *p*-values are uniformly distributed. If the best model is not adequate for some NUTS-3 regions (i.e. the nonstationary model does not fit the SWE maxima reasonably well) and some simulation chains *i*, the residuals will not be normally distributed and it can be expected that the distribution of all the *p*-values will be more weighted towards 0.

### 3.3. Return level estimates

For a given warming level *w*, the 50-year return level of SWE, denoted as $\tilde{RL}_i(w)$, corresponds to the quantile 0.98 of the best model, for a given simulation chain *i*. A single estimate $\tilde{RL}(w)$ can be provided for a NUTS-3 region by taking the average of the return level estimates obtained for the different simulation chains, i.e. the 9 GCM/RCM combinations.

### 3.4. Uncertainty of the return level estimates

Finally, the uncertainty in the return level estimates is evaluated using a nonparametric bootstrap method (Efron and Tibshirani, 1994) with *n* = 100 replicates, as a trade-off between computational time and description of uncertainty.

## 4. Results

### 4.1. Detailed results for two NUTS-3 regions

Fig. 4 illustrates the return levels obtained with the different nonstationary models and simulation chains, for the NUTS-3 Loire-Atlantique where SWE maxima are often equal to zero. For each simulation chain, black points represent the annual SWE maxima as a function of the global warming level reached for the respective years (as shown in Fig. 3). The median SWE maximum (in green) is close to 0 mm w.e. for all simulation chains and all warming levels. The 50-year return levels (in magenta-pink) usually decrease as a function of the warming level, except for a few simulation chains for which large SWE maxima are still simulated at high warming levels, e.g. for the simulation chain MPI-ESM-LR / REMO2009 (fourth line). For this NUTS-3 region, the preferred model is always pIG. Overall, return levels obtained with the different nonstationary models are similar.

Fig. 5 provides the same illustration for the NUTS-3 Osttirol where larger SWE maxima are simulated. The 50-year return levels also tend to decrease as a function of the warming level. In 5 out of 9 cases, the preferred model is GA, in 3 cases, the GEV model is preferred, and the GUM model is selected for one simulation chain. The different models generally produce similar quantiles.

### 4.2. Results at the European scale

Fig. 6 shows the best models for the different NUTS-3 regions when there is a dominant model, i.e. a model that has been selected more often than the other ones among the different simulation chains. Most of the NUTS-3 regions at low elevations (below 1000 m) and low latitudes (below 55°N) experience years without snow accumulation and the
<figure><img src="image_7.png" /><figcaption>Fig. 4. Illustration of return levels for the NUTS-3 Loire-Atlantique, France. Each column corresponds to a nonstationary model, and each row to a simulation chain. Return levels are shown for different annual exceedance probabilities as a function of the global warming level. Red boxes indicate the preferred model. For each simulation chain, annual SWE maxima are shown with black points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

models that have been fitted are either pEXP, pGA, or pIG. Interestingly, the selected models show great spatial coherence. pIG is preferred in Mediterranean regions (e.g. Italy, Greece) with mild climates, and NUTS-3 regions with oceanic climates (west of France, Ireland, United Kingdom, Belgium, Netherlands). pGA model is often preferred in other low-latitude or low-elevation European regions with a more continental climate. GA model is often selected for NUTS-3 regions with high elevations or high latitudes. For these NUTS-3 regions, it is difficult to distinguish a particular spatial pattern for the selected models. However, GA and GG models tend to be preferred below 1000 m and GUM and GEV models are evenly selected for NUTS-3 regions where the mean

elevation is higher than 500 m (see Fig. S1 in the Supplement).

Section S2 in the Supplement shows the results of the Shapiro-Wilk test applied to the normalized residuals of the selected nonstationary model. These goodness-of-fit tests indicate satisfying results and potential overfitting issues with models GG and GEV. In addition, we also provide estimates for the shape parameter ξ of the model GEV, which are generally smaller than 0 (see Fig. S3 in the Supplement), indicating that the distribution of SWE maxima, for a given global warming level, has an upper bound (i.e. a value that cannot be exceeded even for large return periods).

Fig. 7 shows the absolute changes of the proportion of zero SWE
<figure><img src="image_8.png" /><figcaption>Fig. 5. Illustration of return levels for the NUTS-3 Osttirol, Austria. Each column corresponds to a nonstationary model, and each row to a simulation chain. Return levels are shown for different annual exceedance probabilities as a function of the global warming level. Red boxes indicate the preferred model. For each simulation chain, annual SWE maxima are shown with black points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

maxima $p(w)$ between global warming levels of +3 °C and + 1 °C, averaged over the different simulation chains. A positive difference indicates an increasing proportion of years without snow and is greater than 10% for most low-elevation NUTS-3 regions with Mediterranean or oceanic climates. At +3 °C, most of these regions will have between 20 % and 40 % more years without snow compared to the current climate.

Fig. 8 shows the relative difference of 50-year return levels between global warming levels of +3 °C and + 1 °C. Return levels are strongly decreasing under a warmer climate. These decreases are generally moderate (between 0 % and 15 %) for NUTS-3 regions at high latitudes (Scandinavia, Iceland) and/or high elevations (in the Alps) and stronger (between −45 % and − 75 %) at low latitudes, along the Mediterranean coasts. Large decreases are also projected for many medium-altitude mountains (Massif Central in France, Pyrénées in France and Spain, Apennines in Italy, and mountain ranges in Greece). Large relative decreases generally correspond to small absolute differences of less than 40 mm w.e. (Fig. S4 in the Supplement).

4.3. *Uncertainty estimates*

We account for two types of uncertainty. The first source of uncer-

tainty corresponds to the sampling uncertainty. Indeed, the parameters of the different models are estimated using a limited amount of data, and this uncertainty is quantified using a nonparametric bootstrap approach. In the period 1951–2099, 149 years are resampled with replacement a hundred times to provide new samples of SWE maxima and corresponding global warming levels. As indicated in Section 3, this leads to 100 bootstrap estimates of 50-year return levels for each simulation chain. The second type of uncertainty corresponds to the model uncertainty, i.e. the fact that different climate models lead to different future climate conditions. This second type of uncertainty is assessed using different simulation chains. For any quantity (return levels, or relative changes of return levels), an overall estimate of these two uncertainties is provided by the quantiles 0.05 and 0.95 of the 100 × 9 estimates. The difference between these two quantiles summarizes the magnitude of the uncertainties.

Fig. 9 illustrates the uncertainty associated with 50-year return level estimates as a function of the global warming level, for the illustrative NUTS-3 regions. These curves can be compared to the annual SWE maxima (black points) and are expected to be exceeded roughly every 50 years on average (i.e. three times during the period 1951–2099 covered by the simulation chains). For each simulation chain (first nine
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Selected nonstationary models for the mean elevation of each NUTS-3 region, when a model has been selected more times than the others among the nine simulation chains. The NUTS-3 regions are indicated as “No dominant model” otherwise. No model is fitted when the series of SWE maxima contains more than 90 % of zeros (dark gray).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Average difference of the proportion of zero SWE maxima between global warming levels of +3 °C and +1 °C, i.e. $\bar{p}(3) - \bar{p}(1)$ where $\bar{p}(w)$ is the average estimate of $p(w)$ obtained for the different simulation chains. A positive difference indicates an increasing proportion of years without snow.</figcaption></figure>

lines), the thin colored curves correspond to the bootstrap estimates and vary much more for the NUTS-3 Loire-Atlantique than for the NUTS-3 Osttirol. Indeed, in Loire-Atlantique, the number of nonzero SWE maxima is strongly dependent on the GCM/RCM used to produce the simulation. For almost all simulation chains, a few SWE maxima values influence the return level estimates, and some bootstrap estimates seem erratic. For the simulation chain CNRM-CM5 / RCA4 (in brown), one large SWE maximum of about 20 mm w.e. at a global warming level of +3.8 °C leads to two populations of bootstrap estimates. One population where the return levels slowly decrease as a function of the warming level, similar to the estimate obtained without bootstrap, i.e. with all the SWE maxima (thick curve), and one population where return level estimates are close to zero for warming levels above +3 °C. Bootstrap estimates are more stable for the NUTS-3 Osttirol, which are usually close to the best estimate. The last line shows the quantiles 0.05 and 0.95

obtained from the 100 × 9 estimates and the mean return level estimates, averaged over the different simulation chains. For the NUTS-3 Loire-Atlantique, mean return levels are small (below 20 mm w.e.), and the uncertainty interval's width is rather large compared to the NUTS-3 Osttirol.

**Fig. 10** provides the map of the 90 % interquantile range (difference in % between quantiles 0.95 and 0.05) for relative differences of 50-year return levels between global warming levels of +3 °C and +1 °C. These uncertainties are generally smaller than 40 % for NUTS-3 regions where SWE maxima are nonzero (Fig. 7). In low-elevated regions, these uncertainties are larger than 80 %, the corresponding quantiles 0.05 and 0.95 indicating respectively marked decreases (between −60 % and −75 %) and small to large increases (Figs. S4 and S5 in the Supplement). However, it must be noted that these changes correspond to small absolute values, as shown in Fig. S4 and illustrated in Fig. 9.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Relative difference (in %) of 50-year return levels at warming levels of +3 °C and +1 °C, i.e. $100 \times (RL(3)/RL(1) - 1)$. Negative relative differences are shown with brown colors and indicate decreasing 50-year return levels in a warmer climate. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

# 5. Discussion

## 5.1. MTMSI dataset

The results shown in this study rely on the MTMSI dataset that uses refined regional climate projections, adjusted with a pan-European reanalysis, to run the Crocus snowpack model. Monteiro and Morin (2023) recently showed that the MTMSI dataset compares well with observations for the past decades, in terms of mean seasonal snow depth values in the European Alps. Kamir et al. (2024) extends this evaluation to annual snow depth maxima and for a wider geographic domain including low-lying and higher-latitude areas. This former study shows a satisfying agreement between simulated and observed annual maxima of snow depth. The evaluation of SWE simulations is more difficult due to the rarity of in-situ SWE observations. Krinner et al. (2018) shows that Crocus model generally provides a very satisfying reproduction of SWE observations in many regions of the world and, in particular, for sites in France, Switzerland, and Finland.

The SWE maxima used in this study are available at the scale of the NUTS-3 regions and the corresponding mean elevation. This spatial scale has the advantage of being compatible with other socio-economic indicators but has some limitations regarding local climate features (e.g. coastal effects, Chen et al., 2021) or the elevation dependency (Morin et al., 2021). Indeed, at high elevations (e.g. above 1500 m), it is sometimes difficult to find nearby pixels with a close elevation in the reference dataset UERRA. This can create discrepancies between neighboring NUTS-3 regions for a given elevation. This limitation motivates the choice to illustrate return level estimates only for the mean elevation of the NUTS-3 regions and not for higher elevations, even if they are available in the MTMSI dataset.

Another shortcoming concerns the representation of convective events by the RCMs used to produce the MTMSI dataset described in Section 2. Indeed, these regional climate models consider parameterized convection schemes that are not able to produce local and intense convective events. Recently, convection-permitting regional climate models have been shown to produce more realistic climate simulations (more intense precipitation and reduced wet-hour frequency) in mountainous areas (Ban et al., 2021; Monteiro et al., 2022). However, these simulations usually cover short future time windows that do not enable statistical analyses such as those introduced in the present study,

in particular for extreme values analyses.

Last, we mention that, in contrast to the discrepancy identified between RCM and GCM simulated temperature changes over Europe in summer, there is a good overall consistency for the winter season, which is the main seasonal focus of our study (Taranu et al., 2023; Ribes et al., 2022). The direct use of global warming levels applied to the CMIP5/EURO-CORDEX GCM/RCM pairs is deemed an acceptable choice, even if further work is required to ensure better consistency between global and regional warming levels across GCM and RCM model chains.

## 5.2. Statistical framework

In this study, we introduce a statistical framework that applies several nonstationary models to annual maxima of daily SWE values, using global warming levels as covariates for the parameters of different distributions. The statistical framework includes a statistical criterion (AIC) to select the best model, a goodness-of-fit test to check the validity of the selected model, and a bootstrap method to estimate the uncertainty related to the parameter estimation (a.k.a. sampling uncertainty). While there is a rich literature concerning the fitting of precipitation extremes (Papalexiou and Koutsoyiannis, 2013), the statistical assessment of SWE extremes raises additional difficulties:

* In regions where snowfalls are frequent, annual maxima of daily SWE values correspond to the peak of snow accumulations during winter. Daily SWE time series are thus strongly autocorrelated. The GEV distribution is still applicable if annual maxima can be considered as independent (Leadbetter et al., 1983) and it has been successfully applied to snow depth extremes (Blanchet and Lehning, 2010; Schellander and Hell, 2018). However, the GEV distribution, as an approximation to the distribution of very large block maxima, might not be as efficient in these cases (Coles, 2001, Section 5.3.1) and could explain why the model GEV is not the preferred model in most cases.

* In regions where snowfall is rare, other problems are encountered, mainly due to a great majority of zero daily SWE values. Annual SWE maxima are thus the maxima over a few nonzero values and are sometimes equal to zero. Under these conditions, discrete-continuous distributions must be used with the risk of having not many zero SWE maxima for low warming levels or, on the contrary,
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Illustration of the uncertainty associated with 50-year return level estimates as a function of the global warming level, for the illustrative NUTS-3 regions: Loire-Atlantique, France (0 m) on the left column and Osttirol, Austria (2000 m) on the right column. The first nine lines correspond to the different simulation chains. The different thin curves show the 100 bootstrap estimates and the thick curve shows the estimate obtained using all SWE maxima (i.e. not bootstrapped). The last line shows the overall estimates. The interval in gray corresponds to the quantiles 0.05 and 0.95 obtained from the 100 × 9 estimates of 50-year return levels represented by all the thin curves of the first nine lines. The thick black curve is the mean return level estimate, averaged over the different simulation chains. For each simulation chain, annual SWE maxima are shown with black points.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** 90 % interquantile range (difference in % between quantiles 0.95 and 0.05) associated with the relative differences of 50-year return levels at warming levels of +3 °C and +1 °C. This uncertainty combines estimates from 100 bootstrap replicates and the nine simulation chains.</figcaption></figure>

not many positive SWE maxima for large warming levels, which can be problematic for the estimation of the parameter trends.

In addition, in the literature on precipitation extremes, many distributions with three or more parameters have been tested, for example, Log-Pearson type-3, Generalized Logistic (Benyahya et al., 2013), and many others (see Table S1 in Wilks, 1993). However, the flexibility of these distributions comes at the cost of a larger uncertainty for return level estimates (Klemes, 2000a, 2000b). Similarly, simple linear links are considered for the evolution of the distribution parameters, as more flexible relationships (e.g. smoothing splines) often lead to dubious forms of the evolution of the parameters (e.g. local peaks).

### 5.3. Comparison to existing studies

In Europe, countries have adopted different approaches to estimate 50-year return levels of snow load involved in the Eurocodes, in a stationary context (i.e. ignoring the impact of climate change). Several factors can explain the inhomogeneous estimates that are obtained (Croce et al., 2019) and can be summarized as follows:

* **Snow data:** The data used to estimate the 50-year return levels of snow load are usually in-situ snow measurements. Snow depth measurements are often available on longer periods and at a higher spatial density, compared to SWE measurements (see, e.g. Schellander et al., 2021, in Austria). However, the relationship between 50-year return levels of snow depth and snow load is not direct and assumptions about the snow density must be made. For example, in France, a constant snow density of 150 km m⁻³ is assumed which tends to underestimate snow loads (Le Roux et al., 2020). In this study, as indicated above, snow characteristics are provided by the Crocus snowpack model which considers many processes (accumulation, thaw-freeze, melt, compaction, etc.) and thus avoids additional assumptions about snow depth/snow load relationships.

* **Spatialization:** Each country also defines different zones that are considered climatically homogeneous. Inside these large zones, different relationships between snow load and elevation are established which can create differences across borders (in addition to meteorological effects), see Croce et al. (2019). In this study, NUTS-3 regions and mean elevation of each region define the spatial scale and provide a relatively fine resolution over Europe in comparison to the standard climatic zones defined nationally. For example, in

France, 11 climatic zones were considered in the Eurocodes (Biétry, 2005) compared to the 96 NUTS-3 regions in metropolitan France.

* **Statistical method:** Many distributions have been considered for ground snow loads (Gumbel, Log-normal, Weibull, Croce et al., 2019) and the “best” distribution varies according to the elevation and the geographical location. The Gumbel distribution seems to be favored in some recent applications (Croce et al., 2018b, 2021) but several studies also motivate the application of the GEV distribution (which includes the Gumbel distribution as a special case) to extreme snow variables (Blanchet and Lehning, 2010; Gaume et al., 2013; Schellander and Hell, 2018; Le Roux et al., 2020, 2021; Schellander et al., 2021). Since there is no single distribution adapted to all the climatic conditions in Europe, and the various elevations, the statistical method presented in this study provides a relevant compromise where various candidate distributions are favored depending on their performance. However, this approach does not consider an explicit relationship between the snow load maxima and the elevation on the contrary to previous studies focusing on mountainous areas (see, e.g. Schellander et al., 2021).

Most of the aforementioned studies do not consider a changing climate. Croce et al. (2018b, 2021) consider snow outputs from the RCMs and fit Gumbel distributions to future time windows, which has two main limitations: i) snow variable simulations by the RCMs have severe shortcomings (Matiu et al., 2020), and ii) the stationarity assumption for consecutive annual maxima over large time windows (e. g. 40 years in Croce et al., 2021) is questionable because trends inside the time window can be large and are usually not verified (Kharin et al., 2013). On the other side, small time windows (e.g. 20 years) lead to small samples of annual maxima, and thus to uncertain estimates. Uncertain snow load estimates for the reference and the future periods can be problematic when relative changes are estimated (a.k.a. factor of changes Croce et al., 2021), especially if the estimates for the reference periods are close to zero. Nonstationary models proposed in this study can be considered as more robust in this regard since they are fitted to the entire series of maxima (150 years) such that long-term trends are more easily identified.

## 6. Conclusion

In this study, a statistical framework is proposed to estimate 50-year return levels of daily SWE values as a function of global warming levels.
Different nonstationary models are tested according to the presence of zeros in the projections of annual SWE maxima. This flexible approach can handle various climates and areas with a wide range of mean elevations (e.g. low-elevation and mountainous areas). Return levels are provided at the NUTS-3 scale, using a multi-model ensemble of snow cover projections that have been produced on a pan-European domain. These estimates can be transformed (multiplication by the gravitational acceleration constant) into 50-year return levels of snow load, which is the primary component of the Eurocodes concerning the action of snow loads on structures (Sanpaolesi, 1998). In a + 3 °C world, extreme snow loads are projected to decrease in all of Europe compared to the current climate. These decreases are, however, moderate at high latitudes (Scandinavia, Iceland) and at high elevations (e.g. in the Alps). At low elevations and low latitudes (e.g. along the coasts of the Mediterranean area), these decreases can be strong (−50 % to −75 %) but are uncertain due to the sensitivity of these relative changes to small estimates under the current climate. This study also highlights large uncertainties due to disagreements between the climate models used to produce the simulated SWE maxima in a future climate.

Snow load estimates are not harmonized at the European scale and each country provides a different map of snow loads (Croce et al., 2019). This study provides a methodology to obtain homogeneous estimates at a fine regional scale in Europe, taking advantage of regional climate simulations that explicitly account for the effect of elevation, as well as oceans/seas/lakes, and represent the main meteorological processes. The snowpack model used represents the main physical processes operating in the snowpack (e.g. thermal diffusion, phase change, metamorphism, compaction, snowmelt) and is a powerful alternative to the standard approach that consists of assessing snow depth data with a prescribed snow density and snow load/altitude relationships (Croce et al., 2019) or the direct use of snow cover state simulated by the RCMs (Matiu et al., 2020; Croce et al., 2021).

## Code and data availability

Simulated temperature from CMIP5 models can be accessed at https://esgf-node.ipsl.upmc.fr/projects/cmip5-ipsl/. Access to HadCRUT5 dataset is detailed in Morice et al. (2021). The MTMSI dataset (Morin et al., 2021) is available on the Copernicus Data Store following this url (last access: 02 December 2024): https://cds.climate.copernicus.eu/datasets/sis-tourism-snow-indicators?tab=overview under the Copernicus license. All computations were performed with R software version 4.1.1. The codes are available from a repository (https://doi.org/10.5281/zenodo.14259156, Evin, 2024) which includes datasets and scripts to perform the statistical analyses and to create Figs. 3, 4, 5, and 9 of the present article.

## Authorship

We confirm that the manuscript has been read and approved by all named authors including the order of authors listed.

## CRediT authorship contribution statement

G. Evin: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. E. Le Roux: Writing – review & editing, Methodology. E. Kamir: Writing – review & editing, Methodology. S. Morin: Writing – review & editing, Resources, Project administration, Methodology, Funding acquisition, Conceptualization.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests.

Elisa KAMIR reports financial support was provided by European Commission. If there are other authors, they 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 work has been performed in the context of the “Framework Partnership Agreement for Copernicus User Uptake” (FPCUP) action “Estimation of snow load data using Copernicus and in-situ data”. Elisa Kamir has been supported by FPCUP which is financed by the European Commission under the FPA no. 275/G/GRO/COPE/17/10042. The authors thank Raphaëlle Samacoïts from Météo-France for computing and providing the SWE maxima dataset. This work also contributes to the outcomes of the PEPR TRACCS (projects ANR-22-EXTR-0004 and ANR-22-EXTR-0005).

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.coldregions.2025.104424.

## Data availability

Data will be made available on request.

## References

Ban, N., Caillaud, C., Coppola, E., Pichelli, E., Sobolowski, S., Adinolfi, M., Ahrens, B., Alias, A., Anders, I., Bastin, S., Belušić, D., Berthou, S., Brisson, E., Cardoso, R.M., Chan, S.C., Christensen, O.B., Fernández, J., Fita, L., Frisius, T., Gašparac, G., Giorgi, F., Goergen, K., Haugen, J.E., Hodnebrog, O., Kartsios, S., Katragkou, E., Kendon, E.J., Keuler, K., Lavin-Gullon, A., Lenderink, G., Leutwyler, D., Lorenz, T., Maraun, D., Mercogliano, P., Milovac, J., Panitz, H.J., Raffa, M., Remedio, A.R., Schär, C., Soares, P.M.M., Srnek, L., Steensens, B.M., Stocchi, P., Tölle, M.H., Truhetz, H., Vergara-Temprado, J., de Vries, H., Warrach-Sagi, K., Wulfmeyer, V., Zander, M.J., 2021. The first multi-model ensemble of regional climate simulations at kilometer-scale resolution, part I: evaluation of precipitation. Clim. Dyn. 57, 275–302. URL: https://doi.org/10.1007/s00382-021-05708-w.

Benyahia, L., Gachon, P., St-Hilaire, A., Laprise, R., 2013. Frequency analysis of seasonal extreme precipitation in southern Quebec (Canada): an evaluation of regional climate model simulation with respect to two gridded datasets. Hydrol. Res. 45, 115–133. URL: https://doi.org/10.2166/nh.2013.066.

Biétry, J., 2005. Charges de neige au sol en France : proposition de carte révisée. Technical Report. CSTB, Saint-Herblain, France.

Blanchet, J., Lehning, M., 2010. Mapping snow depth return levels: smooth spatial modeling versus station interpolation. Hydrol. Earth Syst. Sci. 14, 2527–2544. URL: http://www.hydrol-earth-syst-sci.net/14/2527/2010/ https://doi.org/10.5194/hess-14-2527-2010.

Burnham, K.P., Anderson, D.R. (Eds.), 2004. Model Selection and Multimodel Inference. Springer, New York, NY. URL: http://link.springer.com/10.1007/b97636 https://doi.org/10.1007/b97636.

Chen, G., Wang, W.C., Cheng, C.T., Hsu, H.H., 2021. Extreme snow events along the coast of the northeast united states: Potential changes due to global warming. J. Clim. 34, 2337–2353. URL: https://journals.ametsoc.org/view/journals/clim/34/6/JCLI-D-20-0197.1.xml https://doi.org/10.1175/JCLI-D-20-0197.1. place: Boston MA, USA Publisher: American Meteorological Society.

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

Cooley, D., 2013. Return periods and return Levels under climate Change. In: AghaKouchak, A., Easterling, D., Hsu, K., Schubert, S., Sorooshian, S. (Eds.), Extremes in a Changing Climate: Detection, Analysis and Uncertainty. Springer Netherlands, Dordrecht. Water Science and Technology Library, pp. 97–114. URL: https://doi.org/10.1007/978-94-007-4479-0 4.

Croce, P., Formichi, P., Landi, F., Marsili, F., 2018a. Climate change: Impact on snow loads on structures. Cold Reg. Sci. Technol. 150, 35–50. URL: http://www.sciencedirect.com/science/article/pii/S0165232X17304834 https://doi.org/10.1016/j.coldregions.2017.10.009.

Croce, P., Formichi, P., Landi, F., Mercogliano, P., Bucchignani, E., Dosio, A., Dimova, S., 2018b. The snow load in Europe and the climate change. Climate. Risk Manage. 20, 138–154. URL: http://www.sciencedirect.com/science/article/pii/S2212096317301754 https://doi.org/10.1016/j.crm.2018.03.001.

Croce, P., Formichi, P., Landi, F., Marsili, F., 2019. Harmonized European ground snow load map: Analysis and comparison of national provisions. Cold Reg. Sci. Technol. 168, 102875. URL: https://www.sciencedirect.
com/science/article/pii/S0165232X19300400 https://doi.org/10.1016/j.coldregions.2019.102875.

Croce, P., Formichi, P., Landi, F., 2021. Extreme Ground Snow Loads in Europe from 1951 to 2100. Climate 9, 133. URL: https://www.mdpi.com/2225-1154/9/9/133 https://doi.org/10.3390/cli0900133. Number: 9
Publisher: Multidisciplinary Digital Publishing Institute.

Efron, B., Tibshirani, R., 1994. An Introduction to the Bootstrap. Chapman and Hall/CRC, New York. https://doi.org/10.1201/9780429246593.

Evin, G., 2024. Estimating Changes in Extreme Snow Load Conditionally to Global Warming Levels. URL: https://doi.org/10.5281/zenodo.14259157. zenodo.

Evin, G., Ribes, A., Corre, L., 2024. Assessing CMIP6 Uncertainties at Global Warming Levels. Climate Dynamics. https://doi.org/10.1007/s00382-024-07323-x. URL:

Gaume, J., Eckert, N., Chambon, G., Naaim, M., Bel, L., 2013. Mapping extreme snowfalls in the French Alps using max-stable processes. Water Resour. Res. 49, 1079–1098. https://doi.org/10.1002/wrcr.20083.

IPCC, 2018. IPCC special report on the impacts of global warming of 1.5 °C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty. [V. Masson-Delmotte, P. Zhai, H. O. Pörtner, D. Roberts, J. Skea, P.R. Shukla, A. Pirani, W. Moufouma-Okia, C. Péan, R. Pidcock, S. Connors, J. B. R. Matthews, Y. Chen, X. Zhou, M. I. Gomis, E. Lonnoy, T. Maycock, M. Tignor, T. Waterfield]. URL: http://www.ipcc.ch/report/sr15/. 151.

IPCC, 2021. In: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S.L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M.I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J.B.R., Maycock, T.K., Waterfield, T., Yelekçi, O., Yu, R., Zhou, B. (Eds.), Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. https://doi.org/10.1017/9781009157896.

IPCC, 2022. High Mountain areas. In: the Ocean and Cryosphere in a changing climate: special Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, pp. 131–202. URL: https://www.cambridge.org/core/books/ocean-and-cryosphere-in-a-changing-climate/high-mountain-areas/D4285C8P6E9C47D040525BD2140575A https://doi.org/10.1017/9781009157964.004.

Kamir, E., Morin, S., Evin, G., Gehring, P., Wichura, B., 2024. Estimation of snow load data using Copernicus and in-situ data. In preparation.

Kharin, V.V., Zwiers, F.W., Zhang, X., Wehner, M., 2013. Changes in temperature and precipitation extremes in the CMIP5 ensemble. Clim. Chang. 119, 345–357. URL: https://doi.org/10.1007/s10584-013-0705-8.

Kim, H., Kim, S., Shin, H., Heo, J.H., 2017. Appropriate model selection methods for nonstationary generalized extreme value models. J. Hydrol. 547, 557–574. URL: https://www.sciencedirect.com/science/article/pii/S0022169417300756 https://doi.org/10.1016/j.jhydrol.2017.02.005.

Klemes, V., 2000a. Tall Tales about Tails of Hydrological Distributions. J. Hydrol. Eng. 5, 227–231. https://doi.org/10.1061/(ASCE)1084-0699(2000)5:3(227).

Klemes, V., 2000b. Tall Tales about Tails of Hydrological Distributions. J. Hydrol. Eng. 5, 232–239. https://doi.org/10.1061/(ASCE)1084-0699(2000)5:3(232).

Kotlarski, S., Keuler, K., Christensen, O.B., Colette, A., Déqué, M., Gobiet, A., Goergen, K., Jacob, D., Lüthi, D., van Meijgaard, E., Nikulin, G., Schär, C., Teichmann, C., Vautard, R., Warrach-Sagi, K., Wulfmeyer, V., 2014. Regional climate modeling on European scales: a joint standard evaluation of the EURO-CORDEX RCM ensemble. Geosci. Model Dev. 7, 1297–1333. URL: https://www.geosci-model-dev.net/7/1297/2014/ https://doi.org/10.5194/gmd-7-1297-2014.

Kotlarski, S., Gobiet, A., Morin, S., Olefs, M., Rajczak, J., Samacoits, R., 2023. 21st Century alpine climate change. Clim. Dyn. 60, 65–86. URL: https://doi.org/10.1007/s00382-022-06303-3.

Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S., Clark, M., Hall, A., Rott, H., Brutel-Vuilmet, C., Kim, H., Ménard, C.B., Mudryk, L., Thackeray, C., Wang, L., Arduini, G., Balsamo, G., Bartlett, P., Boike, J., Boone, A., Chéruy, F., Colin, J., Cuntz, M., Dai, Y., Decharne, B., Derry, J., Ducharne, A., Dutra, E., Fang, X., Fierz, C., Ghattas, J., Gusev, Y., Haverd, V., Kontu, A., Lafaysse, M., Law, R., Lawrence, D., Li, W., Marke, T., Marks, D., Ménégoz, M., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Raleigh, M.S., Schaedler, G., Semenov, V., Smirnova, T.G., Stacke, T., Strasser, U., Svenson, S., Turkov, D., Wang, T., Wever, N., Yuan, H., Zhou, W., Zhu, D., 2018. ESM-SnowMIP: assessing snow models and quantifying snow-related climate feedbacks. Geosci. Model Dev. 11, 5027–5049. URL: https://gmd.copernicus.org/articles/11/5027/2018/ https://doi.org/10.5194/gmd-11-5027-2018. publisher: Copernicus GmbH.

Le Roux, E., Evin, G., Eckert, N., Blanchet, J., Morin, S., 2020. Non-stationary extreme value analysis of ground snow loads in the French Alps: a comparison with building standards. Nat. Hazards Earth Syst. Sci. 20, 2961–2977. URL: https://nhess.copernicus.org/articles/20/2961/2020/ https://doi.org/10.5194/nhess-20-2961-2020. publisher: Copernicus GmbH.

Le Roux, E., Evin, G., Eckert, N., Blanchet, J., Morin, S., 2021. Elevation-dependent trends in extreme snowfall in the French Alps from 1959 to 2019. Cryosphere 15, 4335–4356. URL: https://tc.copernicus.org/articles/15/4335/2021/ https://doi.org/10.5194/tc-15-4335-2021. Publisher: Copernicus GmbH.

Le Roux, E., Evin, G., Eckert, N., Blanchet, J., Morin, S., 2022. A non-stationary extreme-value approach for climate projection ensembles: application to snow loads in the French Alps. Earth Syst. Dynam. 13, 1059–1075. URL: https://esd.copernicus.org/articles/13/1059/2022/ https://doi.org/10.5194/esd-13-1059-2022. Publisher: Copernicus GmbH.

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

Marty, C., Schlögl, S., Bavay, M., Lehning, M., 2017. How much can we save? Impact of different emission scenarios on future snow cover in the Alps. Cryosphere 11, 517–529. URL: https://www.the-cryosphere.net/11/517/2017/ https://doi.org/10.5194/tc-11-517-2017.

Matiu, M., Petitta, M., Notarnicola, C., Zebisch, M., 2020. Evaluating Snow in EURO-CORDEX Regional climate Models with Observations for the European Alps: Biases and their Relationship to Orography, Temperature, and Precipitation Mismatches. Atmosphere 11, 46. URL: https://www.mdpi.com/2073-4433/11/1/46 https://doi.org/10.3390/atmos11010046. Number: 1 Publisher: Multidisciplinary Digital Publishing Institute.

Monteiro, D., Morin, S., 2023. Multi-decadal analysis of past winter temperature, precipitation and snow cover data in the European Alps from reanalyses, climate models and observational datasets. Cryosphere 17, 3617–3660. URL: https://tc.copernicus.org/articles/17/3617/2023/ https://doi.org/10.5194/tc-17-3617-2023. Publisher: Copernicus GmbH.

Monteiro, D., Caillaud, C., Samacoits, R., Lafaysse, M., Morin, S., 2022. Potential and limitations of convection-permitting CNRM-AROME climate modelling in the French Alps. Int. J. Climatol. 42, 7162–7185. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/joc.7637. doi: 10.1002/joc.7637. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/joc.7637.

Morice, C.P., Kennedy, J.J., Rayner, N.A., Winn, J.P., Hogan, E., Killick, R.E., Dunn, R.J. H., Osborn, T.J., Jones, P.D., Simpson, I.R., 2021. An Updated Assessment of Near-Surface Temperature Change from 1850: the HadCRUT5 Data Set. J. Geophys. Res.-Atmos. 126, e2019JD032361. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/2019JD032361. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.2019JD032361.

Morin, S., Samacoits, R., François, H., Carmagnola, C.M., Abegg, B., Demiroglu, O.C., Pons, M., Soubeyroux, J.M., Lafaysse, M., Franklin, S., Griffiths, G., Kite, D., Hoppler, A.A., George, E., Buontempo, C., Almond, S., Dubois, G., Cauchy, A., 2021. Pan-European meteorological and snow indicators of climate change impact on ski tourism. Clim. Servic. 22, 100215. URL: https://www.sciencedirect.com/science/article/pii/S2405880721000030 https://doi.org/10.1016/j.ncliser.2021.100215.

Northrop, P.J., 2024. gamlssx: Generalized Additive Extreme Value Models for Location, Scale and Shape. URL: https://cran.r-project.org/package=gamlssx. r package version 1.0.1.

O'Donnell, F.C., Tingerthal, J.S., White, S., 2020. Estimation of Ground Snow Loads for Low-Latitude, High-Elevation Regions. Journal of Cold Regions Engineering 34, 04020008. URL: https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29CR.1943-5495.0000209, doi: 10.1061/(ASCE)CR.1943-5495.0000209. publisher: American Society of Civil Engineers.

O'Rourke, M.J., 1997. Snow load on buildings: Statistics of snowfall plus small-scale simulations help engineers design buildings that can withstand the loads of rooftop drifts. Am. Sci. 85, 64–70. URL: http://www.jstor.org/stable/27856692.

Ousset, I., Evin, G., Raynaud, D., Faug, T., 2023. Back analysis of a building collapse under snow and rain loads in a Mediterranean area. Nat. Hazards Earth Syst. Sci. 23, 3509–3523. URL: https://nhess.copernicus.org/articles/23/3509/2023/ https://doi.org/10.5194/nhess-23-3509-2023. Publisher: Copernicus GmbH.

Papalexiou, S.M., Koutsoyiannis, D., 2013. Battle of extreme value distributions: a global survey on extreme daily rainfall. Water Resour. Res. 49, 187–201. https://doi.org/10.1029/2012WR012557.

Papalexiou, S.M., Koutsoyiannis, D., Makropoulos, C., 2013. Hydrom extreme is extreme? An assessment of daily rainfall distribution tails. Hydrol. Earth Syst. Sci. 17, 851–862. URL: http://www.hydrol-earth-syst-sci.net/17/851/2013/ https://doi.org/10.5194/hess-17-851-2013.

Ribes, A., Boé, J., Qasmi, S., Duboisson, B., Douville, H., Terray, L., 2022. An updated assessment of past and future warming over France based on a regional observational constraint. Earth Syst. Dynam. 13, 1397–1415. URL: https://esd.copernicus.org/articles/13/1397/2022/ https://doi.org/10.5194/esd-13-1397-2022. Publisher: Copernicus GmbH.

Rigal, A., Azaïs, J.M., Ribes, A., 2019. Estimating daily climatological normals in a changing climate. Clim. Dyn. 53, 275–286. URL: https://doi.org/10.1007/s00382-018-4584-6.

Roston, P., 1995. Remark AS R94: a Remark on Algorithm AS 181: the W-test for Normality. J. R. Stat. Soc.: Ser. C: Appl. Stat. 44, 547–551. URL: https://www.jstor.org/stable/2986146 https://doi.org/10.2307/2986146. Publisher: [Wiley, Royal Statistical Society].

Sadovský, Z., Fasko, P., Mikulová, K., Peko, J., 2012. Exceptional snowfalls and the assessment of accidental snow loads for structural design. Cold Reg. Sci. Technol. 72, 17–22. URL: http://www.sciencedirect.com/science/article/pii/S0165232X11002655 https://doi.org/10.1016/j.coldregions.2011.12.003.

Sanpaoliels, L., 1998. Snow loading: scientific basis, problems and challenges. In: Progress in Structural Engineering and Materials, 1, pp. 443–451. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/pse.2260010414. doi: 10.1002/pse.2260010414. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/pse.2260010414.

Schellander, H., Hell, T., 2018. Modeling snow depth extremes in Austria. Nat. Hazards. https://doi.org/10.1007/s11069-018-3481-y. URL:

Schellander, H., Winkler, M., Hell, T., 2021. Towards a Reproducible Snow Load Map – An Example for Austria. In: Advances in Science and Research. Copernicus GmbH, pp. 135–144. URL: https://asr.copernicus.org.
org/articles/18/135/2021/ https://doi.org/10.5194/asr-18-135-2021 (iSSN: 1992-0628).

Soci, C., Bazile, E., Besson, F., Landelius, T., 2016. High-resolution precipitation re-analysis system for climatological purposes. Tellus A: Dynam. Meteorol. Oceanogr. 68, 29879. URL: https://doi.org/10.3402/tellusa.v68.29879.

Stasinopoulos, D.M., Rigby, R.A., 2008. Generalized Additive Models for Location Scale and Shape (GAMLSS) in R. J. Stat. Softw. 23, 1–46. URL: 10.18637/jss.v023.i07.

Strasser, U., 2008. Snow loads in a changing climate: new risks? Nat. Hazards Earth Syst. Sci. 8, 1–8. URL: https://nhess.copernicus.org/articles/8/1/2008/ https://doi.org/10.5194/nhess-8-1-2008. Publisher: Copernicus GmbH.

Taranu, I.S., Somot, S., Alias, A., Boé, J., Delire, C., 2023. Mechanisms behind large-scale inconsistencies between regional and global climate model-based projections over Europe. Clim. Dyn. 60, 3813–3838. URL: https://doi.org/10.1007/s00382-022-06540-6.

UNFCCC, 2015. Paris Agreement - Paris Climate Change Conference - November 2015 - COP 21. URL: https://unfccc.int/documents/184656.

Verfaillie, D., Déqué, M., Morin, S., Lafaysse, M., 2017. The method ADAMONT v1.0 for statistical adjustment of climate projections applicable to energy balance land

surface models. Geosci. Model Dev. 10, 4257–4283. https://www.geosci-model-dev.net/10/4257/2017/. https://doi.org/10.5194/gmd-10-4257-2017.

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 1500m altitude in the Chartreuse mountain range, Northern French Alps. Cryosphere 12, 1249–1271. URL: https://www.the-cryosphere.net/12/1249/2018/ https://doi.org/10.5194/tc-12-1249-2018.

Wilks, D.S., 1993. Comparison of three-parameter probability distributions for representing annual extreme and partial duration precipitation series. Water Resour. Res. 29, 3543–3549. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/93WR01710. doi: 10.1029/93WR01710. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/93WR01710.

Xavier, A.C.F., Blain, G.C., Morais, M.V.B.D., Sobierajski, G.D.R., 2019. Selecting “the best” nonstationary Generalized Extreme Value (GEV) distribution: on the influence of different numbers of GEV-models. Bragantia 78, 606–621. URL: https://www.scielo.br/j/brag/a/3ZPSYGD6ZXpfSgVSQSWQcNQ/?lang=en https://doi.org/10.1590/1678-4499.20180408. publisher: Instituto Agronômico de Campinas.