Towards a reproducible snow load map – an example for Austria

The European Committee for Standardization defines zonings and calculation criteria for different European regions to assign snow loads for structural design. In the Alpine region these defaults are quite coarse; countries therefore use their own products, and inconsistencies at national borders are a common problem.

A new methodology to derive a snow load map for Austria is presented, which is reproducible and could be used across borders. It is based on (i) modeling snow loads with the specially developed Δsnow model at 897 sophistically quality controlled snow depth series in Austria and neighboring countries and (ii) a generalized additive model where covariates and their combinations are represented by penalized regression splines, fitted to series of yearly snow load maxima derived in the first step. This results in spatially modeled snow load extremes.

The new approach outperforms a standard smooth model and is much more accurate than the currently used Austrian snow load map when compared to the RMSE of the 50-year snow load return values through a cross-validation procedure. No zoning is necessary, and the new map's RMSE of station-wise estimated 50-year generalized extreme value (GEV) return levels gradually rises to 2.2 kN m⁻² at an elevation of 2000 m. The bias is 0.18 kN m⁻² and positive across all elevations. When restricting the range of validity of the new map to 2000 m elevation, negative bias values that significantly underestimate 50-year snow loads at a very small number of stations are the only objective problem that has to be solved before the new map can be proposed as a successor of the current Austrian snow load map.

1 Introduction

The current European standard for structural design – Eurocode EN 1991-1-3 (CEN, 2015) – defines rules to calculate acceptable snow loads on constructions and buildings. It is based on snow load maps for different climate regions, defining snow load zones. Within these zones a linear or quadratic “elevation–snow load relationship” is used to calculate s_k, the so-called characteristic snow load on the ground, which in turn is based on an annual probability of exceedance of 0.02 and referred to as the return level of snow load (in kN m⁻²) of a 50-year return period. EN 1991-1-3 defines, for example, an “Alpine region” reproduced in Fig. 1, where $s_k=(\mathrm{0.642}\cdot Z+\mathrm{0.009})\cdot \left[\mathrm{1}+{\left(\frac{A}{\mathrm{728}}\right)}^{\mathrm{2}}\right]$, with Z being the zone number as defined in the legend of Fig. 1 and A being elevation of the location for which s_k should be calculated (A≤1500 m a.s.l.). All strategies outlined in EN 1991-1-3 for how to consider the various influences of roof shapes, thermal properties, wind drift, etc. are ultimately connected to s_k, making it the most crucial snow-load-related quantity for structural design.

Figure 1 Alpine region snow loads at sea level as provided by the Eurocode EN 1991-1-3 (CEN, 2015). Together with the respective elevation–snow load relationship (see text), the characteristic snow load on the ground (s_k) can be calculated for every location below 1500 m elevation.

The aim of this default snow load definition by the European Committee for Standardization (CEN) is “to eliminate or reduce the inconsistencies of snow load values in CEN member states and at borderlines between countries” (CEN, 2015). It results from comprehensive research work of the late 1990s under contracts of the European Commission (Sanpaolesi et al., 1998, 1999). However, Annex C of EN 1991-1-3, where maps and s_k relations are defined, has an informative character, not a normative. Member states can define their own maps, as has already been done before the Eurocode suggestions were developed. In Austria, ÖNORM B 4013, (ON, 1983), for example, already provides a map with much more attention to detail than the Eurocode map. Snowfall in Austria is characterized by a complex climatic setting (Auer et al., 2006) as well as the country's structured topography, leading to small-scale shading and enhancement effects for precipitation (e.g., Wastl and Zängl, 2010), frequent valley inversions, and spatially and temporally highly diverse snow lines (e.g., Unterstrasser and Zängl, 2006). In 2006, when the Austrian standard was based on Eurocode for the first time (ON, 2006), the zoning was adopted to make use of the abovementioned relation. Still, the richness of detail of the map was kept rather than using the coarse Eurocode suggestion (reprinted in Fig. 1; details on the development of the Austrian snow load standard are given in Winkler et al., 2020). Austria's neighbors also use their own snow load maps (e.g., DIN, 2019; SIA, 2020; Autonome Provinz Bozen Südtirol, 2002). Inconsistencies at borderlines between countries are very common. The underlying methods of how snow load maps are developed are roughly comparable: statistical analysis of yearly snow depth maxima using extreme value theory, assigning snow densities, regionalization into climatic snow load regions, etc. (the final draft for the coming Eurocode – CEN, 2020, Annex A.3 – provides guidance), but in the details they are very different. Not least, some of these methods may be subjective and not reproducible (e.g., the drawing of the zones), and others might not be applied in a transparent way (e.g., how and why certain snow density assumptions are made). In this respect CEN's Eurocode EN 1991-1-3 (CEN, 2015) currently misses its target (cf. Croce et al., 2019). Furthermore, the current Austrian standard (ASI, 2018) suffers from an old database that largely ends in the 1980s. Besides the mentioned methodical inadequacies concerning most standards, this particularly offers motivation to update the Austrian snow load map (see also CEN, 2020, Annex A.5).

