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

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}{\mathrm{kg}\mathrm{m}}^{-\mathrm{2}}$) on average at 3500 m and a decrease of −10 % ($-\mathrm{7}{\mathrm{kg}\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.

1 Introduction

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}\mathrm{\%}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$, i.e., $\mathrm{7}\mathrm{\%}{}^{\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}\mathrm{\%}{}^{\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)

Table 1 Temporal trends in extreme snowfall with respect to elevation in and around the French Alps. Elevations are in meters above sea level (m a.s.l.). S_Nd denotes the annual maximum of snowfall in N consecutive days.

Figure 1 (a–d) Time series of annual maxima of daily snowfall from 1959 to 2019 for the Vanoise massif (purple region in e) clustered by the four ranges of elevations considered. (e) Topography and delineation of the 23 massifs of the French Alps (Durand et al., 2009).

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.

2 Snowfall data

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², 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⁻², 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 Method

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,

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }P(Y\le y)=\left\{\begin{array}{l}\exp \left[-{\left(\mathrm{1}+\mathit{\xi }\frac{y-\mathit{\mu }}{\mathit{\sigma }}\right)}_{+}^{-\frac{\mathrm{1}}{\mathit{\xi }}}\right]\\ \text{ if }\mathit{\xi }\ne \mathrm{0}\text{ and where }{u}_{+}\text{ denotes }max(u,\mathrm{0}),\\ \exp \left[-\exp \left(-\frac{y-\mathit{\mu }}{\mathit{\sigma }}\right)\right]\\ \text{ if }\mathit{\xi }=\mathrm{0}\text{; in other words }Y\sim \text{Gumbel}(\mathit{\mu },\mathit{\sigma }),\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{Y}_{z,t}\sim \text{GEV}\left(\mathit{\mu }\right(z,t),\mathit{\sigma }(z,t),\mathit{\xi }(z\left)\right).\end{array}$$

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.

Table 2 Elevational–temporal models considered rely on the GEV distribution. For the elevational non-stationarity, the location and the scale parameters vary linearly with the elevation z, while the shape is either constant or linear with z. For temporal non-stationary models, the location and/or the scale vary linearly with time t. See Sect. 3.2 for a full description of the terms used in the table.

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₁ 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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }P(Y_{z,t}\le y_p(z,t\left)\right|{\widehat{\mathit{\theta }}}_{\mathcal{M}})=\mathrm{1}-p\leftrightarrow y_p(z,t)=\mathit{\mu }(z,t)\\ \text{(3)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathit{\sigma }(z,t)}{\mathit{\xi }\left(z\right)}}\left[\mathrm{1}-{\left(-\log (\mathrm{1}-p)\right)}^{-\mathit{\xi }\left(z\right)}\right].\end{array}$$

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

Figure 2 Changes in GEV parameters (a–c) and in 100-year return levels (d) with the elevation for the 23 massifs of the French Alps. GEV distributions are estimated pointwise for the annual maxima of daily snowfall every 300 m of elevation.

Figure 3 Elevation gradients for the GEV parameters (a–c) and 100-year return levels (d) for the 23 massifs of the French Alps. GEV distributions are estimated pointwise for the annual maxima of daily snowfall every 300 m of elevation. Elevation gradients are estimated with a linear regression. The R² coefficient is written in black for each massif.

Figure 4 Selected models and shape parameter values for each range of elevations in the 23 massifs of the French Alps. We write the suffix of the name of each selected model on the map; e.g., we write μ_t,σ_t for the model ${\mathcal{M}}_{{\mathit{\mu }}_t,{\mathit{\sigma }}_t}$. We underline the suffix when the model has a significant trend (Sect. 3.3). Hatched grey areas denote missing data, e.g., when the elevation is above the top elevation of the massif. Shape parameter values are computed at the middle elevation for each range, e.g., at 1500 m for the range 1000–2000 m.

Figure 5 Comparison of pointwise distributions with our approach based on piecewise elevational–temporal models for the Vanoise massif. GEV parameters (a–c) and 100-year return levels with their 80 % confidence intervals (d) are shown from 600 to 3600 m of elevation.

Figure 6 Percentages of massifs with significant/non-significant trends in 100-year return levels of daily snowfall for each range of elevation. A massif has an increasing/decreasing trend if the 100-year return level of the selected elevational–temporal model has increased/decreased.

Figure 7 (a) Distributions of changes in 100-year return levels between 1959 and 2019 for one elevation in each range of elevation. The mean and the median are displayed with a green triangle and an orange line, respectively. (b) Same as (a) but for the relative changes. Distributions of changes are computed at the middle elevation for each range, e.g., at 1500 m for the range 1000–2000 m.

