the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Elevation-dependent trends in extreme snowfall in the French Alps from 1959 to 2019

### Guillaume Evin

### Nicolas Eckert

### Juliette Blanchet

### Samuel Morin

Climate change projections indicate that extreme snowfall is expected to increase in cold areas, i.e., at high latitudes and/or high elevation, and to decrease in warmer areas, i.e., at mid-latitudes and low elevation. However, the magnitude of these contrasting patterns of change and their precise relations to elevation at the scale of a given mountain range remain poorly known. This study analyzes annual maxima of daily snowfall based on the SAFRAN reanalysis spanning the time period 1959–2019 and provided within 23 massifs in the French Alps every 300 m of elevation. We estimate temporal trends in 100-year return levels with non-stationary extreme value models that depend on both elevation and time. Specifically, for each massif and four elevation ranges (below 1000, 1000–2000, 2000–3000, and above 3000 m), temporal trends are estimated with the best extreme value models selected on the basis of the Akaike information criterion. Our results show that a majority of trends are decreasing below 2000 m and increasing above 2000 m. Quantitatively, we find an increase in 100-year return levels between 1959 and 2019 equal to +23 % ($+\mathrm{32}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$) on average at 3500 m and a decrease of −10 % ($-\mathrm{7}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$) on average at 500 m. However, for the four elevation ranges, we find both decreasing and increasing trends depending on location. In particular, we observe a spatially contrasting pattern, exemplified at 2500 m: 100-year return levels have decreased in the north of the French Alps while they have increased in the south, which may result from interactions between the overall warming trend and circulation patterns. This study has implications for natural hazard management in mountain regions.

Extreme snowfall can generate casualties and economic damage. For instance, it can cause major natural hazards (avalanche, winter storms) that might be intensified with high winds and freezing rain. Heavy snowfall can also disrupt transportation (road, rail, and air traffic), tourism, electricity (power lines), and communication systems with a significant impact on economic services (Changnon, 2007; Blanchet et al., 2009). Subsequently, snow overloading can lead to the collapse of buildings such as a shed, a greenhouse, or something as large as an exhibition hall (Strasser, 2008). It remains a counterintuitive phenomenon that extreme snowfall can increase in a warming climate, at least transiently, i.e., as long as local temperatures are cold enough (Frei et al., 2018). Therefore, to adapt protective measures, it is crucial to determine temporal trends in extreme snowfall for various areas (regions, elevations) and timescales, and to understand the underlying causes of these trends.

Extreme snowfall stems from extreme precipitation occurring in a range of optimal temperatures slightly below 0 ^{∘}C according to O'Gorman (2014). This optimal range of temperatures favors both high precipitation intensities and percentages of precipitation falling as snow close to 100 %.
Thus, changes in extreme snowfall depend on a trade-off between trends in extreme precipitation and changes in the probability of experiencing temperature in this optimal range.