Here we propose a novel methodology which will lead to an updated snow load map for Austria. It is based on updated snow data and on transparent and reproducible methods described in Sect. 2 that could also be used in other countries. Section 3 provides comparisons to another methodical candidate and the current Austrian snow load standard as well as accuracy measures. Concluding remarks and an outlook on possible future developments are given in Sect. 4.

2 Data and methods

Long-term snow depth data from 2740 stations in Austria and nearby (ca. 50 km) in Germany, Switzerland, Italy and Slovenia were collected. The data were corrected in many ways. Only a few tasks are shortly outlined here for a more comprehensive illustration the reader is referred to the report on the “Schneelast.Reform” project (Winkler et al., 2021b). A common observational issue is to simply add new snow values over several days and use that as snow depth. Such cases were detected by comparing snow depths with empirical elevation-related thresholds. The wrong values were then corrected by assuming an empirical densification rate with the aid of neighboring stations. One of the essential prerequisites of the Δsnow model is a gapless snow depth record. For Austrian stations gaps up to 31 consecutive days were filled using the dual-layer snow model SNOWGRID (Olefs et al., 2013, 2017). To replace missing observations at a station, the corresponding snow depth time series of the closest SNOWGRID grid point (maximum distance 0.71 km) was linearly interpolated to the observations confining the gap. Years with gaps larger than 31 d were rejected. For non-Austrian stations years with data gaps longer than 6 d were not used, up to 6 d the snow depth records were filled by linear interpolation. About 1.8×10⁷ daily snow depth values were used; 8400 (0.05 %) were filled by the procedure described. At two stations a new snow depth maximum for one season was created; at one station the above-described gap filling led to the second-highest value in one season. Finally, snow depth series from 897 stations fulfilled the following requirements:

• at least 30 years of regular daily snow depths between 1960 and 2019;

• latest year with data is 2009 or later;

• no multiple stations within 0.1^∘ latitude or longitude and within 100 m elevation difference.

The stations lie between 118 and 3109 m elevation (median: 655 m; Fig. 2), and the record lengths are between 30 and 59 years, with a median of 49 years. They cover most of the eastern Alps; see Fig. 3.

Figure 2 Relative number of stations in elevation classes, with the actual number on top of each bar (897 stations).

Figure 3 Values of 50-year snow load return from Δsnow-modeled yearly maxima locally fitted to a GEV distribution for all 897 stations. Elevation is given in grayscale from 100 m (dark) to 3800 m (light).

The newly derived set of snow depth data series is the most comprehensive of its kind in Austria. No other snow dataset contains so many records with regular data at daily resolution for so many years – the most important prerequisites for deriving return levels of snow load. Individual single-day snow depth values are less important than completeness and capturing the yearly maxima as accurately as possible. Snow depth data still provide the base for snow load modeling rather than snow load or water equivalent observations, for which only too short, too sparse and too irregular series exist (cf. Winkler et al., 2021a). However, using snow depth as a starting point necessitates the simulation of snow mass by applying a snow (density) model.

2.1 Snow model

