Accounting for Non-stationarity in Extreme Snow Loads: a Comparison with Building Standards in the French Alps

In a context of climate change, trends in extreme snow loads need to be determined to minimize the risk of structure collapse. We study trends in annual maxima of ground snow load (GSL) using non-stationary extreme value models. Trends in return levels of GSL are assessed at a mountain massif scale from GSL data, provided for the French Alps from 1959 to 2019 by a meteorological reanalysis and a snowpack model. Our results indicate a temporal decrease in 50-year return levels from 900 m to 4200 m, significant in the Northwest of the French Alps until 2100 m. Despite this decrease, in half of the 5 massifs, the return level in 2019 at 1800 m exceeds the return level designed for French building standards under a stationary assumption. We believe that this high number of exceedances is due to questionable assumptions concerning the computation of current standards. For example, these were devised with GSL, estimated from snow depth and constant snow density set to 150 kg m−3, which underestimate typical GSL values for the full snowpack.

Snowpack variables related to GSL (snow depth, snow water equivalent) evolve with climate change. As shown in Table 1, literature on past trends in snowpack variables for the Western Alps shows a decreasing trend. Literature on projected trends also 25 points to a decrease (stronger for the second half of the 21st century under a high greenhouse gas emission scenario than with strong reductions in greenhouse gas emissions) for mean winter (December-May) snow water equivalent in the European Alps (IPCC, 2019). However, anthropogenic climate change impacts climatic variables in their averages, but also in their extremes (G Klein Tank and K Nnen, 2003;IPCC, 2012). For instance, annual maxima of snow depth have decreased in Switzerland (Marty and Blanchet, 2012). Projected trends in extreme snowpack variables are prone to strong uncertainties (Strasser, 2008; 30 Beniston et al., 2018) as both mean winter temperature (IPCC, 2019) and winter precipitation extremes (Rajczak and Schär, 2017) are projected to increase in the European Alps.  1931-1999Laternser and Schneebeli (2003 Winter mean (Dec to Feb) Decrease in the North  Table 1. Past trends in snow depth (HS) and snow water equivalent (SWE) according to existing studies in the Western Alps, i.e. in Italy (IT), France (FR), and Switzerland (CH). In the trend column, "North" and "South" refer to the considered country.
The impact of climate change on GSL was not taken into account in current European standards for structural design, a.k.a Eurocodes (Sanpaolesi et al., 1998), which drive French standards (Biétry, 2005). These standards define that structures must withstand their own weight plus a pressure proportional to a characteristic value. The latter is the stationary 50-year return 35 level of GSL, i.e. exceeded once every 50 years in average. Thus, studying trends in 50-year return levels of GSL is needed for updating these standards (Croce et al., 2018). In the literature, past and projected trends in 50-year return levels of GSL have rarely been investigated with the exception of (Rózsás et al., 2016;Il Jeong and Sushama, 2018;Croce et al., 2018). In the French Alps, several studies focused on extreme snow variables (Biétry, 2005;Gaume et al., 2012Gaume et al., , 2013 and their spatial dependence (Nicolet et al., 2015(Nicolet et al., , 2016(Nicolet et al., , 2017(Nicolet et al., , 2018. However, trends in 50-year return levels of GSL remain unexplored. 40 We fill these gaps by studying annual maxima of GSL provided every 300 m of elevation at a mountain massif scale for the 23 French Alps massifs. We rely on a reanalysis produced by the SAFRAN-Crocus chain (Durand et al., 2009a;Vionnet et al., 2012) available for the period 1959-2019. The major advantage of this reanalysis is to take benefit of an advanced snowpack model which provides daily estimates of ground snow load values, while previous studies relied on approximate values directly related to snow depth with a crude estimation of snow density (Biétry, 2005). Our methodology consists in applying stationary 45 and non-stationary extreme value models to annual maxima time series. We select one model by massif and altitude with the AIC statistical criterion, and assess its significance with the likelihood ratio statistical test. Finally, for each massif and altitude, we compute the relative change in 50-year return level between 1960 and 2010, and we compare the non-stationary return level in 2019 with the stationary return level designed for French building standards. Our approach considers only natural snow processes, i.e. we do not account for snow removal throughout the year and consider all processes (accumulation, thaw/freeze, 50 melt, compaction etc.) occurring during the winter season. This paper is organized as follows. Section 2 presents our data. Section 3 describes standards for ground snow load. Then, section 4 explains our methodology. Results, discussion and conclusions are introduced in Sections 5, 6 and 7, respectively.

Ground snow load data
The study area covers the French Alps which are located between Lake Geneva to the north and the Mediterranean Sea to the 55 south ( Fig. 1). The climate is contrasted, colder and wetter in the northern Alps and much drier in the southern Alps (Durand et al., 2009a). This region is typically divided into 23 mountain massifs. In this work, we rely on a reanalysis (Vernay et al., 2019) from the SAFRAN-Crocus chain (Durand et al., 2009a;Vionnet et al., 2012) available from August 1958 to July 2019 at the scale of these massifs, for 300 m elevation bands from 300 m to 4800 m. However, annual maxima are only available from 1959 to 2019. Indeed, in this work, annual maxima denotes the maxima during a year centered on the winter season, e.g.  SAFRAN meteorological reanalysis (Durand et al., 2009a) provides consistent meteorological data (precipitation, temperature, humidity, radiation, wind speed) over the considered mountain massifs and elevations. The Crocus snowpack model (Vionnet et al., 2012) infers snow depth and snow water equivalent (SWE) based on SAFRAN time series. Crocus is a one-processes such as thermal diffusion, phase changes and metamorphism. A large intercomparison of snow models illustrates that Crocus is among the top models to simulate SWE (Krinner et al., 2018). Furthermore, another recent intercomparison of SWE products emphasizes Crocus usefulness as it concludes that ensembles containing Crocus (and/or Crocus driven by ERA-Interim reanalysis) perform better than those that do not (Mortimer et al., 2019).
In this paper, GSL equals the SWE from the SAFRAN-Crocus chain times the gravitational acceleration. We study time 70 series of annual maxima of GSL for each massif from 1959 to 2019 for 300 m elevation bands from 300 m to 4800 m (Fig. 1).

Standards for ground snow load in the French Alps
GSL French standards (Biétry, 2005) are based mostly on Eurocodes (Sanpaolesi et al., 1998) and on prior French standards (AFNOR, 1996). Each French department, and by extension each French Alps massif, is associated with a region (C or E) that  French standards were elaborated with annual maxima time series of snow depth on the ground measured in stations from 1945 to 1992. GSL data were approximated from annual maxima of snow depth and by assuming that snow density equals 150 kg m −3 . Following Eurocodes, characteristic value of GSL is defined as the 50-year return level of a Gumbel distribution (Sect. 4.1). This distribution was fitted using the least squares method and by removing the top annual maximum when considered exceptional (Biétry, 2005).

80
Following extreme value theory, we employ 2 stationary models and 6 non-stationary models for time series of annual maxima of GSL (Sect. 4.1). We select a single model for each time series (i.e. for each massif and altitude) with the AIC statistical criterion, and assess its significance with the likelihood ratio statistical test (Sect. 4.2). Finally, we compute the relative change in 50-year return levels of GSL, and quantify their uncertainty (Sect. 4.3). Climate extremes are generally studied with statistics. As underlined in the IPCC special report on climate extremes, a large amount of statistical literature builds on extreme indices to examine moderate extremes (IPCC, 2012). However, since we focus on rarer extremes, it is recommended to rely on extreme value theory (EVT) (Coles, 2001). Such statistical models provide and hypothesize additional prior information in order to compensate the limited amount of empirical observations that commonly 90 span only several decades. These models can be used to extrapolate beyond the empirical observations, and to estimate return levels (Sect. 4.3).
EVT offers a suitable framework to analyse extreme values, i.e. to model the form of the tail for almost any probability distribution. Asymptotically, as the central limit theorem motivates sample means modelling with the normal distribution, the Fisher-Tippett-Gnedenko theorem (Fisher and Tippett, 1928;Gnedenko, 1943) encourages sample maxima modelling with 95 the GEV distribution. This theorem justifies that the maximum of finite-sized blocks with a large enough block size can be modeled with the GEV distribution. In practice, an annual maximum is thus usually considered as a realization of a GEV distribution. Three parameters define the GEV distribution: a location µ, a scale σ > 0 and a shape ζ (a.k.a extremal index or tail index). The GEV distribution includes three specific types of distributions: Weibull (ζ < 0), Fréchet (ζ > 0) and Gumbel (ζ = 0). Thus, by definition, if Y represents annual maximum of GSL, we can assume that Y follows a GEV distribution, i.e.

100
Y ∼ GEV(µ, σ, ζ), which implies that: In a context of climate change, a large amount of hydrological literature builds on non-stationary modelling (Milly et al., 2008) to assess whether a time series is generated by a unique probability distribution (stationary model), or if the generating probability distribution is changing (non-stationary model). Non-stationary extremes are usually studied with both non-105 stationary modelling and EVT (Katz et al., 2002). Annual maxima are assumed independent but not necessarily identically distributed (Serinaldi and Kilsby, 2015). Such approaches combine a stationary random component (a fixed extreme value distribution) with non-stationary deterministic functions that map each temporal covariate t to the changing parameters of the distribution (Montanari and Koutsoyiannis, 2014). In a non-stationary context, Zhang et al. (2004) showed that tests based on this parametric approach have stronger power of detection when compared with non-parametric methods.

110
In this paper, we consider non-stationarity for both the Gumbel distribution and the more general GEV distribution, since they represent natural extensions of the Gumbel distribution which was used for French building standards (Sect. 3). For any , as the Gumbel distribution correspond to ζ(t) = 0. For a model M, we denote as θ M all parameters for its functions (µ(t), σ(t) and ζ(t)). We focus on simple linear functions due to the limited length of time series (60 years). The linearity starts in 1959 which is the first winter with available data. As shown in Table 2, 115 we consider only models with a constant shape parameter, but where the location and/or the scale parameter can vary linearly with years t. Table 2. Statistical models considered for annual maxima of GSL are based on the Gumbel or the GEV distribution, and are extensions of the stationary Gumbel model. For non-stationary models, the location and/or the scale vary linearly with years t after the starting year 1959.

Model selection
Let y = (y 1959 , ..., y 2019 ) represents a time series of annual maxima of GSL, i.e. for a massif and an elevation band (Sect. 2). First, models are fitted with the maximum likelihood method. For every model M, we compute the maximum likelihood 120 estimator θ M which corresponds to the parameter θ M that maximizes the likelihood: Then, for each y, i.e. for each massif and elevation band, we select the model M N with the minimal AIC value (Akaike, 1974), as it is the best information criterion in a non-stationary context with small sample sizes (Kim et al., 2017). We define that: The selected model M N can be any model from where∼ means distributed under suitable regularity conditions. In practice, the test works as follows. We first choose a 0.05 level of significance. Then, if LR is greater than q χ 2 N , the 135 1 − 0.05 = 0.95 quantile of the χ 2 N distribution, it means that we reject the nested model M 0 in favor of the selected model M N . In this case, if the selected model M N is non-stationary, then we consider the associated trend as significant.

Return level
In a stationary context, the T -year return level, which corresponds to a return period of T years, is the classical metric to quantify hazards of extreme events (Cooley, 2012). For a stationary model, there is a one-to-one relationship between a return 140 level (a quantile exceeded each year with probability p) and a return period (a duration exceeded every T = 1 p years in average). In a non-stationary context, return level and return period concepts (Cooley, 2012) become further ambiguous, prone to misconceptions and can lead to misleading conclusions (Serinaldi, 2015). We focus on the yearly level for a fixed probability of exceedance, a.k.a effective return level (Katz et al., 2002;Cheng et al., 2014), as it conveys best that hazard evolves with time.

145
For the stationary Gumbel model M 0 , the return level z p (θ M0 ) is defined as the level exceeded each year with probability p. In other words, if Y denotes an annual maximum, then P (Y ≤ z p (θ M0 )) = 1 − p. This return level is constant through time and equals z p (θ M0 ) = µ 0 − σ 0 log (− log (1 − p)). In this paper, we set p = 1 50 = 0.02 as it corresponds to the 50-year return period defined by French standards (based on European standards) for the design working life of buildings (Sect. 3).
For the selected model M N , we define our return level as the yearly level for a fixed probability of exceedance p. For any 150 model considered in Table 2, we obtain that z where we set µ 1 , σ 1 or ζ 0 to 0 if they are not defined in the model M N . For example, for the Gumbel model M 0 , the return level is constant: for any year t, (1 − p)). For any considered model, the time derivative of the return level is constant, as . It quantifies the yearly change in return level. Thus, the relative difference of return levels between year t 1 and year t 2 is: In the context of maximum likelihood estimation, uncertainty related to return levels can be derived by the delta method, which quickly provides confidence intervals both in the stationary and non-stationary case (Coles, 2001;Gilleland and Katz, 2016). First, the return level estimator associated to the maximum likelihood estimator simply equals z p ( θ M ). Then, due to the asymptotic normality of the maximum likelihood estimator (MLE), we can assume that, even with a finite number of data, the MLE is normally distributed. Therefore, we have that under regularity conditions, limits of the 1 − α = 95% confidence where q α 2 is the 1α 2 quantile of the standard normal distribution, and v zp is a function that associates to each parameter θ M the variance of the approximate normal distribution associated to its return level z p (θ M ). For a full expression of the function v zp and for details on the delta method, we refer to Theorem 2.4 of Coles (2001). In particular, this theorem explains that the delta method is valid for ζ 0 < 1, which is respected in our case as −0.5 ≤ ζ 0 ≤ 0.5 (Sect. 5).

165
Also, uncertainty of non-stationary return levels z p ( θ M , t) can be obtained by incorporating the covariate t in the function z p .

Application
First, we exclude 4 times series of annual maxima with more than 10% of zeros, i.e. years without GSL. Then, we fit models on time series, and retain only those with −0.5 ≤ ζ 0 ≤ 0.5. This impacts 3 time series. We made this choice because ζ 0 > 0.5 designates distributions with an "exploding" tail which are known to be physically implausible (Martins and Stedinger, 2000). 170 We select one model for each time series (i.e. for each massif and altitude) with the AIC statistical criterion, and assess its significance with the likelihood ratio test (Sect. 4.2). Figure 3 shows that a stationary model, i.e. models M 0 and M ζ0 , is selected for a majority (57%) of time series studied (Sect.

175
2). Models with a linearity in both the location and scale parameters are the most frequently selected non-stationary models (22%). For both stationary and non-stationary models, Gumbel models are always more selected that their corresponding GEV models (Fig. 3, Fig. 4). All in all, we highlight that 40% of selected models are significantly different from the stationary Gumbel model M 0 .  massifs are white-colored illustrating that a (stationary or non-stationary) Gumbel model (i.e. ζ 0 = 0) is selected (Sect. 5). It emphasizes that a Gumbel distribution often explains more succinctly the data than a GEV distribution. Also, with the GEV distribution, the estimated most likely shape parameter ζ 0 is often quite uncertain, i.e. confidence intervals are large, which is the main reason why French standards did not rely on it. This uncertainty comes from the limited length of time series, and longer time series would certainly reduce it, and thus estimate the most likely shape parameter more robustly.

185
On Figure 4, we further observe that non-null shape parameters at low altitudes (≤ 900 m) are always positive (browncolored massifs), i.e. a Fréchet distribution is preferred. On the other hand, for high altitudes (≥ 1500 m) non-null shape parameters are always negative (green-colored massifs), i.e. a Weibull distribution is favored. We hypothesize that this might be due to the different nature of annual maxima of GSL between low and high altitudes. At high altitudes, annual maxima are mainly due to snowpack accumulation during several months, while at low altitudes this accumulation is limited, and thus 190 annual maxima roughly correspond to heavy precipitations.    Figure 6 emphasizes the evolution of decreasing trends between 900 m and 4800 m of altitude. We observe that decreasing trends are significant for more than one-third of the massifs, located in the Northwest of the Alps (Sect. A), until 2100 m (black 200 bars). In half a century, return levels have dropped in average by up to 30% at 900 m. At higher altitudes we observe a decline in the percentage of massifs with a significant decreasing trend, which is non null up to 3300 m. However, for both the relative decrease and the percentage of massifs with a decreasing trend, even if we notice a similar declining pattern, we also detect a slight growth above 3000 m. This echoes results from Lüthi et al. (2019), which found that, in the Alps above 3000 m, the relative decrease for projected winter-mean of fresh snow water equivalent is more marked than at 3000 m (see their Figure 8). 205 We emphasize, however, that most meteorological observations used as input to the SAFRAN-Crocus reanalysis are situated below 2000 m. Therefore, trends beyond 2000 m altitude should be considered with great caution.  Indeed, we found that these annual maxima correspond to snow load accumulated in few days, and thus are mainly driven by heavy precipitations rather than a full season snowpack accumulation. In particular, we hypothesize that the 2 massifs with important increasing trends (red color) might be caused by a local phenomenon called "East return" which is a low pressure system coming from the Mediterranean sea, that might be intensifying with global warming (Garavaglia et al., 2010;Faranda, 2019).

215
To sum up trends in return levels of ground snow load: from 900 m to 4800 m, either no trend or decreasing trend of 50-year return levels of GSL are found (Fig. 5, Fig. 6 and Fig. A1). while at 300 m and 600 m, no clear trends are found (Fig. 7).

Comparison of return levels of ground snow load with French standards
Every 300 m of altitude, from 300 m to 1800 m, we compute 50-year return levels and their uncertainty from GSL data (Sect.   Figure 9 sums up the comparison between French standards and our results. We display the percentage of massifs whose 225 return level estimated from data exceeds return level from standards. The number of massifs with available data is equal to 11 at 300 m, 18 at 600 m, and 23 at 900 m and above. Limits of the confidence intervals (black bars) are computed as the percentage of exceedances for the limits of return levels' 95% confidence interval (black bars on Fig. 8)..
First, if we estimate return levels from data with French standards method (Fig. 9 Left), i.e. with a stationary Gumbel model M 0 , and GSL data approximated with snow depth obtained from reanalysis and ρ SNOW = 150 kg m −3 , then we observe few 230 exceedances (always less than 10%). Thus, in this setting, estimations from our reanalysis are consistent with French standards.
However, if we consider the actual GSL, i.e. computed with the snow water equivalent from reanalysis, then French standards drastically underestimate return levels. Indeed, with a stationary Gumbel model M 0 , then for altitudes above or equal to 900 m, French standards are exceeded for a majority of massifs ( Fig. 9 Center). But, if we consider the selected model M N , i.e.
if we account for the decreasing trend in 50-year return levels, we have less exceedances at all altitudes ( Fig. 9 Right). In the 235 latter case, at worst, i.e. at 1800 m, half of the massifs still exceed French standards return levels.

Methodological considerations
We discuss in depth the statistical models chosen for this study. It is well-known that an annual maximum based approach can be wasteful in terms of data (Coles, 2001). However, since our objective is to estimate 50-year return levels and since 240 we have 60 years of data, we still decide to rely on the annual maximum based approach (with the GEV distribution) rather than on the concurrent approach based on threshold exceedances (with the Generalized Pareto distribution). Also, with the GEV distribution (and its particular case: the Gumbel distribution) our methodology is a direct extension of French building standards (Sect. 3).
For the non-stationary models, we focus on simple deterministic function of time (µ(t), σ(t), ζ(t)) due to the limited length 245 of time series. A linear non-stationarity seems preferable to a non-stationarity based on the Heaviside step function due to the continuous nature of climate change. We start the linear non-stationary at the initial year, i.e. 1959. We tried to start the nonstationarity only after a most likely year  but our results lacked coherence both between close massifs (no clear spatial pattern for the most likely years) and w.r.t. the altitude (no clear trend in most likely years as altitude rises).
We decided to consider non-stationarity only for the location and scale parameter. Indeed, in the literature, a linear non-250 stationarity is considered sometimes only for the location parameter (Fowler et al., 2010;Tramblay and Somot, 2018) but more often both for 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). Another reason for which we considered a non-stationarity 13 https://doi.org/10.5194/nhess-2020-81 Preprint. Discussion started: 27 March 2020 c Author(s) 2020. CC BY 4.0 License.
for both parameters is that we did not find the scale parameters to be proportional to the location parameters, which could have otherwise simplified our parametrization. Also, the shape parameter is typically considered constant, with few exceptions in 255 the literature and we follow this approach.
For time series containing zeros, French standards rely on a mixed discrete-continuous distribution. We did not rely on this choice because for time series containing zeros considered in this study, i.e. with less than 10% of zeros (Sect. 5), we almost obtain the same 50-year return levels with this distribution than with our approach (absolute difference remains lower than 0.1 kN m −2 ). 260 6.2 On the limitation to approximate annual maxima of ground snow load with annual maxima of snow depth Snow water equivalent (SWE) times the gravitational constant equals GSL. However, most countries do not measure SWE but only have access to snow depth (HS) (Haberkorn et al., 2019). In that case, snow density is required to obtain SWE (and subsequently GSL) from HS (Sect. 1). In particular, French standards approximate annual maxima of GSL with annual maxima of HS and by assuming a constant snow density, equal to ρ SNOW = 150 kg m −3 . We highlight limitations of such approaches 265 with our reanalysis (Sect. 2) that provides, for the whole snowpack, daily values of SWE, HS, and thus of snow density.
We find that annual maxima of GSL is always underestimated by French standards' approximation ( Fig. 10 Left). Main reason is that, when annual maxima of GSL is reached, snow density is in average largely superior to ρ SNOW = 150 kg m −3 ( Fig. 10 Center). Indeed, we observe that at the time of the annual maxima of GSL the snow density is around ≈ 350 kg m −3 in average at 2700 m, and close to ≈ 250 kg m −3 in average at 900 m. Finally, despite high variations along the years, we also 270 notice that, when annual maxima of GSL is reached, snow depth can be much lower than the annual maxima of snow depth ( Fig. 10 Right), which is another argument against the use of snow depth maxima as a proxy for GSL maxima. Based on both a reanalysis and a snowpack model, we detect an overall temporal decreasing trend w.r.t to 50-year return levels of ground snow load (GSL) between 900 m and 4200 m, which is significant until 2100 m in the Northwest of the French Alps.

275
This confirms other studies in the Western Alps which also found overall decreasing trends in linked snowpack variables: snow water equivalent and snow depth. Despite decreasing return levels, in half of the massifs the 50-year return level in 2019 at 1800 m exceeds the stationary return level designed for French standards.
We hypothesize that this amount of exceedances might be due to an underestimation of GSL by French standards. Indeed, these standards were devised with GSL estimated from snow depth and constant snow density equal to 150 kg m −3 , which 280 underestimate typical GSL values for the full snowpack. Another reason for these exceedances might be ill-designed relationships between altitude and snow load. As shown on Fig. 2, French standards return levels augment linearly by parts w.r.t the altitude. Indeed, French standards (Biétry, 2005) follow previous national standards ( AFNOR, 1996) that advised for a linear relationship between altitude and snow load instead of relying on European standards' results that showed a quadratic relationship for the Alpine Region (Sanpaolesi et al., 1998). Thus, at higher altitudes, French standards might underestimate actual 285 return levels which might explain the augmenting percentage of exceedance that we observe w.r.t the altitude (Fig. 9 Right).
Many potential extensions of this work could be considered. First, our methodology could be extended with more advanced definitions of non-stationary return levels (Rootzén and Katz, 2013;Serinaldi, 2015). Also, instead of considering time series of annual maxima as independent, we believe that our analysis may benefit from an explicit modelling of the spatial dependence between extremes. Then, reanalyses are increasingly available at the European scale (e.g. Soci et al. (2016)), which could be 290 used for extending this work at a wider geographical scale. This requires, however, to remain cognizant of the limitations of such reanalyses, in particular (i) the temporal heterogeneity of the meteorological data input to these reanalyses (Vidal et al., 2010), (ii) the lack of observations at high elevations, requiring caution in analyzing trends for high elevation locations and (iii) model errors (e.g., snowpack model errors) which need to be taken into account when analyzing the results.
Finally, even if, according to our analysis, GSL exceeds French standards return levels in the French Alps, (Fig. 9 Right), 295 few destructions related to snow loads actually occurred. Several reasons might explain that. First, French standards consider a coefficient that maps GSL return levels to roof snow load return level, i.e. multiplication by a coefficient that summarizes several roof features: shape, exposure and thermal transmission (Sanpaolesi et al., 1998)). This coefficient might be overprotective. Also, following European standards, to ensure roofs reliability, designers actually build roofs that withstand return levels of roof snow load multiplied by a safety coefficient (1.5). Above all, French standards do not take into account that, after intense 300 days of snowfall, the snow accumulated on the roof either slides off or is removed. In that case, the main risk lies in extreme snow events that might accumulate in few days enough snow to exceed French standards. Undeniably, most known snow load destructions resulted from such intense snow events, sometimes combined with liquid precipitation that often heavily increase snow load. The response of these short but extreme and complex snow events to climate change might be an interesting topic for future research. https://doi.org/10.5194/nhess-2020-81 Preprint. Discussion started: 27 March 2020 c Author(s) 2020. CC BY 4.0 License.