On a global scale, extreme precipitation is expected to increase with the augmentation of global mean temperature.
Specifically, the most intense precipitation rates are theoretically expected to roughly increase at a rate of $\mathrm{7}\phantom{\rule{0.125em}{0ex}}\mathrm{\%}\phantom{\rule{0.125em}{0ex}}{}^{\circ}{\mathrm{C}}^{-\mathrm{1}}$,
i.e., $\mathrm{7}\phantom{\rule{0.125em}{0ex}}\mathrm{\%}\phantom{\rule{0.125em}{0ex}}{}^{\circ}{\mathrm{C}}^{-\mathrm{1}}$ of global mean warming,
due to an increase in maximum atmospheric water vapor content according to the Clausius–Clapeyron relationship (O'Gorman and Muller, 2010).
In practice, the observed global mean temperature scaling for annual maxima of 1 d precipitation is $\mathrm{6.6}\phantom{\rule{0.125em}{0ex}}\mathrm{\%}\phantom{\rule{0.125em}{0ex}}{}^{\circ}{\mathrm{C}}^{-\mathrm{1}}$ (Sun et al., 2021).
On the other hand, the probability of experiencing temperature in the optimal range for extreme snowfall is expected to decrease in warm areas, i.e., mid-latitude and low-elevation regions, as temperatures are expected to shift away from 0 ^{∘}C.
However, this probability may increase in cold areas, i.e., high-latitude and/or high-elevation regions, where temperatures are expected to shift toward 0 ^{∘}C while remaining below 0 ^{∘}C (Frei et al., 2018).

In the European Alps, past observations show both that the warming rate is larger than the global warming rate and that trends in extreme precipitation depend on the season and on the region.
Indeed, past trends in mean annual surface temperature point to a temperature increase in high mountain regions of central Europe, with a warming rate ranging from 0.15 to 0.35 ^{∘}C per decade since 1960 (Fig. 2.2 of IPCC, 2019)
against a range from 0.08 to 0.14 ^{∘}C per decade since 1951
for the global warming rate (IPCC, 2013).
Furthermore, past trends in daily maxima of precipitation largely depend on the season (Fig. 7 of Ménégoz et al., 2020).
In winter, daily maxima precipitation (which may generate extreme snowfall) has trends that vary between −40 % to +40 % per century depending on the location.
On the other hand, projected trends in winter precipitation in the European Alps indicate mostly positive trends in 100-year return levels (Fig. 12 of Rajczak and Schär, 2017). For instance, an increase between 5 % and 30 % is expected under a high greenhouse gas emission scenario (comparing 2070–2099 to 1981–2010 for RCP8.5).

In and around the French Alps, studies analyzing extreme snowfall are rare (Beniston et al., 2018). A few papers describe the trends in extreme snowfall depending on elevation (Table 1). On one hand, past observations for Swiss stations below 1800 m present either a majority of decreasing trends or insignificant changes in the mean annual maximum of snowfall (Marty and Blanchet, 2012; Scherrer et al., 2013). On the other hand, climate projections for the Pyrenees and the Alps under a high greenhouse gas emission scenario (SRES A2 and RCP8.5, respectively) show that the mean seasonal maximum of snowfall is expected to decrease below a transition elevation and increase above it. López-Moreno et al. (2011) estimated a transition elevation of around 2000 m for the Pyrenees (comparing 2070–2100 to 1960–1990), while Frei et al. (2018) estimated a transition of around 3000 m for the Alps (comparing 2070–2099 to 1981–2010). In the French Alps, despite the existence of sufficient snowfall records and of previous studies that exploited them in an explicit extreme value framework, temporal trends in extreme snowfall remain poorly described. Indeed, earlier works focused rather on the spatial non-stationarity (with respect to latitude and longitude) of 3 d maxima of snowfall with max-stable processes (Davison et al., 2012). For instance, Gaume et al. (2013) estimated conditional 100-year return level maps at a fixed elevation of 2000 m, while Nicolet et al. (2016) found that the spatial dependence range of extreme snowfall has been decreasing.

Marty and Blanchet (2012)Scherrer et al. (2013)López-Moreno et al. (2011)de Vries et al. (2014)Frei et al. (2018)This study addresses the gap identified above, by assessing past temporal trends in the 23 massifs of the French Alps, with special emphasis on the 100-year return levels of daily snowfall. We rely on the SAFRAN reanalysis (Durand et al., 2009) available for the period 1959–2019, which provides, among other variables, time series of daily snowfall (from which annual maxima can be computed) for each massif and every 300 m of elevation between 600 and 3600 m (Vernay et al., 2019). In order to properly account for the specific statistical nature of maximal daily snowfall, our methodology relies on non-stationary extreme value models that depend on both elevation and time. Specifically, for each massif and four ranges of elevations (below 1000, 1000–2000, 2000–3000, and above 3000 m), temporal trends in 100-year return levels are estimated with a model selected on the basis of the Akaike information criterion.

We study annual maxima of daily snowfall in the French Alps, which are located between Lake Geneva to the north and the Mediterranean Sea to the south (Fig. 1).
This region is typically divided into 23 mountain massifs of about 1000 km^{2}, which correspond to 23 spatial units covering the French Alps (Vernay et al., 2019), the climate being considered homogeneous inside each massif for a given elevation.

The SAFRAN reanalysis (Durand et al., 2009; Vernay et al., 2019) combines large-scale reanalyses and forecasts with in situ meteorological observations to provide daily snowfall data, i.e., snow water equivalent of solid precipitation measured in kg m^{−2}, available for each massif from August 1958 to July 2019. We consider annual maxima of daily snowfall centered on the winter season; e.g., an annual maximum for the year 1959 corresponds to the maximum from 1 August 1958 to 31 July 1959.
Thus, we study annual maxima from 1959 to 2019.

The SAFRAN reanalysis focuses on the elevation dependency of meteorological conditions. Indeed, this reanalysis is not produced on a regular grid but provides data for each massif every 300 m of elevation. As illustrated in Fig. 1, we consider four ranges of elevation: below 1000 m, between 1000 and 2000 m, between 2000 and 3000 m, and above 3000 m. For instance, the maxima for the range “below 1000 m” correspond to the maxima at 600 m and the maxima at 900 m. We note that for each massif, we do not have any maxima above the top elevation of the massif.

The SAFRAN reanalysis has been evaluated both directly with in situ temperature and precipitation observations and indirectly with various snow depth observations compared to snow cover simulations of the model Crocus driven by SAFRAN atmospheric data (Durand et al., 2009; Vionnet et al., 2016; Quéno et al., 2016; Revuelto et al., 2018; Vionnet et al., 2019). Specifically, in Vionnet et al. (2019), the SAFRAN reanalysis has been evaluated for snowfall against two numerical weather prediction (NWP) systems for winter 2011–2012. The authors find that the seasonal snowfall averaged over all the massifs of the French Alps reaches 546 mm in SAFRAN, 684 mm in the first NWP, and 737 mm in the second NWP. In detail, they find that SAFRAN significantly differs from the two NWP systems in (i) areas of high elevation, probably due to the limited number of high-elevation stations and gauge undercatch, and (ii) areas on the windward side of the different mountain ranges due to the assumption of climatological homogeneity within each SAFRAN massif. In Ménégoz et al. (2020), the SAFRAN reanalysis has been compared to the regional climate model MAR which uses ERA-20C as forcing. The authors found that the vertical gradient of the annual mean of total precipitation of SAFRAN is generally smaller than those simulated by MAR.

## 3.1 Statistical distribution for annual maxima

Following the block maxima approach from extreme value theory (Coles, 2001),
we model annual maxima of daily snowfall with the generalized extreme value (GEV) distribution. Indeed theoretically, as the central limit theorem motivates asymptotically sample means modeling with the normal distribution, the Fisher–Tippett–Gnedenko theorem (Fisher and Tippett, 1928; Gnedenko, 1943) encourages asymptotically sample maxima modeling with the GEV distribution.
In practice, if *Y* is a random variable representing an annual maximum, we can assume that $Y\sim \text{GEV}(\mathit{\mu},\mathit{\sigma},\mathit{\xi})$.
Then, if *y* denotes an annual maximum,

where the three parameters are the location *μ*, the scale *σ*>0, and the shape *ξ*.
The GEV distribution encompasses three sub-families of distribution called reversed Weibull, Gumbel and Fréchet, which correspond to *ξ*<0, *ξ*=0, and *ξ*>0, respectively.

## 3.2 Elevational–temporal models

We consider non-stationary models that depend on both elevation and time.
Such models combine a stationary random component (a fixed extreme value distribution, e.g., GEV distribution) with non-stationary deterministic functions that map each covariate to the changing parameters of the distribution (Montanari and Koutsoyiannis, 2014).
Specifically, for each massif and each range of elevation (below 1000, 1000–2000, 2000–3000, and above 3000 m), if *Y*_{z,t} represents an annual maximum at the elevation *z* (within one of the four ranges of elevation) for the year *t* (between 1959 and 2019), we assume that

As illustrated in Table 2, we consider eight models that verify Eq. (2).
For a model ℳ, we denote as *θ*_{𝓜} the set of parameters of $\mathit{\mu}(z,t),\mathit{\sigma}(z,t),\text{and}\mathit{\xi}\left(z\right)$.
Following a preliminary analysis with pointwise distributions (Sect. 4.1), we consider models with a location and a scale parameter that vary linearly with respect to elevation.
Then, the shape parameter is either constant or linear with respect to elevation. Finally, the location and/or the scale can vary linearly with time.
As shown in Table 2,
we assume that the temporal and elevational effects are separable inside each range of elevation.
Thus, we do not consider models with cross terms, i.e., terms involving both the elevation and the years such as *z*×*t*. We discuss this assumption in Sect. 5.1.

First, models are fitted with the maximum likelihood method.
Let $\mathit{y}=({y}_{{z}_{\mathrm{1}},{t}_{\mathrm{1}}},\mathrm{\dots},{y}_{{z}_{\mathrm{1}},{t}_{M}},\mathrm{\dots},{y}_{{z}_{N},{t}_{\mathrm{1}}},\mathrm{\dots},{y}_{{z}_{N},{t}_{M}})$ represent a vector of annual maxima from year *t*_{1} to *t*_{M} and for the range of elevations containing ${z}_{\mathrm{1}},\mathrm{\dots},{z}_{N}$ for a given massif (Sect. 2).
We classically assume that maxima are conditionally independent given *θ*_{𝓜}.
For each model ℳ, we compute the maximum likelihood estimator ${\widehat{\mathit{\theta}}}_{\mathcal{M}}$ which corresponds to the parameter *θ*_{𝓜} that maximizes the likelihood *p*(** y**|

*θ*_{𝓜}), where $p\left(\mathit{y}\right|{\mathit{\theta}}_{\mathcal{M}})={\prod}_{z}{\prod}_{t}\frac{\partial P({Y}_{z,t}\le {y}_{z,t})}{\partial {y}_{z,t}}$.

Then, we select the model with the minimal Akaike information criterion (AIC) for each massif and range of elevations. Indeed, the AIC is the best criterion in a non-stationary context with small sample sizes (Kim et al., 2017). AIC equals $\mathrm{2}\times [\mathit{\#}{\mathit{\theta}}_{\mathcal{M}}-p(\mathit{y}\left|{\widehat{\mathit{\theta}}}_{\mathcal{M}}\right)]$,
where *#**θ*_{𝓜} is the number of parameters for the model ℳ.
Thus, minimizing the AIC corresponds to selecting models that both have few parameters, i.e., low *#**θ*_{𝓜}, and fit the data well, i.e., high $p\left(\mathit{y}\right|{\widehat{\mathit{\theta}}}_{\mathcal{M}})$.
Goodness of fit is assessed with Q–Q plots which show a good fit for the selected models (Appendix A).

## 3.3 Return levels

The *T*-year return level, which corresponds to a
quantile exceeded each year with probability $p=\frac{\mathrm{1}}{T}$, is the classical metric to quantify hazards of extreme events (Coles, 2001; Cooley, 2012).
We set $p=\frac{\mathrm{1}}{\mathrm{100}}=\mathrm{0.01}$ as it corresponds to the 100-year return period which is widely used for hazard mapping and the design of defense structure in France, notably for snow-related hazards (Eckert et al., 2010). Let ℳ denote a model from Table 2
and ${\widehat{\mathit{\theta}}}_{\mathcal{M}}$ the corresponding maximum likelihood estimator. Then, the associated return levels *y*_{p}, which depend on the elevation *z* and the year *t*, can be computed as follows:

We study trends in return levels.
For any considered model, the time derivative of the return level $\frac{\partial {y}_{p}(z,t)}{\partial t}$ is constant and quantifies the yearly change in return level.
Thus, for each range of elevations, a massif is said to have an increasing trend if the associated return level has increased, i.e., if $\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}>\mathrm{0}$. A massif has a decreasing trend if $\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}<\mathrm{0}$.
In the Result section, we display changes in 100-year return levels between 1959 and 2019, i.e., over the last 60 years, which equal $\mathrm{60}\times \frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}$.
If $\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}\ne \mathrm{0}$, i.e., if the selected model is temporally non-stationary, we compute the significance of the trend with a semi-parametric bootstrap resampling approach (Appendix B). We generate *B*=1000 bootstrap samples using the parameter ${\widehat{\mathit{\theta}}}_{\mathcal{M}}$. For each bootstrap sample *i*, we compute the time derivative of the return level ${\frac{\partial {y}_{p}(z,t)}{\partial t}}^{\left(i\right)}$.
Finally, a massif with an increasing trend is said to have a significant trend if
$\widehat{p}(\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}>\mathrm{0}|{\widehat{\mathit{\theta}}}_{\mathcal{M}})=\frac{\mathrm{1}}{B}{\sum}_{i=\mathrm{1}}^{B}{\mathbf{1}}_{{\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}}^{\left(i\right)}>\mathrm{0}}>\mathrm{1}-\mathit{\alpha}$, where *α*=5 % is the significance level.
In other words, the trend is significantly increasing when the percentage of bootstrap samples for which the return levels are increasing is above the threshold 1−*α*.
Likewise, a massif has a significant decreasing trend if $\widehat{p}(\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}<\mathrm{0}|{\widehat{\mathit{\theta}}}_{\mathcal{M}})>\mathrm{1}-\mathit{\alpha}$.