Former developments of snow load standards used constant snow densities or parameterized elevation, snow type and/or age dependencies. An overview is given in Sanpaolesi et al. (1998). The methods are very heterogeneous. European standards also suggest snow densities (e.g., CEN, 2015; ASI, 2018). Given the complex layered structure and viscoelastic behavior of snow as a porous material very close to its pressure melting point, these simple density assumptions fall short. Also “empirical regression models” (cf. Winkler et al., 2021a) like Jonas et al. (2009), Sturm et al. (2010) and Guyennon et al. (2019) are not applicable since they fail to model yearly maxima sufficiently well (Winkler et al., 2021a). On the other hand, sophisticated “thermodynamic snow models” like SNOWPACK (Lehning et al., 2002) or Crocus (Vionnet et al., 2012) need meteorological forcing that is not available for those long-term snow data.

To simulate proper snow loads for the snow depth records described above, Δsnow was developed (Winkler et al., 2021a; Schellander and Winkler, 2020). It is a semi-empirical layer-resolving snow model which exclusively needs regular snow depth series as input to model respective series of snow loads. Winkler et al. (2021a) show that yearly snow load maxima can be modeled with a very low bias and a root mean square error of 0.36 kN m⁻², which is much more accurate than empirical regression models or using categorized or constant snow densities. The Δsnow model was applied to the daily snow depth series described above, leading to 897 snow load series used for spatial interpolation.

2.2 Spatial extreme value statistics

The characteristic snow load on the ground (s_k) can be calculated by assuming that yearly s_k maxima follow a generalized extreme value (GEV) distribution (Coles, 2001). Among many other environmental parameters, this has been already proven to be reasonable for snow depth (e.g., Marty and Blanchet, 2012; Nicolet et al., 2018). Given the strong correlation between snow depth and snow load, the GEV has also successfully been used for the estimation of extreme snow loads (Le Roux et al., 2020). However, there are also other choices like the lognormal distribution (Mo et al., 2016; DeBock et al., 2017).

Eurocode EN 1991-1-3 (CEN, 2015) demands a spatial representation of s_k. Different kriging variants were already used to interpolate snow load extremes for hazard mapping in Canada (Hong and Ye, 2014) and Croatia (Perčec Tadić et al., 2015). There exist, however, more appropriate methods for the spatial modeling of extremes, all trying to model the three parameters location (μ), scale (σ) and shape (ξ) of the GEV with spatially varying covariates. Another typical choice for modeling spatial extremes are max-stable processes (Haan, 1984), which were used by, for example, Schellander and Hell (2018) for snow depth. With max-stable processes also the spatial dependence of extremes has to be modeled. However, for the accuracy of the marginal distribution this has no noticeable effect (Schellander and Hell, 2018), discouraging max-stable processes for the simple case of mapping 50-year snow load return levels. Another approach which models GEV parameters with smooth functions in space (termed smooth spatial modeling, or SSM, hereafter) was already used to calculate snow depth extremes in Switzerland (Blanchet and Lehning, 2010) and Austria (Schellander and Hell, 2018). It is based on the maximization of the sum of the station-wise log likelihoods. As an example, snow water equivalents modeled with the Δsnow model were also spatially interpolated with a simple smooth model (see Appendix C in Winkler et al., 2021a).

In this paper we try to overcome known problems with the smooth modeling approach, where responses are modeled with a linear combination of parametric terms, which results in difficulties; e.g., a linear elevation dependence of the location and scale parameters leads to too large return levels in the mountains (Schellander and Hell, 2018; Winkler et al., 2021a). To that end we propose to model the three GEV parameters by fitting a generalized additive model (GAM; Hastie and Tibshirani, 1990) spatially to Δsnow-derived snow load maxima. GAMs are an extension of generalized linear models but use smooth functions rather than parametric terms for response estimation. They have already widely been used, e.g., for extreme wind speeds (Etienne et al., 2010) or daily rainfall (Sezer et al., 2016). We make three assumptions: (i) like in many applications (e.g., Blanchet and Lehning, 2010; Gaume et al., 2013; Schellander and Hell, 2018) the GEV parameters are modeled assuming a spatial non-stationarity in this study as well. For an easily reproducible snow load map, very deliberately, only topographical parameters (longitude, latitude and altitude) are used to explain the spatial variability in μ, σ and ξ. Adding more covariates (e.g., mean snow depth) might explain more spatial effects, leading to a better fit. However, it also results in a higher complexity of the model and the approach itself since all covariates have to be available on the possibly very fine grid on which the final map is to be computed. (ii) We also assume spatially independent snow load extremes. While this might not be the case in reality, Schellander and Hell (2018) showed that for the pure purpose of mapping marginals it does not have to be accounted for. (iii) It has been shown that 50-year snow load return levels in the northwestern French Alps have significantly decreased between 1960 and 2010 (Le Roux et al., 2020). However, in this methodological proposal we do not consider temporal non-stationarity of snow load extremes. As proposed by the forthcoming Eurocode (CEN, 2020, Annex A.5), we first adapt the current snow load map to current climatic conditions and shift consideration of temporal trends to subsequent studies.