4 Result

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² 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}{\mathrm{kg}\mathrm{m}}^{-\mathrm{2}}/\mathrm{100}\mathrm{m}$ and $\mathrm{0.39}{\mathrm{kg}\mathrm{m}}^{-\mathrm{2}}/\mathrm{100}\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² 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² 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 ℳ₀ 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}{\mathrm{kg}\mathrm{m}}^{-\mathrm{2}}$) on average at 3500 m and a decrease of −10 % ($-\mathrm{7}{\mathrm{kg}\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}{\mathrm{kg}\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⁻²), 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.

Figure 8 Changes in 100-year return levels of daily snowfall between 1959 and 2019 for each range of elevations. The corresponding relative changes are displayed on the map. Hatched grey areas denote missing data, e.g., when the elevation is above the top elevation of the massif. Changes in return levels are computed at the middle elevation for each range, e.g., at 1500 m for the range 1000–2000 m. Massifs with non-significant trends are indicated with a pattern of whites dots.

Figure 9 The 100-year return levels in 2019 of daily snowfall for each range of elevations. The 100-year return levels are both illustrated with colors (elevation-range-dependent scale) and written on the map. Hatched grey areas denote missing data, e.g., when the elevation is above the top elevation of the massif. Return levels are computed at the middle elevation for each range, e.g., at 1500 m for the range 1000–2000 m.

5 Discussion

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}{\mathrm{kg}\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.

6 Conclusions and outlook

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}{\mathrm{kg}\mathrm{m}}^{-\mathrm{2}}$) on average at 3500 m and a decrease of −10 % ($-\mathrm{7}{\mathrm{kg}\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).

Appendix A: Quantile–quantile plots

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_{\text{GEV}\to \text{Standard Gumbel}}(Y_{z,t};{\mathit{\theta }}_{\mathcal{M}})=\\ \text{(A1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\left(z\right)}}\log \left(\mathrm{1}+\mathit{\xi }\left(z\right){\displaystyle \frac{Y_{z,t}-\mathit{\mu }(z,t)}{\mathit{\sigma }(z,t)}}\right)\sim \text{Gumbel}(\mathrm{0},\mathrm{1}).\end{array}$$

Figure A1 Q–Q plots of the selected elevational–temporal models for the Vanoise massif for the four ranges of elevations considered (see Fig. 1 for the time series). (a) Below 1000 m. (b) Between 1000 and 2000 m. (c) Between 2000 and 3000 m. (d) Above 3000 m.

The transformed observations, a.k.a. residuals, are denoted as ${\mathit{ϵ}}_{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{ϵ}}_{z_{\mathrm{1}},t_{\mathrm{1}}},\mathrm{\dots },{\mathit{ϵ}}_{z_{\mathrm{1}},t_M},\mathrm{\dots },{\mathit{ϵ}}_{z_N,t_{\mathrm{1}}}\mathrm{\dots },{\mathit{ϵ}}_{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).

Appendix B: Semi-parametric bootstrap method

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{ϵ}}_{z_{\mathrm{1}},t_{\mathrm{1}}},\mathrm{\dots },{\mathit{ϵ}}_{z_{\mathrm{1}},t_M},\mathrm{\dots },{\mathit{ϵ}}_{z_N,t_{\mathrm{1}}}\mathrm{\dots },{\mathit{ϵ}}_{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{̃}}{\mathit{ϵ}}}_{z_{\mathrm{1}},t_{\mathrm{1}}},\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{ϵ}}}_{z_{\mathrm{1}},t_M},\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{ϵ}}}_{z_N,t_{\mathrm{1}}}\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{ϵ}}}_{z_N,t_M}$. Afterwards, we transform these bootstrapped residuals into bootstrapped annual maxima as follows: ${\stackrel{\mathrm{̃}}{y}}_{z,t}=f_{\text{GEV}\to \text{Standard Gumbel}}^{-\mathrm{1}}({\stackrel{\mathrm{̃}}{\mathit{ϵ}}}_{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{̃}}{y}}_{z_{\mathrm{1}},t_{\mathrm{1}}},\mathrm{\dots },{\stackrel{\mathrm{̃}}{y}}_{z_{\mathrm{1}},t_M},\mathrm{\dots },{\stackrel{\mathrm{̃}}{y}}_{z_N,t_{\mathrm{1}}}\mathrm{\dots },{\stackrel{\mathrm{̃}}{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)}}$.

Appendix C: Sensitivity to the return period

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.

Figure C1 Percentages of massifs with significant/non-significant trends in return levels of daily snowfall for each range of elevation and for a return period equal to (a) 10, (b) 50, or (c) 100 years. A massif has an increasing/decreasing trend if the return level of the selected elevational–temporal model has increased/decreased.

Appendix D: Trends in 100-year return levels of winter precipitation

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.

Figure D1 Changes in GEV parameters (a–c) and in 100-year return levels (d) with the elevation for the 23 massifs of the French Alps. GEV distributions are estimated pointwise for the annual maxima of daily winter precipitation every 300 m of elevation.

Figure D2 Changes in 100-year return levels of daily winter precipitation between 1959 and 2019 for each range of elevations. The corresponding relative changes are displayed on the map. Hatched grey areas denote missing data, e.g., when the elevation is above the top elevation of the massif. Changes in return levels are computed at the middle elevation for each range, e.g., at 1500 m for the range 1000–2000 m. Massifs with non-significant trends are indicated with a pattern of white dots.

Appendix E: Seasons when the annual maxima of daily snowfall occurred

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.

Figure E1 Seasons when the annual maxima of daily snowfall occurred for elevation range 1 (below 1000 m), elevation range 2 (between 1000 and 2000 m), elevation range 3 (between 2000 and 3000 m), and elevation range 4 (above 3000 m).