## 3.4 Workflow

In Sect. 4.1, we analyze changes in pointwise distribution of annual snowfall maxima with elevation in the 23 massifs of the French Alps, which helped us define the eight elevational–temporal models considered (Sect. 3.2). Pointwise distribution stands for a distribution fitted on the annual maxima from a single elevation of one massif. Specifically, we fit a pointwise GEV distribution with the maximum likelihood method for each massif every 300 m of elevation from 600 to 3600 m. We exclude physically implausible distributions; i.e., $\mathit{\xi}\notin \left[-\mathrm{0.5},\mathrm{0.5}\right]$ (Martins and Stedinger, 2000). Then, we compute elevation gradients for the three GEV parameters and the 100-year return level with a linear regression.

In Sect. 4.2, we compare pointwise distributions with our approach based on piecewise elevational–temporal models.
Piecewise models stand for models fitted on the annual maxima from a range of elevation of one massif.
We present the elevational–temporal models selected in each massif for each range of elevations (below 1000, 1000–2000, 2000–3000, and above 3000 m) obtained with the methodology described in Sect. 3.2.
First, we fit the eight models from Table 2 with the maximum likelihood method. Then, we select one model with the AIC. Finally,
if the selected model is temporally non-stationary, we assess the significance of the trend using a semi-parametric bootstrap resampling approach with a significance level *α*=5 % (Sect. 3.3).

In Sect. 4.3, we present for each massif and range of elevations the temporal trends in 100-year return levels obtained from selected models. We compute 100-year return levels in 2019 and their changes between 1959 and 2019 with the selected elevational–temporal models (Sect. 3.3). For each range of elevations, a massif has an increasing (decreasing) trend if the associated return level has increased (decreased).

## 4.1 Pointwise distribution for each elevation.

According to pointwise fits, the location and scale parameters increase linearly with elevation (Fig. 2a and b). *R*^{2} coefficients are always larger than 0.8, except for the Bauges massif for the scale parameter (Fig. 3b).
The average elevation gradient for the location and the scale parameters is equal to
$\mathrm{2.1}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}/\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$ and $\mathrm{0.39}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}/\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$, respectively (Fig. 3a and b).
In particular, this linear augmentation is also valid for any range of elevations considered to fit elevational–temporal models.
Thus, as explained in Sect. 3.2, we assume for elevational–temporal models location and scale parameters that vary linearly with respect to elevation.
On the other hand, changes in the shape parameter rarely follow a linear relationship with elevation between 600 and 3600 m. Indeed, only seven massifs have *R*^{2} coefficients larger than 0.5 (Fig. 3c).
However, as illustrated in Fig. 2c, it does not preclude the shape parameter from varying linearly with the elevation
locally, i.e., within an elevation range.
Therefore, as explained in Sect. 3.2, we assume for elevational–temporal models that the shape parameter is either constant or linear with respect to elevation.

In Fig. 2d we show the change in 100-year return levels with elevation, while in Fig. 3d we display their elevation gradients.
Return levels augment linearly with elevation, which is confirmed by *R*^{2} coefficients always larger than 0.8.
The largest return levels and elevation gradients correspond to the Mercantour (Southern Alps) and Haute-Maurienne massif (eastern part of the French Alps).

## 4.2 Elevational–temporal models for each range of elevations

Figure 4 illustrates selected models for each massif and each range of elevations.
The most selected model is a temporally non-stationary models ${\mathcal{M}}_{{\mathit{\mu}}_{t},{\mathit{\sigma}}_{t}}$, which is selected for 54 % of the massifs and is significant for 32 % of the massifs. Further, the temporally stationary model ℳ_{0} and ${\mathcal{M}}_{{\mathit{\xi}}_{z}}$ have been selected for 28 % and 1 % of the massifs.
The remaining temporally non-stationary models have been selected for 17 % of massifs.
We notice that the four-most-represented temporally non-stationary models have a linearity with respect to the years *t* for the scale parameter (which impacts both the mean and the variance of the GEV distribution), potentially indicating that often both the intensity and the variance of maxima are changing over time.

Furthermore, we observe that the shape parameter values remain between −0.4 and 0.4, which is a physically acceptable range (Martins and Stedinger, 2000). Except a majority of Weibull distributions (massifs displayed in green in Fig. 4) at 500 m in the northwest of the French Alps and a clear majority of Fréchet distributions (yellow massifs) at 2500 m, we do not find any clear spatial/elevation patterns for the shape parameter.

Figure 5 exemplifies the differences between pointwise distribution and our approach based on piecewise elevational–temporal models. We consider annual maxima from 600 to 3600 m of elevation for the Vanoise massif (Sect. 2). First, our approach makes it possible to interpolate GEV parameter values (and thus to deduce 100-year return levels) for each range of elevations (blue line) rather than having point estimates (green dots). Furthermore, it reduces confidence intervals for return levels (shaded areas) which were computed with an approach based on semi-parametric bootstrap resampling (Appendix B). Finally, our approach accounts for temporal trends. For example, for the Vanoise massif above 3000 m, the selected model is a temporally non-stationary model (Fig. 4) with an increasing trend in return levels (Fig. 8). This explains why return levels in 2019 estimated from the elevational–temporal model exceed return levels estimated from the pointwise distribution.

## 4.3 Temporal trends in return levels

Figure 6 shows that both increasing and decreasing trends in 100-year return levels are found for all elevation ranges. We also observe that a majority of trends are decreasing below 2000 m and increasing above 2000 m. If we analyze only significant trends, the elevation pattern remains the same. On one hand, more than 30 % of massifs are significantly decreasing below 2000 m: 40 % below 1000 m and 30 % for the range 1000–2000 m. On the other hand, roughly 30 % of massifs are significantly increasing above 2000 m: slightly less than 30 % for the range 2000–3000 m and slightly less than 40 % above 3000 m. We note that the sign and the significance of the trends (summarized with the percentages in Fig. 6) remain largely similar for the trends in events of 10- and 50-year return periods (Appendix C).