A simple methodology was used to find the best suitable GAM for the purpose of mapping snow load return values. The key idea was to incrementally add complexity to a null model and find the point where additional complexity turns into overfitting, which was checked by RMSE and the Akaike information criterion (AIC; Akaike, 1974; see Appendix for details). For comparison a smooth spatial model was also fitted to Δsnow-derived snow load maxima; 50-year return levels of the generalized additive and smooth models were then compared with each other and with the current Austrian standard's map (ASI, 2018). As references, 50-year snow load return values, i.e., Δsnow-modeled yearly maxima locally fitted to a GEV distribution using maximum likelihood, were used (see Appendix for details).

Figure 4 RMSE (a) and bias (b) of 50-year snow loads for the different models and the Austrian standard (AS) classified by elevation. GAM stands for the generalized additive model, SSM for the smooth spatial model. The corresponding number of values within the elevation classes can be found in Fig. 2. Note that the current Austrian standard is not valid for elevations above 1500 m.

3 Results

Figure 3 shows 50-year return levels of Δsnow-modeled yearly maxima locally fitted to a GEV distribution for all 897 stations as described in Sect. 2. More than 75 % of the values are below 4 kN m⁻², with 50 % being smaller than 2.2 kN m⁻², especially in the lower elevations in the north, east and southeast. As expected, the highest values occur in the snowy northern and southern Alps, where topography enhances precipitation. Only eight stations between 1740 and 3109 m exhibit a 50-year snow load return level larger than 15 kN m⁻².

As outlined in Sect. 2.2 and explained in more detail in the Appendix, the best-suited GAM was found via a 10-fold cross-validation by fitting a GEV distribution to Δsnow-modeled maxima at the available 897 stations. In the notation of the R package mgcv (Wood, 2017) it is

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathit{\mu }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{120})+s(\mathrm{alt},k=\mathrm{9})\\ & \mathit{\sigma }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{120})+s(\mathrm{alt},k=\mathrm{9})\\ & \mathit{\xi }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{120})+s(\mathrm{alt},k=\mathrm{9}).\end{array}\end{array}$$

Thereby ∼ signals a formula; “maxima”, “long”, “lat” and “alt” refer to yearly snow load maxima, longitude, latitude and altitude at a station; “s” marks an isotropic smoother with a thin plate regression spline as a basis function; and “k” defines the dimension of the basis function, i.e., the degrees of freedom available for this term. Note that a possible improvement of the suggested model would mean a significant effort, which is associated with increased hardware resources. To test the goodness of fit of the best GAM (Eq. 1) the Anderson–Darling test statistics were used as suggested by Le Roux et al. (2020). Using the R package gnFit Saeb (2018), the test was applied to quantiles retrieved from the distribution modeled with Eq. (1) to assess whether they are likely to follow a GEV distribution at the 5 % significance level. This is true for 88 % of all used stations.

The best SSM for comparison was found via an AIC-based model selection procedure, which led to the following model for the three GEV parameters (see Appendix for details):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathit{\mu }\sim \mathrm{long}+\mathrm{lat}+\mathrm{alt}\\ & \mathit{\sigma }\sim \mathrm{long}+\mathrm{lat}+\mathrm{alt}+\mathit{\mu }\\ & \mathit{\xi }\sim \mathrm{long}+\mathrm{alt}+\mathit{\mu }+\mathit{\sigma }.\end{array}\end{array}$$