In Fig. 7, we show distributions of changes and relative changes in 100-year return levels between 1959 and 2019. We find a temporal increase in 100-year return levels between 1959 and 2019 equal to +23 % ($+\mathrm{32}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$) on average at 3500 m and a decrease of −10 % ($-\mathrm{7}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$) on average at 500 m. For intermediate elevations, i.e., between 1000 and 3000 m, we observe that the distributions of changes and of relative changes remain roughly negative (decrease) at 1500 m and positive (increase) at 2500 m. This result holds for all massifs, for the subset of massifs with a selected model temporally non-stationary, and for the subset with a selected model that is temporally non-stationary and significant.

In Fig. 8, we display the change in 100-year return levels between 1959 and 2019 for each range of elevations. At 500 m, we observe that eight massifs have a stationary trend and five massifs located in the center of the French Alps have a significant decreasing trend. We also note that two massifs located in the western French Alps have an increasing significant trend, with an absolute change in the 100-year return level close to $+\mathrm{20}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$. At 1500 m, six massifs in the center of the French Alps have decreasing trends. At 2500 m, we observe a spatially contrasting pattern: most decreasing trends are located in the north, while most increasing trends are located in the south. We discuss this pattern in Sect. 5.4. At 3000 m, we observe that the six massifs with a significant increasing trend are located in the south of the French Alps.

In Fig. 9, we display the return level in 2019 for each range of elevations.
Combining Figs. 8 and 9 enables us to pinpoint massifs both with high return levels in 2019 and with an increasing or decreasing trend.
For instance, at 2500 m, the Haute-Maurienne massif has one of the highest 100-year return levels in 2019 (185 kg m^{−2}), but it is decreasing with time. On the other hand, at 2500 and 3500 m, most massifs in the south have high 100-year return levels values in 2019 and increasing trends.

## 5.1 Methodological considerations

We discuss the statistical models considered to estimate temporal trends in 100-year return levels of daily snowfall.

For small-size time series of annual maxima, e.g., a few decades, return level uncertainty largely depends on the uncertainty in the shape parameter (Koutsoyiannis, 2004). In order to reduce the uncertainty in the GEV parameters, models are often fitted using the information from several time series e.g., using regionalization methods (Hosking et al., 2009) or using scaling relationships between different aggregation durations (Blanchet et al., 2016). In this work, we fit models to time series from several elevations. In the literature, elevation is often not treated as a covariate but rather accounted for by a spatial distance (Blanchet and Davison, 2011; Gaume et al., 2013; Nicolet et al., 2016). In practice, as illustrated in Fig. 5, we compute the uncertainties in return levels for all the massifs and find that elevational–temporal models effectively decrease confidence intervals compared to pointwise distributions fitted to one time series.

Then, we fit the models to at least two time series, i.e., from at least two elevations. As illustrated in Fig. 1, for each range of elevations, we always have at least two time series, i.e., more than 100 maxima to estimate 100-year return levels. However, in practice annual maxima from consecutive elevations are often dependent. At low elevations, four time series contain zeros, i.e., years without any snowfall, which may lead to misestimation for the models. A potential solution would be to rely on a mixed discrete–continuous distribution: a discrete distribution for the probability of a year without any snowfall and a continuous GEV distribution for the annual maxima of snowfall. However, this would require at least one additional parameter for the discrete distribution and even more if we wish to model some non-stationarity. In practice, to avoid overparameterized models, we removed one time series which had more than 10 % of zeros.

Afterwards, for different ranges of elevation containing at most three consecutive elevations, we fit the models using all corresponding time series.
For each model, this ensures that the temporal non-stationarity can be assumed to not depend on the elevation.
Indeed, initially we intended for each massif to fit a single model to time series from all elevations.
However, to account for decreasing trends at low elevations and increasing trends at high elevations, this led to complex overparameterized models that often did not fit well.
We decided to consider a piecewise approach, i.e., simpler models fitted to
ranges of consecutive elevations at most separated by 900 m
(Figs. 1 and 5).
This ensures that we can assume that the temporal non-stationarity (no trend or decreasing/increasing trend) is shared between all elevations from the same range.
Otherwise, we observe that annual maxima from consecutive altitudes are dependent (Fig. 1). We did not account for this dependence in our statistical model. Instead, we assumed that maxima are conditionally independent given the vector of parameters *θ*_{𝓜} (Sect. 3.2).

Finally, for each range of elevations, we consider models with a temporal non-stationarity only for the location and scale parameter. Indeed, in the literature, a linear non-stationarity is considered sometimes only for the location parameter (Fowler et al., 2010; Tramblay and Somot, 2018) but more often for both the location and the scale (or log-transformed scale for numerical reasons) parameters (Katz et al., 2002; Kharin and Zwiers, 2004; Marty and Blanchet, 2012; Wilcox et al., 2018). Otherwise, the shape parameter is typically considered temporally stationary in the literature, and we followed this approach.

## 5.2 Implication of the temporal trends in 100-year return levels

In Fig. 8, we emphasize massifs with a strong increase in 100-year return levels between 1959 and 2019, e.g., massifs filled with medium/dark red correspond to increase $\ge +\mathrm{50}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$. Settlements in these massifs should ensure that the design of protective measures and building standards against extreme snow events are still adequate after such an increase. Hopefully, this might concern few settlements as such strong increases are always located above 2000 m. However, this might impact ski resorts which should ensure that the design of avalanche protection measures takes this change into account. We note that extreme snow events can sometimes be triggered by one snowfall event but often depend on other factors such as accumulated snow or wind. Therefore, to update structure standards for ground snow load (Biétry, 2005), we should account for both this increasing trend in annual maxima of daily snowfall and trends in annual maxima of ground snow loads (Le Roux et al., 2020). Indeed, most known snow load destructions result from such intense and short snow events, sometimes combined with liquid precipitation, which is not considered in this study. In general, in mountainous regions around the French Alps, if the past trends continue into the future, for extreme snowfall we can expect decreasing trends below 1000 m and increasing trends above 3000 m, as this agrees with both our results and the literature (Table 1).

## 5.3 Data considerations

Following the evaluation of the SAFRAN reanalysis cited in Sect. 2, we can conclude that this reanalysis most likely underestimates high-elevation precipitation (above 2000 m), which probably leads to an underestimation of high-elevation snowfall. This deficiency does not affect the main contribution of this article, i.e., a majority of decreasing trends below 2000 m and a majority of increasing trends above 2000 m. However, this deficiency affects the value of extreme snowfall, i.e., the 100-year return level and the scale of its changes, which may be underestimated above 2000 m. For future works, we note that a direct evaluation of extreme snowfall could help to better pinpoint the locations where return levels might be biased.

## 5.4 Hypothesis for the contrasting pattern for changes in 100-year return levels

In Fig. 8, at 2500 m, we find a spatially contrasting pattern for changes in 100-year return levels of snowfall: most decreasing trends are located in the north, while most increasing trends are in the south.

These changes contradict expectations based on the climatological differences between the north and south of the French Alps.
Indeed, since the north is climatologically colder than the south in both winter and summer (Fig. 5 of Durand et al., 2009), we would have expected the inverse pattern for an intermediate elevation, i.e., increasing trends in the north and decreasing trends in the south.
Indeed, extreme snowfall stems from extreme precipitation occurring in a range of optimal temperatures slightly below 0 ^{∘}C (O'Gorman, 2014; Frei et al., 2018).
Thus, under global warming, we would have expected for an intermediate elevation
that the probability of being in the range of optimal temperatures would be increasing in the north (because mean temperatures would shift toward the optimal range)
while decreasing in the south (because temperatures would be shifting away from the optimal range).
Therefore, we would have observed an increase in extreme snowfall in the north and a decrease in the south.

Thus, this spatial pattern of changes cannot be solely explained by the spatial pattern of mean temperature. In particular, we believe that dynamical changes, i.e., heterogeneous changes in extreme precipitation in the French Alps, may have contributed to generate this contrasting pattern. In Appendix D, we apply our methodology (Sect. 3) on daily winter precipitation from the SAFRAN reanalysis for the period 1959–2019. We focus on winter precipitation because winter is the season when most annual maxima of daily snowfall occur below 3000 m (Appendix E). We observe that at 2500 m (and at all elevations) changes in 100-year return levels of winter precipitation show the same contrasting pattern as 100-year return levels of snowfall

This contrasting pattern is observed at all elevations for changes in 100-year return levels of winter precipitation. The literature also confirms that changes in extreme precipitation are not homogeneous. For instance, we observe for the period 1903–2010 that trends in daily maxima of winter precipitation are stronger in the south (+20 %–40 % per century) compared to the north (from −10 % to +20 % per century) of the French Alps (Fig. 7 of Ménégoz et al., 2020). We also observe for the period 1958–2017 that the 20-year return level of winter precipitation has decreased in the north of the French Alps and has slightly increased or remained the same in the south (Fig. 8 of Blanchet et al., 2021). This observation might be due to a stronger increasing trend in extreme precipitation for the Mediterranean circulation than for the Atlantic circulation. Indeed, precipitation maxima in the north of the French Alps are frequently triggered by the Atlantic circulation, while maxima in the south are often due to the Mediterranean circulation (Blanchet et al., 2020). Furthermore, increasing trends in extreme snowfall have already been observed in the proximity of the Mediterranean Sea. For example, Faranda (2020) identified a certain number of Mediterranean countries showing positive changes in snowfall maxima. D'Errico et al. (2020) propose a physical explanation of this phenomenon: the Mediterranean Sea is warming faster than any other ocean, which enhances convective precipitation and favors heavy snowfalls during cold-spell events.

In practice, this increasing trend in extreme snowfall in the south should be temporary. Indeed, with climate change, temperatures are expected to shift further away from the optimal range of temperatures for extreme snowfall. Thus, in the long run, extreme snowfall is expected to decrease as the increase in extreme precipitation shall not compensate for the decreasing probability of being close to the optimal range of temperatures.

We estimate temporal trends in 100-year return levels of daily snowfall for several ranges of elevation based on the SAFRAN reanalysis available from 1959 to 2019 (Durand et al., 2009). Our statistical methodology relies on non-stationary extreme value models that depend on both elevation and time. Our results show that a majority of trends are decreasing below 2000 m and increasing above 2000 m. Quantitatively, we find an increase in 100-year return levels between 1959 and 2019 equal to +23 % ($+\mathrm{32}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$) on average at 3500 m and a decrease of −10 % ($-\mathrm{7}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{m}}^{-\mathrm{2}}$) on average at 500 m. For the four investigated elevation ranges, we find both decreasing and increasing trends depending on location. In particular, we observe a spatially contrasting pattern, exemplified at 2500 m: 100-year return levels have decreased in the north of the French Alps while they have increased in the south. In the discussion, we highlight that this pattern might be related to known increasing trends in extreme snowfall in the proximity of the Mediterranean Sea.

Many potential extensions of this work could be considered.
First, reanalyses are increasingly available at the European scale (e.g., Soci et al., 2016), which could be used for extending this work to a wider geographical scale.
In this case, instead of considering close massifs as spatially independent, we believe that our methodology may benefit from an explicit modeling of the spatial dependence (Padoan et al., 2010).
Then, climatic projections could enable us to explore temporal trends up to the end of the twenty-first century.
In these circumstances, it might be more relevant to use global mean surface temperature as a temporal covariate to combine ensembles of climate models.
Finally, future research should focus not solely on mountain regions but also on lowland regions such as around the Mediterranean Sea. Indeed, such regions are often heavily impacted by snow-related hazards because they are ill-equipped for such rare events. For instance, extreme snowfall over Roussillon, a Mediterranean coastal lowland, caused major damage in 1986 (Vigneau, 1987), while in 2021 heavy snowfall over Spain caused at least EUR 1.4 billion of damage (The New York Times, 2021).
In these regions, temperatures below the rain–snow transition temperature, i.e., roughly below 0 ^{∘}C, may tend to be rare in the future.
Therefore, in these cases, in addition to directly studying trends in snowfall extremes, we should focus on trends in the compound risk of cold wet events (De Luca et al., 2020).

A quantile–quantile (Q–Q) plot is a standard diagnosis
tool
based on the comparison of empirical quantiles (computed from the empirical distribution) and theoretical quantiles (computed from the expected distribution).
For non-stationary extreme value models, the approach is two-fold (Coles, 2001; Katz, 2012).
First, we transform observations into residuals with
a probability integral transformation *f*_{GEV→Standard Gumbel}.
Then, we construct a Q–Q plot to assess if the residuals follow a standard Gumbel distribution. If the Q–Q plot reveals a good fit, it means that the non-stationary extreme value model has a good fit as well.

We start by transforming the observations ${y}_{{z}_{\mathrm{1}},{t}_{\mathrm{1}}},\mathrm{\dots},{y}_{{z}_{\mathrm{1}},{t}_{M}},\mathrm{\dots},{y}_{{z}_{N},{t}_{\mathrm{1}}}\mathrm{\dots},{y}_{{z}_{N},{t}_{M}}$ into residuals. Let ${Y}_{z,t}\sim \text{GEV}\left(\mathit{\mu}\right(z,t),\mathit{\sigma}(z,t),\mathit{\xi}(z\left)\right)$ with parameter *θ*_{𝓜}. By definition of the probability integral transformation *f*_{GEV→Standard Gumbel}, we obtain that

The transformed observations, a.k.a. residuals, are denoted as ${\mathit{\u03f5}}_{z,t}={f}_{\text{GEV}\to \text{Standard Gumbel}}({y}_{z,t};{\mathit{\theta}}_{\mathcal{M}})$.
Afterwards, we construct a Q–Q plot to assess if the residuals follow a standard Gumbel distribution.
On one hand, the *N*×*M* empirical quantiles correspond to the ordered values of the residuals ${\mathit{\u03f5}}_{{z}_{\mathrm{1}},{t}_{\mathrm{1}}},\mathrm{\dots},{\mathit{\u03f5}}_{{z}_{\mathrm{1}},{t}_{M}},\mathrm{\dots},{\mathit{\u03f5}}_{{z}_{N},{t}_{\mathrm{1}}}\mathrm{\dots},{\mathit{\u03f5}}_{{z}_{N},{t}_{M}}$.
On the other hand, we compute the corresponding *N*×*M* theoretical quantiles, which are the quantile $\frac{i}{N\times M+\mathrm{1}}$ of the standard Gumbel distribution, where $i\in \mathit{\{}\mathrm{1},\mathrm{\dots},N\times M\mathit{\}}$.

In Fig. A1, we display Q–Q plots for the selected model for the time series displayed in Fig. 1. We observe that they show a good fit, as the points stay close to the line. In general, most retained models show a good fit. Furthermore quantitatively, if we rely on an Anderson–Darling statistical test with a 5 % significance level to assess if the residuals follow a standard Gumbel distribution, we find that the largest parts of the tests are not rejected (not shown).