Table 1 shows RMSE and bias of the GAM, SSM and the current Austrian snow load map (ASI, 2018) with respect to 50-year return values of Δsnow-modeled yearly maxima locally fitted to a GEV distribution. The scores for the GAM and SSM were obtained by a leave-one-out cross-validation with the models (Eqs. 1 and 2).

Table 1 Error measures for 50-year return levels of snow load for the generalized additive model (GAM), smooth spatial model (SSM) and Austrian standard (AS), calculated using all 897 stations. (687 stations within Austria and below 1500 m for AS.) The scores are shifted in favor of the GAM if only stations within the Austrian border and below 1500 m are used. The RMSE value of the GAM is only 0.85 kN m⁻² then (not shown).

While both models and also the Austrian standard mostly overestimate reference snow loads, the GAM has the smallest RMSE and a very small, slightly positive bias. This is also confirmed in Fig. 4. Although the SSM uses more covariates, it is not able to model larger snow loads well. The GAM stays below an RMSE of 2.2 kN m⁻² up to 2000 m and rises to 6.9 kN m⁻² for elevations beyond 3000 m. This is accompanied by an increasing bias, which could at least partly be traced back to the small sample of only 14 stations in the highest elevation class (see Fig. 2). The newly proposed GAM approach shows smaller errors throughout all elevations when compared to the Austrian standard, even above 1500 m, but this must be considered with caution as the number of stations with elevations above 1500 m is only 57 (see Fig. 4).

Figure 5 Values of 50-year snow load return for a region around Windischgarsten in the center of Austria, taken from the map of the current Austrian snow load standard (a, AS; valid up to 1500 m) and the newly proposed GAM approach (b, plotted up to 2000 m). Pixel size is 50 m by 50 m. Main roads are depicted in the background for easier orientation.

The elevation class 1500 to 2000 m exhibits two severe negative outliers in the right panel of Fig. 4. Their locally estimated 50-year snow load return values are 14.1 and 15.2 kN m⁻². Both stations are situated in very snowy regions at elevations of 1905 and 1740 m, respectively. Modeling snow loads smoothly in space eliminates another shortcoming of current standards. As outlined in the introduction, European snow load standards rely on zonings. This drawback is outlined in the left panel of Fig. 5, which shows a region in the center of Austria around the village of Windischgarsten. Several strangely curved features can be noticed, resulting from sudden transitions between zones. This leads to jumps of snow load values. Furthermore, the current Austrian standards only allow for calculating s_k up to 1500 m elevation (white regions in the left panel of Fig. 5). The GAM approach (right panel in Fig. 5) is able to model smooth transitions while showing acceptable accuracy also at higher elevations (up to 2000 m in Fig. 5).

4 Conclusions and outlook

A transparent way towards a reproducible snow load map, conformable to CEN standards (CEN, 2015), is proposed: (i) simulating snow loads with the Δsnow model at intensively quality-checked, long-term snow depth records and (ii) using a generalized additive model (GAM) with penalized regression splines to spatially interpolate the generalized extreme value (GEV) distribution of the simulated snow load records. This procedure results in a map with smooth values of characteristic snow loads (s_k).

The presented methodology is proposed for the successor of the current Austrian snow load standard (ASI, 2018). The potential of the two-stage Δsnow – GAM approach could be shown. The RMSE and bias with respect to 50-year snow load references at 897 stations are much smaller than for smooth modeling or the current standard. However, for very high elevations above 2000 m the error measures reach unacceptably large values. Only longitude, latitude and altitude were used as explaining covariates. This not only ensures relatively low model complexity but also keeps the application of the method easy.

A similar approach, combining the methods of snow load simulation and appropriate spatial extreme value modeling, could be applied in other countries as well since long-term snow records are the only prerequisite. This would minimize inconsistencies of snow load standards at national borders, a target expressed by the European Committee for Standardization (CEN, 2015). The proposed approach avoids zoning and mitigates altitudinal restrictions. It provides a base to update snow load maps to current climatic conditions as suggested by the coming Eurocode (CEN, 2020, Annex A.5). On this fundament, finally, efforts of implementing possible climate-change-related trends of snow load extremes (like Croce et al., 2018; Le Roux et al., 2020) could be built on for a future generation of snow load standards.