In the context of maximum likelihood estimation, uncertainty related to return levels can be evaluated with the delta method, which quickly provides confidence intervals in both the stationary and the non-stationary cases (Coles, 2001; Gilleland and Katz, 2016). However, due to the dependence between maxima from consecutive elevations (Fig. 1), we decided to compute confidence intervals with a bootstrap resampling method (Efron and Tibshirani, 1993). This resampling method allows us to estimate the uncertainties resulting from in-sample variability. In this article, we rely on a semi-parametric bootstrap resampling method adapted to non-stationary extreme models (Kharin and Zwiers, 2004; Sillmann et al., 2011).

We generate *B*=1000 bootstrap samples using the parameter ${\widehat{\mathit{\theta}}}_{\mathcal{M}}$. For each bootstrap sample *i*, the semi-parametric bootstrap method is four-fold. First, as explained in Appendix A, we compute the residuals ${\mathit{\u03f5}}_{{z}_{\mathrm{1}},{t}_{\mathrm{1}}},\mathrm{\dots},{\mathit{\u03f5}}_{{z}_{\mathrm{1}},{t}_{M}},\mathrm{\dots},{\mathit{\u03f5}}_{{z}_{N},{t}_{\mathrm{1}}}\mathrm{\dots},{\mathit{\u03f5}}_{{z}_{N},{t}_{M}}$.
Then, from these residuals we draw with replacement a sample of size *M*×*N*. We denote these bootstrapped residuals as ${\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}}_{{z}_{\mathrm{1}},{t}_{\mathrm{1}}},\mathrm{\dots},{\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}}_{{z}_{\mathrm{1}},{t}_{M}},\mathrm{\dots},{\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}}_{{z}_{N},{t}_{\mathrm{1}}}\mathrm{\dots},{\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}}_{{z}_{N},{t}_{M}}$.
Afterwards, we transform these bootstrapped residuals into bootstrapped annual maxima as follows:
${\stackrel{\mathrm{\u0303}}{y}}_{z,t}={f}_{\text{GEV}\to \text{Standard Gumbel}}^{-\mathrm{1}}({\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}}_{z,t};{\mathit{\theta}}_{\mathcal{M}})$.
Finally, we estimate the parameter ${\widehat{\mathit{\theta}}}_{\mathcal{M}}^{\mathbf{\left(}\mathit{i}\mathbf{\right)}}$ of model ℳ with the bootstrapped annual maxima
${\stackrel{\mathrm{\u0303}}{y}}_{{z}_{\mathrm{1}},{t}_{\mathrm{1}}},\mathrm{\dots},{\stackrel{\mathrm{\u0303}}{y}}_{{z}_{\mathrm{1}},{t}_{M}},\mathrm{\dots},{\stackrel{\mathrm{\u0303}}{y}}_{{z}_{N},{t}_{\mathrm{1}}}\mathrm{\dots},{\stackrel{\mathrm{\u0303}}{y}}_{{z}_{N},{t}_{M}}$.
To sum up, this bootstrap procedure provides a set
$\mathit{\{}{\widehat{\mathit{\theta}}}_{\mathcal{M}}^{\mathbf{\left(}\mathbf{1}\mathbf{\right)}},\mathrm{\dots},{\widehat{\mathit{\theta}}}_{\mathcal{M}}^{\mathbf{\left(}\mathit{i}\mathbf{\right)}},\mathrm{\dots},{\widehat{\mathit{\theta}}}_{\mathcal{M}}^{\mathbf{\left(}\mathit{B}\mathbf{\right)}}\mathit{\}}$
of *B* parameters for the model ℳ.

In practice, we rely on this set of GEV parameters to obtain 80 % confidence intervals for 100-year return levels (Fig. 5) or for time derivatives of 100-year return levels (Sect. 3.3). For instance, in the latter case we have $\widehat{p}(\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}>\mathrm{0}|{\widehat{\mathit{\theta}}}_{\mathcal{M}})=\frac{\mathrm{1}}{B}{\sum}_{i=\mathrm{1}}^{B}{\mathbf{1}}_{{\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}}^{\left(i\right)}>\mathrm{0}}$, where ${\frac{\partial {y}_{\mathrm{0.01}}(z,t)}{\partial t}}^{\left(i\right)}$ is the time derivative of 100-year return levels for the parameter ${\widehat{\mathit{\theta}}}_{\mathcal{M}}^{\mathbf{\left(}\mathit{i}\mathbf{\right)}}$.

The 100-year return period was chosen because it is the largest return period considered in the Eurocodes for building structures (Cabrera et al., 2012). We believe that this return period is the most familiar return period for non-experts as it corresponds to a centennial event. For smaller return periods (5–10 years), our results also apply. In Fig. C1, we illustrate our results for the 10-, 50-, and 100-year return periods. We observe that the overall distribution of increasing/decreasing trend for the return levels is almost insensitive to the choice of the return period. For instance, the only noticeable difference between the 10- and 100-year return periods is that for the elevation range 1000–2000 m and for the elevation range 2000–3000 m, we observe that one massif shows an increasing trend for the return period of 10 years, while it is decreasing for the return period of 100 years.

We apply the same methodology as in our study (Sect. 3) to daily winter (December to February) precipitation obtained with the SAFRAN reanalysis and spanning the period 1959–2019. First, a preliminary analysis with pointwise fits indicates that a linear parameterization with respect to the elevation for the location and scale parameters is also valid for the winter precipitation (Fig. D1). In Fig. D2, we illustrate changes in the 100-year return level of winter precipitation. We observe a spatially contrasting pattern at all elevations, i.e., increase in the south and decrease in the north. This underlines that the spatially contrasting pattern observed for changes in the 100-year return level of snowfall at 2500 m may result from the circulation patterns of precipitation.

In Fig. E1, we study the seasons when the annual maxima of daily snowfall occurred. For elevation range 1 (below 1000 m) and for elevation range 2 (between 1000 and 2000 m), we observe that the annual maxima mainly occurred (>60 %) between December and February, i.e., the coldest part of the snow season. For elevation range 3 (between 2000 and 3000 m) more than 40 % of maxima occurred in winter, while slightly less than 30 % occurred in autumn and in spring. For elevation range 4 (above 3000 m) the seasons of occurrence are more spread, even if we observe that more than 40 % of maxima occurred in autumn. In conclusion, we find that below 3000 m, most annual maxima of daily snowfall occur in winter, while above 3000 m they mostly occur in autumn.

The full SAFRAN reanalysis on which this study is grounded is freely available on AERIS (Vernay et al., 2019, https://doi.org/10.25326/37).

ELR, GE, and NE designed the research. ELR performed the analysis and drafted the first version of the manuscript. All authors discussed the results and edited the manuscript.

The authors declare that they have no conflict of interest.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We are grateful to Eric Gilleland for his “extRemes” R package and Ali Saeb for his “gnFit” R package. INRAE, CNRM, and IGE are members of Labex OSUG.

Erwan Le Roux holds a PhD grant from INRAE.

This paper was edited by Jürg Schweizer and reviewed by J. Ignacio López-Moreno, Sven Kotlarski, and one anonymous referee.

Beniston, M., Farinotti, D., Stoffel, M., Andreassen, L. M., Coppola, E., Eckert, N., Fantini, A., Giacona, F., Hauck, C., Huss, M., Huwald, H., Lehning, M., López-Moreno, J.-I., Magnusson, J., Marty, C., Morán-Tejéda, E., Morin, S., Naaim, M., Provenzale, A., Rabatel, A., Six, D., Stötter, J., Strasser, U., Terzago, S., and Vincent, C.: The European mountain cryosphere: a review of its current state, trends, and future challenges, The Cryosphere, 12, 759–794, https://doi.org/10.5194/tc-12-759-2018, 2018. a

Biétry, J.: Charges de neige au sol en France: proposition de carte révisée, Tech. Rep. 0, groupe de travail “Neige” de la commission de normalisation BNTEC P06A, 2005. a

Blanchet, J. and Davison, A. C.: Spatial modeling of extreme snow depth, Ann. Appl. Stat., 5, 1699–1725, https://doi.org/10.1214/11-AOAS464, 2011. a

Blanchet, J., Marty, C., and Lehning, M.: Extreme value statistics of snowfall in the Swiss Alpine region, Water Resour. Res., 45, 5, https://doi.org/10.1029/2009WR007916, 2009. a

Blanchet, J., Ceresetti, D., Molinié, G., and Creutin, J. D.: A regional GEV scale-invariant framework for Intensity–Duration–Frequency analysis, J. Hydrol., 540, 82–95, https://doi.org/10.1016/j.jhydrol.2016.06.007, 2016. a

Blanchet, J., Creutin, J.-D., and Blanc, A.: Retreating Winter and Strengthening Autumn Mediterranean Influence on Extreme Precipitation in the Southwestern Alps over the last 60 years, Environ. Res. Lett., 16, https://doi.org/10.1088/1748-9326/abb5cd, 2020. a

Blanchet, J., Blanc, A., and Creutin, J.-D.: Explaining recent trends in extreme precipitation in the Southwestern Alps by changes in atmospheric influences, Weather and Climate Extremes, 100356, https://doi.org/10.1016/j.wace.2021.100356, 2021. a

Cabrera, A. T., Heras, M. D., Cabrera, C., and Heras, A. M. D.: The Time Variable in the Calculation of Building Structures. How to extend the working life until the 100 years?, 2nd International Conference on Construction and Building Research, Universitat Politècnica de València Valencia, Spain, 14–16 November 2012, pp. 1–6, http://oa.upm.es/22914/1/INVE_MEM_2012_152534.pdf, 2012. a

Changnon, S. A.: Catastrophic winter storms: An escalating problem, Climatic Change, 84, 131–139, https://doi.org/10.1007/s10584-007-9289-5, 2007. a

Coles, S. G.: An introduction to Statistical Modeling of Extreme Values, vol. 208, Springer, London, https://doi.org/10.1007/978-1-4471-3675-0, 2001. a, b, c, d

Cooley, D.: Return Periods and Return Levels Under Climate Change, in: Extremes in a Changing Climate – Detection, Analysis and Uncertainty, pp. 97–114, Springer Science and Business Media, available at: https://www.springer.com/gp/book/9789400744783 (last access: 7 September 2021), 2012. a

Davison, A. C., Padoan, S. A., and Ribatet, M.: Statistical Modeling of Spatial Extremes, Stat. Sci., 27, 161–186, https://doi.org/10.1214/11-STS376, 2012. a

De Luca, P., Messori, G., Faranda, D., Ward, P. J., and Coumou, D.: Compound warm–dry and cold–wet events over the Mediterranean, Earth Syst. Dynam., 11, 793–805, https://doi.org/10.5194/esd-11-793-2020, 2020. a

de Vries, H., Lenderink, G., and van Meijgaard, E.: Future snowfall in western and central Europe projected with a high-resolution regional climate model ensemble, Geophys. Res. Lett., 41, 4294–4299, https://doi.org/10.1002/2014GL059724, 2014. a

D'Errico, M., Yiou, P., Nardini, C., Lunkeit, F., and Faranda, D.: A dynamical and thermodynamic mechanism to explain heavy snowfalls in current and future climate over Italy during cold spells, Earth Syst. Dynam. Discuss. [preprint], https://doi.org/10.5194/esd-2020-61, in review, 2020. a

Durand, Y., Laternser, M., Giraud, G., Etchevers, P., Lesaffre, B., and Mérindol, L.: Reanalysis of 44 year of climate in the French Alps (1958-2002): Methodology, model validation, climatology, and trends for air temperature and precipitation, J. Appl. Meteorol. Clim., 48, 429–449, https://doi.org/10.1175/2008JAMC1808.1, 2009. a, b, c, d, e, f

Eckert, N., Naaim, M., and Parent, E.: Long-term avalanche hazard assessment with a Bayesian depth-averaged propagation model, J. Glaciol., 56, 563–586, https://doi.org/10.3189/002214310793146331, 2010. a

Efron, B. and Tibshirani, R. J.: An introduction to the bootstrap, Chapman and Hall, London, New York, available at: http://www.ru.ac.bd/stat/wp-content/uploads/sites/25/2019/03/501_02_Efron_Introduction-to-the-Bootstrap.pdf (last access: 7 September 2021), 1993. a

Faranda, D.: An attempt to explain recent changes in European snowfall extremes, Weather Clim. Dynam., 1, 445–458, https://doi.org/10.5194/wcd-1-445-2020, 2020. a

Fisher, R. A. and Tippett, L. H. C.: Limiting forms of the frequency distribution of the largest or smallest member of a sample, Math. Proc. Cambridge, 24, 180–190, https://doi.org/10.1017/S0305004100015681, 1928. a

Fowler, H. J., Cooley, D., Sain, S. R., and Thurston, M.: Detecting change in UK extreme precipitation using results from the climateprediction.net BBC climate change experiment, Extremes, 13, 241–267, https://doi.org/10.1007/s10687-010-0101-y, 2010. a

Frei, P., Kotlarski, S., Liniger, M. A., and Schär, C.: Future snowfall in the Alps: projections based on the EURO-CORDEX regional climate models, The Cryosphere, 12, 1–24, https://doi.org/10.5194/tc-12-1-2018, 2018. a, b, c, d, e

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

Gilleland, E. and Katz, R. W.: extRemes 2.0: An Extreme Value Analysis Package in R, J. Stat. Softw., 72, 1–39, https://doi.org/10.18637/jss.v072.i08, 2016. a

Gnedenko, B.: Sur la distribution limite du terme maximum d'une série aléatoire, Ann. Math., 44, 423–453, https://doi.org/10.2307/1968974, 1943. a

Hosking, J. R. M., Wallis, J. R., Hosking, J. R. M., and Wallis, J. R.: Regional frequency analysis, Cambridge University Press, Cambridge, https://doi.org/10.1017/cbo9780511529443.003, 2009. a

IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK and New York, NY, USA, 1535 pp, available at: https://www.researchgate.net/profile/Abha_Chhabra2/publication/271702872_Carbon_and_Other_Biogeochemical_Cycles/links/54cf9ce80cf24601c094a45e/Carbon-and-Other-Biogeochemical-Cycles.pdf (last access: 7 September 2021), 2013. a

IPCC: Special Report: Special report on the ocean and cryosphere in a changing climate, available at: https://www.ipcc.ch/report/srocc/ (last access: 7 September 2021), 2019. a

Katz, R. W.: Statistical methods for nonstationary extremes, in: Extremes in a Changing Climate – Detection, Analysis and Uncertainty, pp. 15–38, Springer Science and Business Media, London, https://doi.org/10.1007/978-94-007-4479-0, 2012. a

Katz, R. W., Parlange, M. B., and Naveau, P.: Statistics of extremes in hydrology, Adv. Water Resour., 25, 1287–1304, https://doi.org/10.1016/S0309-1708(02)00056-8, 2002. a

Kharin, V. V. and Zwiers, F. W.: Estimating extremes in transient climate change simulations, J. Climate, 18, 1156–1173, https://doi.org/10.1175/JCLI3320.1, 2004. a, b

Kim, H., Kim, S., Shin, H., and Heo, J. H.: Appropriate model selection methods for nonstationary generalized extreme value models, J. Hydrol., 547, 557–574, https://doi.org/10.1016/j.jhydrol.2017.02.005, 2017. a

Koutsoyiannis, D.: Statistics of extremes and estimation of extreme rainfall: I. Theoretical investigation, Hydrolog. Sci. J., 49, 575–590, https://doi.org/10.1623/hysj.49.4.575.54430, 2004. a

Le Roux, E., Evin, G., Eckert, N., Blanchet, J., and Morin, S.: 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, https://doi.org/10.5194/nhess-20-2961-2020, 2020. a

López-Moreno, J. I., Goyette, S., Vicente-Serrano, S. M., and Beniston, M.: Effects of climate change on the intensity and frequency of heavy snowfall events in the Pyrenees, Climatic Change, 105, 489–508, https://doi.org/10.1007/s10584-010-9889-3, 2011. a, b

Martins, E. S. and Stedinger, J. R.: Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data, Water Resour. Res., 36, 737–744, https://doi.org/10.1029/1999WR900330, 2000. a, b

Marty, C. and Blanchet, J.: Long-term changes in annual maximum snow depth and snowfall in Switzerland based on extreme value statistics, Climatic Change, 111, 705–721, https://doi.org/10.1007/s10584-011-0159-9, 2012. a, b, c

Ménégoz, M., Valla, E., Jourdain, N. C., Blanchet, J., Beaumet, J., Wilhelm, B., Gallée, H., Fettweis, X., Morin, S., and Anquetin, S.: Contrasting seasonal changes in total and intense precipitation in the European Alps from 1903 to 2010, Hydrol. Earth Syst. Sci., 24, 5355–5377, https://doi.org/10.5194/hess-24-5355-2020, 2020. a, b, c

Montanari, A. and Koutsoyiannis, D.: Modeling and mitigating natural hazards: Stationarity is immortal!, Water Resour. Res., 50, 9748–9756, https://doi.org/10.1002/2014WR016092, 2014. a

Nicolet, G., Eckert, N., Morin, S., and Blanchet, J.: Decreasing spatial dependence in extreme snowfall in the French Alps since 1958 under climate change, J. Geophys. Res., 121, 8297–8310, https://doi.org/10.1002/2016JD025427, 2016. a, b

O'Gorman, P. A.: Contrasting responses of mean and extreme snowfall to climate change, Nature, 512, 416–418, https://doi.org/10.1038/nature13625, 2014. a, b

O'Gorman, P. A. and Muller, C. J.: How closely do changes in surface and column water vapor follow Clausius–Clapeyron scaling in climate change simulations?, Environ. Res. Lett., 5, 0–12, https://doi.org/10.1088/1748-9326/5/2/025207, 2010. a

Padoan, S. A., Ribatet, M., and Sisson, S. A.: Likelihood-based inference for max-stable processes, J. Am. Stat. Assoc., 105, 263–277, https://doi.org/10.1198/jasa.2009.tm08577, 2010. a

Quéno, L., Vionnet, V., Dombrowski-Etchevers, I., Lafaysse, M., Dumont, M., and Karbou, F.: Snowpack modelling in the Pyrenees driven by kilometric-resolution meteorological forecasts, The Cryosphere, 10, 1571–1589, https://doi.org/10.5194/tc-10-1571-2016, 2016. a

Rajczak, J. and Schär, C.: Projections of Future Precipitation Extremes Over Europe: A Multimodel Assessment of Climate Simulations, J. Geophys. Res.-Atmos., 122, 773–710, https://doi.org/10.1002/2017JD027176, 2017. a

Revuelto, J., Lecourt, G., Lafaysse, M., Zin, I., Charrois, L., Vionnet, V., Dumont, M., Rabatel, A., Six, D., Condom, T., Morin, S., Viani, A., and Sirguey, P.: Multi-criteria evaluation of snowpack simulations in complex alpine terrain using satellite and in situ observations, Remote Sens.-Basel, 10, 1–32, https://doi.org/10.3390/rs10081171, 2018. a

Scherrer, S. C., Wüthrich, C., Croci-Maspoli, M., Weingartner, R., and Appenzeller, C.: Snow variability in the Swiss Alps 1864–2009, Int. J. Climatol., 33, 3162–3173, https://doi.org/10.1002/joc.3653, 2013. a, b

Sillmann, J., Mischa, C. M., Kallache, M., and Katz, R. W.: Extreme cold winter temperatures in Europe under the influence of North Atlantic atmospheric blocking, J. Climate, 24, 5899–5913, https://doi.org/10.1175/2011JCLI4075.1, 2011. a

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

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

Sun, Q., Zhang, X., Zwiers, F., Westra, S., and Alexander, L. V.: A global, continental, and regional analysis of changes in extreme precipitation, J. Climate, 34, 243–258, https://doi.org/10.1175/JCLI-D-19-0892.1, 2021. a

The New York Times: Madrid Mayor Says Snowstorm Caused Nearly $2 Billion in Damage, The New York Times, https://www.nytimes.com/2021/01/14/world/europe/madrid-snowstorm-damage.html (last access: 7 September 2021), 2021. a

Tramblay, Y. and Somot, S.: Future evolution of extreme precipitation in the Mediterranean, Climatic Change, 151, 289–302, https://doi.org/10.1007/s10584-018-2300-5, 2018. a

Vernay, M., Lafaysse, M., Hagenmuller, P., Nheili, R., Verfaillie, D., and Morin, S.: The S2M meteorological and snow cover reanalysis in the French mountainous areas (1958–present), AERIS [data set], https://doi.org/10.25326/37, 2019. a, b, c, d

Vigneau, J.-P.: 1986 dans les Pyrénées orientales: deux perturbations méditerranéennes aux effets remarquables, Revue géographique des Pyrénées et du Sud-Ouest, 58, 23–54, https://doi.org/10.3406/rgpso.1987.4969, 1987. a

Vionnet, V., Dombrowski-Etchevers, I., Lafaysse, M., Quéno, L., Seity, Y., and Bazile, E.: Numerical weather forecasts at kilometer scale in the French Alps: Evaluation and application for snowpack modeling, J. Hydrometeorol., 17, 2591–2614, https://doi.org/10.1175/JHM-D-15-0241.1, 2016. a

Vionnet, V., Six, D., Auger, L., Dumont, M., Lafaysse, M., Quéno, L., Réveillet, M., Dombrowski-Etchevers, I., Thibert, E., and Vincent, C.: Sub-kilometer Precipitation Datasets for Snowpack and Glacier Modeling in Alpine Terrain, Front. Earth Sci., 7, 1–21, https://doi.org/10.3389/feart.2019.00182, 2019. a, b

Wilcox, C., Vischel, T., Panthou, G., Bodian, A., Blanchet, J., Descroix, L., Quantin, G., Cassé, C., Tanimoun, B., and Kone, S.: Trends in hydrological extremes in the Senegal and Niger Rivers, J. Hydrol., 566, 531–545, https://doi.org/10.1016/J.JHYDROL.2018.07.063, 2018. a

- Abstract
- Introduction
- Snowfall data
- Method
- Result
- Discussion
- Conclusions and outlook
- Appendix A: Quantile–quantile plots
- Appendix B: Semi-parametric bootstrap method
- Appendix C: Sensitivity to the return period
- Appendix D: Trends in 100-year return levels of winter precipitation
- Appendix E: Seasons when the annual maxima of daily snowfall occurred
- Data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Snowfall data
- Method
- Result
- Discussion
- Conclusions and outlook
- Appendix A: Quantile–quantile plots
- Appendix B: Semi-parametric bootstrap method
- Appendix C: Sensitivity to the return period
- Appendix D: Trends in 100-year return levels of winter precipitation
- Appendix E: Seasons when the annual maxima of daily snowfall occurred
- Data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References