Appendix A: Model selection for the generalized additive model

The R package mgcv (Wood, 2017) is used to fit a GEV spatially to snow load maxima. Smooth terms are represented by penalized regression splines, and the smoothing parameters are estimated during the fitting process. To find a suitable GAM for the spatial representation of snow load maxima, the following steps were performed: firstly a null model was defined in the notation of the R package mgcv, using the covariates longitude (long), latitude (lat) and altitude (alt):

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& \mathit{\mu }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{29})+s(\mathrm{alt},k=\mathrm{9})\\ & \mathit{\sigma }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{29})+s(\mathrm{alt},k=\mathrm{9})\\ & \mathit{\xi }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{29})+s(\mathrm{alt},k=\mathrm{9}).\end{array}\end{array}$$

Thereby “maxima” refers to yearly snow load maxima at a station, and the term “s” stands for a smooth term with a thin plate regression spline (Wood, 2003) as a basis function. Note that the intuitive standard, namely using one smooth term for each of the three covariates, leads to worse AIC and RMSE values. For the null model (Eq. A1), the dimension of the basis “k” is estimated during the fitting process. For the available 897 stations, the values were estimated to k=29 for the first term, s(long, lat), and k=9 for the second term, s(alt).

Secondly, the complexity of the GAM was incrementally increased for both terms by allowing more degrees of freedom, i.e., increasing k. Therefore the dimension k for the first term, s(long, lat), was set to 60, 120, 200 and 400 and in only one step to 20 for the second term, s(alt). The most complex model was then

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& \mathit{\mu }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{400})+s(\mathrm{alt},k=\mathrm{20})\\ & \mathit{\sigma }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{400})+s(\mathrm{alt},k=\mathrm{20})\\ & \mathit{\xi }:\mathrm{maxima}\sim s(\mathrm{long},\mathrm{lat},k=\mathrm{400})+s(\mathrm{alt},k=\mathrm{20}).\end{array}\end{array}$$

To find the best suitable GAM and to avoid overfitting, a 10-fold cross-validation for each model was performed. The dataset of 897 stations was randomly split into 10 subsets (Fig. A1). In 10 passes, 9 subsets were always used as a training dataset and the remaining one as testing dataset.

Figure A2 shows the evolution of the RMSE and AIC for training and testing datasets calculated from the first five models and the locally estimated 50-year return levels of Δsnow-derived snow loads using maximum likelihood for parameter estimation. The more complex models with k=20 for the s(alt) term led to larger RMSE values for both the training and testing datasets and are not shown.

With increasing k, i.e., degrees of freedom, the RMSE of both the training and testing datasets always decreases. But for k>120 the testing dataset no longer benefits from increasing complexity, which is the point where overfitting sets in. As expected, the AIC also decreases (numbers in Fig. A2). Therefore a model (Eq. A1) with k=120 for the first term and k=9 for the second term was selected as the best GAM.

Figure A1 Elevation histograms of the 10 randomly found subsets of all 897 stations. Respective number of members is given in brackets.

Figure A2 Lines show the RMSE based on a 10-fold cross-validation of all 897 stations for testing and training datasets for models with increasing degrees of freedom. The numbers refer to the largest AIC of the training datasets out of 10 passes for each degree of freedom.

Appendix B: Model selection for the smooth spatial model

To fit a smooth spatial model (SSM) to the snow load maxima we followed the approach of Gstöhl (2017) as a blueprint. As covariates for the SSM, longitude, latitude and altitude are used, but the location and scale parameter are also allowed as explaining covariates. To select the best regression formula for each of the three GEV parameters, a generalized linear regression between the GEV parameter (e.g., μ) and covariates added stepwise was performed. The formula between a null model (e.g., μ∼1) and a given full model (e.g., location parameter μ ∼ all covariates) leading to the smallest AIC was selected.