The annual exceedance probability of the once-in-M years event is 1/M. Thereby, M often is referred to as the return period. It is important to note that this term may lead to the impression that there is a regular basis for the occurrence of this event. As can be seen on the example of table 2, which summarizes some probabilities for different ranges of inter-arrival times for the once-in-fifty years event, the inter-arrival times show a large random scatter. As consequence, calibration of structural design rules against existing successful design practice is a kind of self-deception.

Figure 4 shows a typical time series of the random exceedances of the once-in-fifty-years event for 40 consecutive 50 year periods. The largest observed inter-arrival interval in this example is 215 years.

<table><caption>**Table 2:** Examples of inter-arrival times for exceedance of the once-in-fifty years event</caption><tr><td>T</td><td>&lt; 10</td><td>&lt; 20</td><td>&gt; 50</td><td>&gt; 100</td><td>&gt; 200</td></tr><tr><td>f(T)</td><td>0.181</td><td>0.330</td><td>0.368</td><td>0.135</td><td>0.018</td></tr></table>

<figure><img src="image_94.png" /><figcaption>**Figure 4:** Typical time series of the random exceedance of the once-in-fifty-years event for 40 consecutive 50 years periods</figcaption></figure>

In a general way reliability is defined as the probability that a system or component performs its required functions under stated conditions for a specified period of time. In structural engineering, system or component mean the structure or structural elements. The required functions are structural safety, serviceability, durability and robustness. The stated conditions are the actions which are occurring during the construction and anticipated use. Finally, the specified period is the design or projected working life.

Structural safety is the ability of a system or component to withstand all actions as planned. Serviceability is the ability of a system or component to remain fit for the use it is required. These two demands are considered with the Ultimate Limit State (ULS) and the Serviceability Limit State (SLS). Durability is the ability of a system or component to guarantee the structural safety and serviceability during the specified lifetime. Robustness is the ability of a system or component to withstand accidental events like fire, explosions, impact or consequences of human error, with damages which are not disproportional to the triggering event.

The complementary to reliability is the failure probability. For R being the capacity or the resistance of a system and E being the stress or action effect, the failure probability can be obtained from the convolution of the two contributing probability densities. Failure occurs if E becomes larger than R, i.e. the limit state function for a linear system is simply R-E = 0. The failure probability p<sub>f</sub> then is obtained by the so-called failure integral as:

$$p_f = \int_{R=0}^{\infty} f_R(R) \cdot \int_{E=R}^{\infty} f_E(E) \, dE \, dR \quad (10)$$

with f<sub>R</sub> and f<sub>E</sub> - probability densities of the resistance and the action effect

The second integration in equation (10) can be replaced by a term including the cumulative probability distribution of E:

$$p_f = \int_{R=0}^{\infty} f_R(R) \cdot [1 - F_E(R)] \, dR \quad (11)$$

with F<sub>E</sub> - cumulative probability distribution of E

It is important to note that the failure integral leads to the average failure probability. The probability distribution of the failure rate usually shows large positive skewness.

Not all of today's Standards do allow considering the projected working life as an explicit design variable; instead, they often consider implicitly a general value of 50 years. Typical values for the projected working life are 1 - 5 years for temporary structures, 10 - 40 years for industrial structures and buildings, 60 - 80 years for residential buildings and 100 - 150 years for bridges.
Action effects are internal forces like bending moments or normal forces, or stresses from a combination of internal forces, or deflections or any other structural response which has to be considered in the design process. In the general case, several actions with different origin will occur simultaneously and will test the designed structure in regard to a sufficient resistance. Beside snow actions, usually at least actions induced by dead load and wind have to be considered. It is important to note that a consistent evaluation of the reliability requires knowledge of the static system and basic structural behaviour.

Since the solution of the failure integral is too complex for the every-day engineering practice, Standards have introduced the so-called operative failure probability and the corresponding reliability index. The latter is obtained as the corresponding negative value of the reduced variate of the inverse of the standardized normal distribution Φ, i.e.:

$$ 
\begin{aligned}
p_f &= \Phi(-\beta) \\
\beta &= \Phi^{-1}(1-p_f)
\end{aligned} 
\qquad (12) $$

In a simplified approach the design target is met by comparing the design value of the resistance R to the design value of the action effect E in the following inequation:

$$ E_{des} \le R_{des} \qquad (13) $$

If E has no physical upper limit, failure may occur for any given level of R, and the major accumulation of the failure probability extends over a large range of R and/or E, respectively. In other words: failure may occur for events E which are below the design value of the action effects.

The design values R<sub>des</sub> and E<sub>des</sub> can be obtained from statistics for the respective random variables if corresponding target values for the exceedance probabilities are specified. Target values of the reliability can be specified in dependence of the consequences of failure. In the Eurocode EN-1990, three consequence classes (CC) are distinguished. The complete set of the European Action Standards EN-1991 basically aims to provide the required information for CC2 which is described with medium consequences for loss of human life; economic, social or environmental consequences are considerable. Typical examples are residential and office buildings and public buildings where consequences of failure are medium. CC3 covers failures with high consequences for loss of human life, or economic, social or environmental consequences are very great. Bridges, stadia and concert halls belong to this class. The lowest consequence class CC1 is characterized by low consequences for loss of human life, and economic, social or environmental consequences which are small or negligible. A typical example is given with commercial greenhouses where normally no person enters.

The target values for the reliability associated with the ULS are specified in EN 1990 as minimum values of the

reliability index β with reference to 50 years and 1 year. They are summarized in table 3. Additionally in table 4, the corresponding approximate target failure probabilities are given.

<table><caption>**Table 3:** Target reliability values in EN 1990</caption><thead><tr><th>Reliability Class</th><th>1 year</th><th>50 years</th></tr></thead><tbody><tr><td>RC3</td><td>5.2</td><td>4.3</td></tr><tr><td>RC2</td><td>4.7</td><td>3.8</td></tr><tr><td>RC1</td><td>4.2</td><td>3.3</td></tr></tbody></table>

<table><caption>**Table 4:** Approximate target failure probabilities in EN 1990</caption><thead><tr><th>Reliability Class</th><th>1 year</th><th>50 years</th></tr></thead><tbody><tr><td>RC3</td><td>1/10,000,000</td><td>1/100,000</td></tr><tr><td>RC2</td><td>1/1,000,000</td><td>1/10,000</td></tr><tr><td>RC1</td><td>1/75,000</td><td>1/2,000</td></tr></tbody></table>

In a simplified approach, EN-1990 specifies the design values for the resistance and the action effects based on a unique pair of coordinates which are

$$ R_{des} = -0.8 \cdot \beta; \quad E_{des} = 0.7 \cdot \beta \qquad (14) $$

For CC2, the design value of the resistance corresponds to the 0.1%-fractile value, i.e.:

$$ p(R \le R_{des}) = 0.001 \qquad (15) $$

With reference to a single year, the design value for the action effect has an exceedance probability of 1/2000:

$$ p(E \le E_{des}) = 0.0005 / \text{year} \qquad (16) $$

Instead of specifying this design values directly, most Standards still use the concept of partial factors, i.e. the target value R<sub>des</sub> will be replaced by the ratio of a characteristic value R<sub>k</sub> and a corresponding partial factor γ<sub>M</sub>, i.e. the design value is obtained as:

$$ R_{des} = \frac{R_k}{\gamma_M} \qquad (17) $$

Different materials have different partial factors, thus reflecting the different levels of scatter on the resistance side. Basically, the larger the scatter, the large must be the partial factor.

For the design value of the action effects the corresponding equation is:

$$ E_{des} = \gamma_E \cdot E_k \qquad (18) $$

with γ<sub>E</sub> – partial factor for the action effect
The choice of the characteristic value is free, since for every value of e.g. E<sub>k</sub> there is a corresponding value γ<sub>E</sub> which leads to the target value E<sub>des</sub>. In recent years, most Standards, however, prefer to choose an annual exceedance probability of 1/50, thus suggesting that the return period is coupled with the implicitly assumed design working life of 50 years. In the German annex to EN-1990, the partial factor for unfavourable variable loads is uniformly set to 1.5, thus suggesting that the wind and the snow climate show the same scatter all over Germany and furthermore equal the scatter of imposed loads.

A more transparent and consistent approach specifies explicitly the target probabilities for the design values of the actions considering different structural classes, e.g. ISO 22111 (2007) Bases of design of structures – General requirements [1] or ISO 4354 (2009) Wind actions on structures [2].

A refined classification of the consequence classes and the importance level of buildings and structures is given in table 5. The corresponding tentative target values for the exceedance probabilities of the design wind and snow loads are given in table 6 with reference to the projected working life and in terms of an annual exceedance probability for the example of L = 50 years. Basically, the required annual exceedance probabilities for any other projected working life can be obtained from the inverse of equation (5) as:

$$p_{\text{year}} = 1 - (1-p_{\text{lifetime}})^{1/L} \quad (19)$$

with L – projected working life in years

<table><caption>**Table 5:** Refined classification of consequence classes</caption><thead><tr><td>Class</td><td>Building or structure type</td></tr></thead><tbody><tr><td>0</td><td>Buildings or structures presenting no hazard to life and other property in the case of failure, e.g. free field solar parks</td></tr><tr><td>1</td><td>Buildings or structures presenting a low degree of hazard to life and other property in the case of failure, e.g. agricultural buildings without regular human occupancy, greenhouses without public access</td></tr><tr><td>2</td><td>Buildings or structures not included in class 0, 1, 3 and 4, 5, e.g. office buildings, commercial buildings, factories, residential buildings</td></tr><tr><td>3</td><td>Buildings or structures that are designed to contain a large number of people, e.g. high-rise buildings, stadia and concert halls</td></tr><tr><td>4</td><td>Buildings or structures that are essential to post-disaster recovery or associated with hazardous facilities, e.g. hospitals, bridges</td></tr><tr><td>5</td><td>Buildings or structures that are essential to safe the lives of people in case of extreme natural hazards, or have inestimable cultural value, e.g. shelter buildings, monuments</td></tr></tbody></table>

<table><caption>**Table 6:** Tentative values for the exceedance probabilities of the design wind and snow loads</caption><thead><tr><td>Class</td><td>p<sub>lifetime</sub></td><td>p<sub>year</sub> (L=50 years)</td></tr></thead><tbody><tr><td>0</td><td>0.40</td><td>1:100</td></tr><tr><td>1</td><td>0.20</td><td>1:200</td></tr><tr><td>2</td><td>0.10</td><td>1:500</td></tr><tr><td>3</td><td>0.05</td><td>1:1000</td></tr><tr><td>4</td><td>0.025</td><td>1:2000</td></tr><tr><td>5</td><td>0.01</td><td>1:5000</td></tr></tbody></table>

### 3. Simulation as design tool

The core module of any simulation is the random number generator. In 1970, Robert Coveyou created the famous words that *the generation of random numbers is too important to be left to chance*. Strictly speaking the numerically generated results have to be called pseudorandom numbers since they are generated by deterministic algorithms. If the starting or seed value for the procedure remains the same, a once produced sequence can be generated again and again. This is the basic difference to a true sequence of random numbers.

The basic demand to the random numbers is that they are uniformly distributed in the range from 0 to 1. Then, the random numbers can be interpreted as probabilities, and can be translated to any target probability distribution e.g. by the inverse transformation method as illustrated in figure 5.

<figure><img src="image_95.png" /><figcaption>**Figure 5:** Example of using random numbers r<sub>i</sub> to generate random values x<sub>i</sub> = F<sup>-1</sup>(r<sub>i</sub>)</figcaption></figure>

The basic demands for the generated random numbers are:

* The period length of the sequence should be large enough to avoid any repetition of the sequence in the actual simulation problem.

* There should be no internal correlation in the sequence.

* The demand to a uniform distribution is met for n dimensions, e.g. for two dimensions it is required that x<sub>i-1</sub> is uniformly distributed for whatever small range is specified for x<sub>i</sub>.
A bad random number generator can ruin a simulation. An illustrative example is shown in figure 6, applying the routine ran1 from Numerical Recipes [3]. In the version of 1992, the authors state that in fact we do not know of any statistical test that ran1 fails to pass, except when the number of calls starts to become on the order of the period, say $10^8$. The example generates for a complete year 365.25·24 = 8766 hourly mean wind speeds following a Weibull distribution which is given as:

$$F_v(v) = 1 - \exp\left(-\left(\frac{v}{v_0}\right)^k\right) \qquad (20)$$

with $v_0 = 4.36$ m/s and $k = 1.78$

<figure><img src="image_96.png" /><figcaption>**Figure 6:** Example of applying a bad random number generator</figcaption></figure>

From each set, the maximum wind speed is sampled together with the hourly mean value before and after the extreme. Altogether, $10^7$ runs are performed. The first feature to notice is that the largest generated value is fixed with a value of 21.157 m/s thus indicating that the largest value being smaller than 1 is $1-1/2^4$. Since the demand of uniform distribution for more than one dimension is not met, the wind speeds in the hour before and after the maximum show strange probability distributions. Basically, the neighbouring hours should have the same probability distribution. Since the values are limited by the actual maximum value, the distribution

should show an under-linear characteristic; however, they should increase in a strict monotonic way. The bad ran1-generator leads to fixed values for the neighbouring hours, thus ruining the simulation completely.

Much better results are obtained with the class of Mersenne Twister random number generators. Their names derive from the length of the period which is given by a Mersenne prime. The first version has been published in 1997 by Matsumoto and Nishimura [4] and provides the unbelievable period of $2^{19937}-1$. The demand to uniform distribution is met for 623 dimensions. The code is available for free in the internet. The largest number smaller than 1 is $1-1/(2^{32}-1)$. The improved version SFMT [5] is faster than the original version, the dimension of uniform distributions is better and the period is pushed even further to $2^{216091}-1$. The largest number below 1 is $1-1/(2^{35}-1)$. Another high-quality random number generator [6] has been introduced as KISS (Keep It Simple Stupid). The 64 bit version has a period of $2^{121}+2^{63}-1$. The largest number below 1 is $1-1/(2^{32}-1)$.

Basically, the largest number below 1 becomes important and may cause problems when the target probability distribution extends to infinity and when the ensembles become large. Then, all the simulations end up with the same largest number. A typical example is obtained for traffic-induced loads for bridges which often are assumed to follow a normal distribution. Today's design scenarios [7] consider 10,000 trucks per day and 200 busy days per year. Each truck has 5 axles. With a projected working life of 100 years, the required number of random weights is $5 \cdot 10,000 \cdot 200 \cdot 100 = 10^9$. The simulation with the KISS-algorithm leads to only $2^{32} = 4.295 \cdot 10^9$ different random numbers, i.e. repeated runs over the projected working life produce on average in 1 of 4.295 runs the same largest number, which in terms of the reduced variate of the standardized normal distribution is 6.2303.

A refined approach fills each gap between the original integer numbers between $-2^{31}$ and $2^{31}-1$ with $2^{21}$ uniformly distributed values, thus ending up with $2^{53}$ random numbers which matches a double precision floating point number. The final random number is obtained from two calls of the KISS-random number generator:

$$r = (r_{1,KISS} + 2^{31} + \frac{r_{2,KISS}/2^{11} + 2^{20}}{2^{21}-1})/(2^{32}-1) \qquad (21)$$

with r uniformly distributed in [0, 1]

The largest value smaller 1 then becomes $1-1/(2^{53}-1)$, which for a normal distribution corresponds to a reduced variate of 8.2095. If the gaps are filled with the complete available range of $2^{32}$ numbers, a quadruple precision format is required thus extending the range of numbers smaller than 1 to $1-1/(2^{64}-1)$. The corresponding reduced variate of the standard normal distribution becomes 9.080.

Finally, a simulation is only as good as are the assumed models to describe the contributing processes. An inconsistent model may ruin the simulation results as may a bad random number generator.
## 4. Probability models for the snow load on the ground

There are two options to analyse the extreme values of the observed snow loads. The straight forward method is based on annual extremes; the corresponding theoretical probability distribution is the generalized extreme value distribution (GEVD) which is given as:

$$F(x) = \exp \left[ - \left( f_1 - \text{sign}(\tau) \cdot f_2 \cdot \frac{x-m}{\sigma} \right)^{1/\tau} \right] \quad (22)$$

m - mean value,

σ - standard deviation,

τ - shape parameter

The coefficients $f_1$ and $f_2$ depend on the actual shape parameter $\tau$ and are given as follows:

$$f_1 = \Gamma(1+\tau), \quad f_2 = \sqrt{\Gamma(1+2\cdot\tau) - f_1^2} \quad (23)$$

$\Gamma$ - Gamma function

For positive shape parameters, type III distributions are obtained. This family of distributions has a finite upper tail, which cannot be exceeded with a probability of 1. The maximum value is given as:

$$x_{\max} = m + \sigma \cdot \frac{f_1}{f_2} \quad (24)$$

Negative shape parameters lead to the family of the type II distributions. These distributions have a finite lower tail. For $\tau = 0$, the extreme distribution type I is obtained, i.e. in the limit equation (22) becomes:

$$F(x) = \exp \left[ -\exp \left( - \left[ \gamma + \frac{\pi}{\sqrt{6}} \frac{x-m}{\sigma} \right] \right) \right] \quad (25)$$

$\gamma$ - Euler constant = 0.577216

The distributions of type II and III each form a family of curves with specific characters. Compared to the type I distribution, they show as special feature a certain curvature when plotted in Gumbel probability paper (figure 7). Generally, a probability paper is a graph paper with one axis specially ruled to transform the distribution function of a specified theoretical expression to a straight line when it is plotted against the variate as the abscissa. While the curves for the type II distribution bend in a concave shape in respect to the axis of the reduced variate, the curves corresponding to type III show a distinct convex character. These two types are separated by the type I distribution which appears in the plot as a straight line.

Some studies use the normal or log-normal distribution to describe the trace of the non-exceedance probability of annual extremes. In figure 8, the respective

traces are shown. Obviously, these approaches can be covered with corresponding type III or type II distributions.

The alternative approach is based on extremes over a threshold value. This strategy makes sense in snow climates which do not lead to a continuous snow cover through the whole winter period, and/or do not have snow in each winter. As additional parameter, the number of independent events per twelve-month cycle has to be considered. The appropriate theoretical expression for the non-exceedance probability per event is given with the generalized Pareto distribution (GPD) as:

$$F(x) = 1 - \left[ 1 - \frac{x - x_s}{s} \cdot k \right]^{1/k} \quad (26)$$

$x_s$ - threshold value

s - scale parameter

k - shape parameter

In case of a positive shape factor, the distribution has a finite upper tail which cannot be exceeded. The corresponding largest value is given as:

$$x_{\max} = x_s + s/k \quad (27)$$

For $k = 0$, the exponential distribution is obtained:

$$F(x) = 1 - \exp \left( - \frac{(x - x_s)}{s} \right) \quad (28)$$

The annual non-exceedance probabilities are obtained from combining the probabilities of the number of events and non-exceedance probabilities as follows:

$$p(x \le x_{\text{ref}} | \text{year}) = \sum_{N=0}^{\infty} p(N) \cdot p(x \le x_{\text{ref}})^N \quad (29)$$

If the number of events per twelve month follows the Poisson distribution, the trace from equation (29) can be expressed by the GEVD as given in equation (22). The Poisson distribution is given as:

$$p(K) = \frac{\lambda^K}{K!} \exp(-\lambda) \quad (30)$$

with $\lambda$ - average number of events per 12 month cycle

For very rare events, i.e. for $p(x \le x_{\text{ref}})$ close to unity, equation (29) can be simplified to:

$$p(x \le x_{\text{ref}} | \text{year}) = p(x \le x_{\text{ref}})^{\lambda} \quad (31)$$

Translating the GPD in combination with the Poisson distribution to the corresponding annual GEVD leads to the same exponents:

$$k = \tau \quad (32)$$
The basic relation between the scale in the GPD and the standard deviation in the GEVD is given as:

$$ \sigma = (c_0 + c_1 \cdot \ln \lambda + c_2 \cdot \ln^2 \lambda) \cdot sc \quad (33) $$

The coefficients $c_0$, $c_1$ and $c_2$ depend on the actual $\tau$-values. For positive $\tau$-values the following approximations apply:

$$ \begin{aligned} c_0 &= 1 - 0.1 \cdot \tau - 0.1 \cdot \exp(1) \cdot \tau^2 \\ c_1 &= 0.02672 - 1.07426 \cdot \tau + 0.25264 \cdot \tau^2 \\ c_2 &= -0.05419 \cdot \tau + 0.6 \cdot \tau^2 \end{aligned} \quad (34) $$

For negative $\tau$-values, the coefficients become:

$$ \begin{aligned} \ln c_0 &= 0.25242 - 1.10219 \cdot \tau + 3.3628 \cdot \tau^2 \\ c_1 &= -0.78782 \cdot \tau + 5.15605 \cdot \tau^2 \\ \ln c_2 &= -7.9336 - 31.11332 \cdot \tau + 38.14 \cdot \tau^2 \end{aligned} \quad (35) $$

The corresponding mean value is obtained as:

$$ m = x_{th} + \frac{sc}{\tau} - \sigma \cdot \frac{f_1}{f_2} \quad (36) $$

The limiting case $\tau = 0$ leads to:

$$ \begin{aligned} m &= x_{th} + (\gamma + \ln \lambda) \cdot sc \\ \sigma &= \frac{\pi}{\sqrt{6}} \cdot sc \end{aligned} \quad (37) $$

<figure><img src="image_97.png" /><figcaption>**Figure 7:** Traces of extreme value distributions in Gumbel probability paper</figcaption></figure>

<figure><img src="image_98.png" /><figcaption>**Figure 8:** Traces of the normal and log-normal distribution in Gumbel probability paper
1 - normal distribution, 2 - 5 log-normal distribution with V = 0.1, 0.3, 0.5 and 1.0, 0 - extreme value distribution type I</figcaption></figure>

For the estimation of the parameter triple (m, $\sigma$, $\tau$) or ($\lambda$, sc, k), different methods can be applied:

* the method of moments,

* the maximum likelihood method,

* the least-square approach

* the method of L-moments

* the Best Unbiased Linear Estimator

For some of the above methods, the plotting position of the individual observations is required which is estimated as:

$$ p_i (x \le x_i) = \frac{i \cdot \alpha}{N + (1 - 2 \cdot \alpha)} \quad (38) $$

with i - position of the observation $x_i$ in the ascending sorted list

Typical values for $\alpha$ are in the range from 0 to 0.5. In Gumbel's original concept [8], $\alpha$ is set to zero. Later, Gringorten [9] recommended that for plots in the Gumbel probability paper $\alpha = 0.44$ should be used. Basically, the estimated values for the design snow load will decrease if $\alpha$ increases.

For a confined ensemble, each of the above methods leads to different values of the describing parameters. The statistical stability of the identified parameters remains poor even for large observation periods of several decades. This is shown in the following on the
example of observed snow loads at Dresden, which has been published already in [10].

Figure 9 shows the trace of the observed non-exceedance probabilities for 59 independent events in 28 winter periods at Dresden and the identified theoretical distribution. The corresponding parameters are $sc = 0.2105 \text{ kN/m}^2$, $k = 0.163$ and $\lambda = 2.107$. The identified curve seems to be a quiet convincing fit of the observed trace, i.e. the identified GPD seems to be an appropriate model. The naïve design value for a projected working life of 50 years becomes $0.975 \text{ kN/m}^2$ for a class 3 structure.

<figure><img src="image_99.png" /><figcaption>**Figure 9:** Trace of the observed non-exceedance probability of the snow load on the ground and identified theoretical distribution (Dresden 1978/79 - 2005/06)</figcaption></figure>

The analysis of the statistical stability of the identified parameter-triple ($\lambda$, sc, k) is based on simulations. In figure 10, the 90%-confidence intervals are shown. The example illustrates that the naïve estimation of the design value purely based on the identified statistical parameters is not recommended. Obviously, an ensemble size of even 100 years is not large enough to specify the design value in a sufficient small two-sided confidence interval. Therefore, EN 1990 recommends applying a one-sided confidence interval to avoid an overestimation of the resistance and an underestimation of the action effects. Since large confidence values lead to over-conservative design, as compromise between safety and economy the 75% confidence is chosen.

The problem of confined ensembles also effects the estimation of the design values for wind and temperature. The question if the extremes follow type I, II or III often can hardly be answered. A pragmatic approach therefore assumes that extremes follow the extreme value distribution type I.

If the applied estimation method does lead to the desired confidence level, an adjusting factor is required. A consistent approach is discussed in [10].

<figure><img src="image_100.png" /><figcaption>**Figure 10:** Influence of the ensemble size on the statistical stability of the estimated parameter triple ($\lambda$, sc, k) [9]</figcaption></figure>

## 5. Outliers

The general definition defines an outlier as an observed value which lies considerably outside the range of the other observations. Outliers may occur by chance or may indicate erroneous observations. The latter can be removed from the ensemble. For the randomly observed outlier, a consistent and generally agreed strategy is still missing. An outlier may influence or distort the estimation of the describing parameters of the assumed probability distribution, thus leading to biased design values. This is illustrated in figure 11 for the extreme value distribution type I varying the ensemble size from 8
to 64 and considering exceedance probabilities 1/R for the randomly observed event from R = 100 to 2000.

The fitting procedure for the (N-1) values leads to a straight line for the identified theoretical distribution, i.e. the true observations can be replaced by equidistant points along the straight line. Then, it is easy to calculate the bias due to the outlier on the estimated design value which here is chosen with p (x ≤ xdes) = 0.999. The smaller the ensemble and the larger the return period R of the outlier, the larger will be the distorting effect of the outlier.

<figure><img src="image_101.png" /><figcaption>Figure 11: Relative position of outliers to the lsq-fitted trace for the reduced sample</figcaption></figure>

<figure><img src="image_102.png" /><figcaption>Figure 12: Biasing factor for the estimated design value with p(x ≤ xdes) = 0.999 / per year</figcaption></figure>

Unfortunately, the real world problem is a bit more complex due to the fact that for a confined ensemble the detection of an outlier is not always successful since the true describing parameters may not be correctly identified. Basically, the evaluation of an ensemble in regard to an outlier results in one of four cases:

* true 1: the ensemble contains no outlier and the detection method identifies no outlier

* true 2: the ensemble contains an outlier and the detection method identifies an outlier

* error 1: the ensemble contains no outlier but the detection method identifies the largest observation as outlier

* error 2: the ensemble contains an outlier but the detection method fails to identify the largest observation as outlier

In simulations, it is possible to distinguish these four classes; however, in a real world problem, the results are only sorted in two classes, namely:

* no: the detection method has not identified an outlier

* yes: the detection method has identified an outlier

Based on simulations, for each of the four cases, the corresponding probabilities can be obtained. The final evaluation is based on the probability of case yes which is the sum of true 2 and error 1.

The following example illustrates the uncertainties in the outlier detection based on an ensemble size of N = 32 observations. The return period is chosen as R = 500. Then, on average, an outlier is expected in 32 observations with a probability of (1- (1-1/500)³²) = 6.2%. The mean value and the standard deviation are obtained from a least-square fit in the probability paper applying the plotting position i/(N+1). It is important to note that this agrees with the strategy of identifying exceptional snow fall events in Europe. The outlier is identified based on the so-called k-method:

$$x_{\text{out}} > m_{N-1} + k \cdot \sigma_{N-1} \quad (39)$$

with mN-1 and σN-1 – statistical parameters for the ensemble without the largest value

The weighting factor k is obtained as the reduced variate of the respective probability distribution for the specified target exceedance probability 1/R of an outlier. For R = 500, the weighting factor k becomes 4.395 for the Gumbel distribution, which is the inverse of equation (25) for F(x) = 0.998.

Since the chosen plotting position overestimates the standard deviation, the detection probability for the case true 2 is only 2.23%, i.e. only about 1 of 3 true outliers is detected. Erroneous outliers are obtained with a probability of 0.74%, i.e. about 1 of 4 detected outliers is not really an outlier. Increasing the return periods R for the definition of an outlier leads to a decreasing detection rate for true outliers and an increasing rate of erroneous outliers.

The main conclusion at this stage has to be that the exceptional snow loads in the Eurocode are not based on scientific sound statistics. Figure 13 shows the probability traces for randomly estimated characteristic values with pannual (x > xk = 1/50) and design values with pannual (x > xk = 1/1000) for the four classes true1, true2, error1, error2 and for the resulting combinations no and yes. If erroneously an outlier is removed, the estimated characteristic and design values are far on the unsafe side, as can be seen with the traces for error 1. Failing to identify an outlier on the other hand leads to drastic overestimations (error 2). However, the observed trace does not reveal if or if not the detected outlier is a true outlier. Therefore, the final evaluation has to be based on the traces for the cases that the detection method identifies the largest value to be an outlier (yes) or not
(no). For the case that no outlier is detected, the traces reveal that the characteristic value is under-estimated with a probability of 30.8%. The estimation of the design value leads to an error probability of 29.7%. Both values are not that far from the EN-recommended target-value of 25%.

However, if the detection method identifies the largest value as an outlier, the error probabilities increase drastically to 52.9% for the characteristic value and

51.5% for the design value. These high levels of under-design are not acceptable. For all stations, where exceptional snow loads have occurred, the estimated characteristic values are on the unsafe side with a probability of about 50%. It is strongly recommended to increase the respective snow loads as soon as possible applying an appropriate adjusting factor for the respective ensemble sizes.

<figure><img src="image_103.png" /><figcaption>**Figure 13:** Traces of the estimated characteristic and design value for different error-cases and their combinations</figcaption></figure>

The question arise, how the detection strategy can be improved. Changing the plotting position to the recommended value by Gringorten increases the detection rate to about 45%, the rate of false outliers

increases to 33%. The method of moments leads to a detection rate of true outliers of 38.6%; however, the rate of wrong outliers is increased by a factor of about 18 compared to the initial strategy, i.e. for 10% of the
random ensembles, an outlier is identified with about only 1 of 4 being a true outlier.

For the application of BLUE the problem has to be solved how to obtain the weighting coefficients for a censored ensemble. In [13] tables are presented, however only for ensemble size up to 30. In [14] a series approximation is proposed for the calculation of the BLUE-weighting coefficients for larger ensemble sizes. This routine can be re-written to obtain weighting coefficients for censored ensembles as well. The detection rate with BLUE increases to 69%, the rate of false detections is 42%. It is important to note that both cases yes and no require an adjusting factor for the estimation of the design values to meet the 75% confidence level.

## 6. The influence of wind

Basically, the wind speed during a snow fall event will influence the amount of snow on the roof. The stronger the wind blows, the lesser snow will accumulate on the roof. During and after the snow fall event, wind may lead to a re-distribution of the snow masses on the roof. A probabilistic model for the snow on gable roofs is presented in the JCSS probabilistic model code [15], specifying probability distributions for the shape factor μ and the re-distribution coefficient c_r which depend on the mean wind speed. The roof snow loads for the upwind and downwind part of the roof is specified as follows:

$$s_1 = (\eta - c_r) \cdot s_{gr} \qquad (40)$$

$$s_2 = (\eta + c_r) \cdot s_{gr}$$

In [16], a refined proposal for the design values of the shape factor and the re-distribution coefficient is presented. The basic strategy convolves the probability densities of the shape factor, the re-distribution coefficient and the snow loads on the ground. The identified values for the shape factor explicitly depend on the mean wind speed at roof height. Basically, beside the non-uniform snow load distribution it is always advised to check the balanced snow load with $s_1 = s_2 = \mu \cdot s_{gr}$ as well. Figure 14 summarizes the results for the assumption that the snow loads follow an extreme value distribution type I.

The actual EN-document on snow actions recommends a shape-factor of 0.8. It is important to note that this value requires mean wind speeds in the range from 5 to 10 m/s.

For the final design of a structure, the combined action effects have to be considered. Today's standards follow the pragmatic strategy of specifying accompanying load amplitudes, which are modelled by combination factors, i.e. the two design scenarios are full wind load plus a reduced snow load and full snow load plus a reduced wind load. In figure 15, the trace of the non-exceedance probabilities of the largest hourly wind speed in individual storms is plotted for Dresden. The average number of storms per twelve month cycle is 1.19,

the fit of the GPD leads to a threshold value of $vt_h = 15.1$ m/s, a scale value of $v_{sc} = 2.170$ m/s and a shape parameter of 0.2092.

<figure><img src="image_104.png" /><figcaption>Figure 14: Design values for the snow loads on gable roofs [15]</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>Figure 15: Probabilistic model for hourly wind speeds in storms (Dresden: 1973 – 2007)</figcaption></figure>

The analysis of the wind-snow climate has to be based on the observed snow heights which are available with one value per day. The snow loads are only available for
three times per week. In a simplified approach, all 24 hourly mean wind speeds are assigned to the height of snow cover of that specific day.

In figure 16, the wind-snow climate is analysed for Dresden. For each level of snow on the ground, the mean value of all observed hourly wind speeds is shown. Additionally, the largest observed value is shown. Obviously, the probability of getting a storm when there is snow reduces with increasing snow height. The mean value of the wind speed decreases for increasing snow height. Fitting a linear trend allows evaluating the expected wind conditions for the level of characteristic and the design snow loads. Therefore, in figure 16 tentative levels are marked for the characteristic and the design snow load. For the design snow level, the mean value of the expected wind speed will be in the range of 2 m/s. The corresponding velocity pressure is 0.0025 kN/m², assuming the air density is 1.25 kg/m³. This wind load is small enough to justify that the design snow load requires no accompanying wind load, i.e. the respective combination factor is 0.

<figure><img src="image_106.png" /><figcaption>**Figure 16:** Wind conditions for different levels of the snow height on the ground (Dresden 1973-2007)</figcaption></figure>

For the analysis of the joint-probability of storm wind speeds and accompanying snow only 17 observations are available, which does not allow developing a sophisticated model for the joint-probability of snow and wind. The observations however suggest that the larger the wind speed, the smaller will be the snow height. The largest observed wind speed with 2 cm snow on the ground is 17.4 m/s.

A refined approach tries to extrapolate from the basic population of the observed wind speeds for different snow height levels. In figure 17, the histogram of the observed wind speeds is shown for a snow height of 9 cm. Obviously, the histogram can be fitted to a Weibull distribution which is given as:

$$f(v) = \frac{k}{v_0} \cdot \left(\frac{v}{v_0}\right)^{k-1} \cdot \exp\left(-\left(\frac{v}{v_0}\right)^k\right) \quad (41)$$

with $v_0$ – scale parameter, k – shape parameter

The fitting procedure is based on a least-square fit of the cumulative probability distribution. This gives larger

weight to higher values than fitting the density. Figure 18 shows the dependency of the two describing parameters $v_0$ and k on the height of the snow cover for a range up to 25 cm. Beyond this value, the number of observed wind speeds is no longer large enough to obtain stable estimates of the Weibull distribution.

<figure><img src="image_107.png" /><figcaption>**Figure 17:** Histogram of the observed hourly wind speeds for a snow height of 9 cm and fitted Weibull distribution (Dresden 1978-2007: $v_0 = 4.87$ m/s, $k = 1.73$)</figcaption></figure>

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

<figure><img src="image_109.png" /><figcaption>**Figure 18:** Change of the identified scale and shape parameter of the Weibull distribution for wind speeds for different levels of snow height on the ground (Dresden 1978-2007)</figcaption></figure>

Both parameters show a clear decreasing trend with increasing snow height. Hence, it seems to be justified to assume that for a storm at Dresden no snow loads have to be considered. Then, the corresponding combination factor can be assumed as 0. Many stations in the non-mountainous regions of Germany show a similar snow-wind climate, i.e. the DIN-EN-recommended com-
bination factors are far on the safe side with 0.6 for wind and 0.5 for snow for altitudes not above 1000 m.

Figure 19 presents the wind-snow climate for Rostock. The threshold value for storm is 17 m/s; 13 of 61 observed snow events with snow heights above 5 cm occur in storms. Since the amount of snow on the roof depends on the wind conditions, obviously, two scenarios have to be considered. The first scenario includes all snow events with wind speeds below storm level. For this ensemble, the observed wind speeds again show a trend to smaller values for increasing snow height. For the second scenario the interdependence of snow and wind is shown in figure 20.

<figure><img src="image_110.png" /><figcaption>Figure 19: Wind conditions for different levels of the snow height on the ground (Rostock 1954-2007)</figcaption></figure>

<figure><img src="image_111.png" /><figcaption>Figure 20: Snow heights for largest storm wind speeds (Rostock 1954-2007)</figcaption></figure>

A linear fit identifies an increasing trend, the correlation coefficient is 0.225; however, the statistical stability is poor with a probability of 0.483. Therefore, it seems to be reasonable to consider that storm wind speeds and snow heights are not correlated.

For Rostock, the two scenarios occur with different probabilities: 80% for snow and no storm, and 20% for snow and storm. Then, the exceedance probability of example given a bending moment in a structure becomes:

$$p(M > M_{\text{des}}) = 1 - (1 - p_1 \cdot p_{\text{ex},1}) \cdot (1 - p_2 \cdot p_{\text{ex},2}) \quad (42)$$

with $p_1, p_2$ – probability of scenario 1 and 2

$p_{\text{ex},1}, p_{\text{ex},2}$ – probability of $M_1$ or $M_2$ larger $M_{\text{des}}$

$M_1, M_2$ – bending moment in scenario 1 or 2

Since there is only one conditional equation, which is the target exceedance probability as specified in table 6, there is an infinite number of combinations to distribute the individual probabilities $p_{\text{ex},1}$ and $p_{\text{ex},2}$. One reasonable approach is to demand that both individual contributions are the same. Assuming a class 3 structure and 50 years projected working life specifies as target an annual exceedance probability of 1/1000. The individual targets then become:

$$ \begin{aligned} p_{\text{ex},1} &= 1 - \frac{1 - p_{\text{target}}^{1/2}}{p_1} \\ &= 1 - \frac{1 - 0.999^{1/2}}{0.2} = 0.9975 = 1 - 1/400 \end{aligned} \quad (43) $$

$$ \begin{aligned} p_{\text{ex},2} &= 1 - \frac{1 - p_{\text{target}}^{1/2}}{p_2} \\ &= 1 - \frac{1 - 0.999^{1/2}}{0.8} = 0.99937 = 1 - 1/1600 \end{aligned} $$

The final question deals with the combination factors. Strictly speaking, for the specification of the combined loads, the structural behaviour has to be considered. Then, instead of two separate load cases only a single combined load is obtained. This load can be identified based on the LRC-method [17] which is used in wind engineering to specify the effective load distribution. The following input values are required: mean values and standard deviations of the snow and wind actions and their correlation $\rho_{\text{sw}}$. Additionally, weighting factors $a_s$ and $a_w$ have to be considered which are the considered structural responses due to unit snow load and unit wind load. Finally, an appropriate target peak factor $g$ is required. The maximum structural response $R_{\text{max}}$ then is obtained from the following load combination:

$$R_{\text{max}} = a_s \cdot (\bar{s} + g \cdot \rho_{\text{Rs}} \cdot \sigma_s) + a_w \cdot (\bar{w} + g \cdot \rho_{\text{Rw}} \cdot \sigma_w) \quad (44)$$

with $\bar{x}$ – mean value and $\sigma_x$ – standard deviation

The two weighting coefficients are specified as:

$$ \rho_{\text{Rs}} = \frac{a_s \cdot \bar{s} + a_w \cdot \rho_{\text{sw}} \cdot \sigma_w}{\sigma_R} \quad (45) $$

$$ \rho_{\text{Rw}} = \frac{a_w \cdot \bar{w} + a_s \cdot \rho_{\text{sw}} \cdot \sigma_s}{\sigma_R} $$

The standard deviation of the structural response $\sigma_R$ is given with:

$$ \sigma_R = \sqrt{a_s \cdot \sigma_s^2 + a_w \cdot \sigma_w^2 + 2 \cdot \rho_{\text{sw}} \cdot a_s \cdot a_w \cdot \sigma_s \cdot \sigma_w} \quad (46) $$

The load amplitudes for the two combined loads can be expressed as a reduction factor times the individual design loads for wind and snow. The reduction factors especially depend on the ratio of the structural weighting factors. In the following example, it is assumed that
extremes of snow and wind follow the Gumbel distribution. The peak factor is set to 4.936, corresponding to an annual exceedance probability of 1/1000. The variation coefficient of the wind loads is assumed as 0.3, the corresponding value for snow is 0.45. Figure 21 shows the corresponding reduction factors. The results clearly show that the consistent codification of combined wind and snow actions is not trivial.

<figure><img src="image_112.png" /><figcaption>**Figure 21:** Reduction factors for the combined wind and snow actions</figcaption></figure>

For two uncorrelated Gumbel distributions, the accompanying load is given with the mode, which is obtained as follows:

$$x_{\text{mode}} = x_{\text{mean}} - \frac{\sqrt{6}}{\pi} \cdot x_{\text{sdev}} \cdot \gamma \quad (47)$$

Then, the two combination factors are given as:

$$\psi_{\text{sw}} = \frac{w_{\text{mode}}}{w_{\text{des}}}; \quad \psi_{\text{ws}} = \frac{s_{\text{mode}}}{s_{\text{des}}} \quad (48)$$

The two load combination then become:

$$\text{LF 1: } s_{\text{des}} \otimes \psi_{\text{sw}} \cdot w_{\text{des}} \quad (49)$$

$$\text{LF 2: } w_{\text{des}} \otimes \psi_{\text{ws}} \cdot s_{\text{des}} \quad (49)$$

For the example values, the two combination factors become $\psi_{sw} = 0.35$ and $\psi_{ws} = 0.25$. The quality of this approach can be evaluated by the ratio of the maximum structural response based on the EN-approach to the real maximum which is shown in figure 22. Additionally, the obtained exceedance probability in a projected working life of 50 years is shown. The target exceedance probability is met for cases where either the snow or the wind load is of no importance. Unacceptable solutions are obtained if the target exceedance probability of about 5% is exceeded by a factor of 2. Then, in this example, the range for $a_w/a_s$ from 0.65 to 3.50 leads to clear under-design. For the recommended DIN-EN combination factors, which are rather based on tradition than on research, the results are closer to the target. The example

points out that a revision of combination factors should not be based on a pure analysis of the actions.

<figure><img src="image_113.png" /><figcaption>**Figure 22:** Ratio of the structural response based on load combinations using combination factors to the largest response and corresponding exceedance probability of the design value</figcaption></figure>

## 7. Analysis of possible climate change

As has been shown in chapter 4, observation periods of only a few decades of snow loads on the ground do not allow identifying the underlying process. In the worst case, the extreme value statistics for an individual station will even not answer the question, if the extremes follow the type I, II or III extreme value distribution. A consistent approach considering that the true extreme snow climate remains unknown has been presented in [10] for the estimation of the design value considering a one-sided target confidence interval of 75%. The basic idea is to identify in a first step all snow climates which randomly may lead to the observed parameter triple ($\lambda_{obs}$, $sc_{obs}$, $k_{obs}$). Simulations then identify the corresponding relative frequencies for each triple ($\lambda_i$, $sc_i$, $k_i$). The final probabilities that the snow climate j leads to the observed values is given as:

$$p_j = \frac{\text{rel } f_j}{\sum_i \text{rel } f_i} \quad (50)$$

Basically, the number of events and the scale and shape parameter can be treated as independent since they do not occur in the same probability distribution. Then, the simulation effort is reduced. Figure 23 shows on the example of Dresden the huge range of possible Poisson distributions randomly leading to the observed 59 events in 28 winter periods. The obtained distribution is slightly skewed, the mean value is 60, the standard deviation is 7.753. The possible range for $\lambda$ extends from 35 to 85. In figure 24, the joint probability distribution for the possible combinations of sc and k randomly leading to the observed values is shown.

For each possible triple ($\lambda_i$, $sc_i$, $k_i$), a design snow load can be calculated. Sorting the possible design values in ascending order and adding simultaneously the corresponding $p_i$ values leads to a trace of the non-
exceedance probabilities of probable design values, assuming a specific value for the projected working life and specifying the target structural class. With the target confidence level of 75% the best estimate of the design snow load is obtained. Table 7 compares some results from the naïve approach (the observed values are used directly to estimate the design values) to the consistent probabilistic approach.

<figure><img src="image_114.png" /><figcaption>**Figure 23:** Probability density of observing 59 values in 28 winter periods for different $\lambda$-values with $\lambda = N_{\text{true}}/28$</figcaption></figure>

<figure><img src="image_115.png" /><figcaption>**Figure 24:** Joint-probability density sc-k of observing 59 values leading randomly to sc = 0.2105 kN/m², k = 0.163</figcaption></figure>

<table><caption>**Table 7:** Comparison of design values for the snow load in [kN/m²] on the ground for Dresden (1978-2007), structural class 3</caption><thead><tr><th>L</th><th>10</th><th>20</th><th>50</th><th>80</th></tr></thead><tbody><tr><td>best</td><td>0.984</td><td>1.064</td><td>1.164</td><td>1.212</td></tr><tr><td>naïve</td><td>0.863</td><td>0.915</td><td>0.975</td><td>1.000</td></tr><tr><td>ratio</td><td>1.14</td><td>1.16</td><td>1.19</td><td>1.20</td></tr></tbody></table>

Figure 25 shows the number of independent events and the corresponding maximum snow load over the time. Both traces suggest a trend. The number of events increases, while the intensity decreases over time. These findings seem to be plausible with the scenario of milder

winters. Former long periods of continuous snow cover now are interrupted by fairly warm weather; the amount of snow decreases. Fitting a linear trend to the number of events leads to 0.03 more events per year. A linear trend fitted to the snow load amplitudes identifies a decrease in the intensity of 0.002 [kN/m²] per event. It is important to note that the classical approach is based on the average value which is appropriate for the centre of the observation period. The extrapolation of the trend for the number of events to the year 2017 leads to an expected number of events of 2.886, instead of the observed average of 2.107. The expected mean value of the snow loads on the ground has to be obtained considering the increasing number of events. The expected value for the year 2017 is 0.111 kN/m². Extending the linear trend further into the future, however, soon leads to snow loads below zero. This scenario obviously is violating the basic features of the process. A physically more reasonable approach fits a linear trend to the logarithms of the observed snow loads. Then, in the limit, snow loads become zero in the future. From a structural engineer's point of view this approach seems fairly bold.

<figure><img src="image_116.png" /><figcaption>**Figure 25:** Suggested trends in the number of events per year and the amplitude of snow loads (Dresden 1978-2007)</figcaption></figure>

Once all possible combination of the parameter quintuple ($\lambda$, sc, k, $a_N$, $a_{sc}$) have been identified, joint probability distributions can be obtained. This approach has been successfully applied for the specification of the design wind speed in [18]. Figure 26 shows the joint-probability for the linear trend in the number of events per year and the true average number per year. The
random possible trend ranges from -0.05 to 0.11. The average number of observed events per year ranges between 35/28 to 95/28. The joint distribution is slightly skewed for a_N and considerably skewed for λ. For the trends in the scale the joint probabilities have three dimensions: sc, k and a_sc. The joint probability of a_sc-k for a specific sc-value is almost symmetric in regard to the value of the observed trend; in regard to k, the densities are skewed.

Figure 27 shows the joint probability for a_sc and sc accumulated over all k-values. This joint-probability is irrelevant for the identification of the best estimated of the design snow load. However, the figure illustrates that for each true scale value sc and true positive trends a_sc there is a considerable probability that the observed negative trend is obtained. For these cases, the true design snow load will be larger than the design snow load obtained with the naïve approach.

<figure><img src="image_117.png" /><figcaption>**Figure 26:** Joint probability of a_N-λ randomly leading to the observed average number of events per year and the observed linear trend (Dresden 1978-2007)</figcaption></figure>

<figure><img src="image_118.png" /><figcaption>**Figure 27:** Joint probability of a_sc-sc randomly leading to the observed scale and observed linear trend (Dresden 1978-2007)</figcaption></figure>

Similar to the approach before, for each quintuple (λ, sc, k, a_N, a_sc) a relative frequency is found based on simulations. For each combination, a design snow load is obtained. It is important to note that the exceedance probability of the design snow load is no longer constant over time, but changes with each year. The exceedance probability in the lifetime then is obtained as 1 minus the product of the non-exceedance probabilities in each individual year over the lifetime:

$$p_L(s > s_{\text{des}}) = 1 - \prod_{i=M}^{L+M} (1 - p_i (s > s_{\text{des}})) \quad (51)$$

with L – lifetime in years

p_L – exceedance probability after L years

p_i – individual exceedance probability in year i

M – starting year of exposure

The design value has to be obtained iteratively. It is important to note that the accumulated exceedance probability in the lifetime depends on the start of the exposure period.

Sorting the design snow loads for each scenario (lifetime in years, starting year of exposure and target structural class) in ascending order and summing the relative frequencies finally again leads to the trace of the non-exceedance probability of the possible design snow loads. The best estimate for each scenario is obtained with the non-exceedance probability 75%, thus meeting the confidence target.

In figure 28, the respective curves are shown for an assumed lifetime of 50 years for a class 3 structure varying the starting years from 1986 to 2026. In this example, the decreasing trend in the snow load amplitudes has more weight than the increasing trend in the number of events, i.e. the later the start of the exposure, the smaller are the design snow loads.

Table 8 summarizes the best estimates of the design snow load for exposures starting from 1986 to 2026. All five values are larger than the value obtained by applying the naïve approach, i.e. if long-term trends are considered then the snow loads increase, although the identified dominant trend is negative.

<figure><img src="image_119.png" /><figcaption>**Figure 28:** Probably design snow loads for L = 50 years and structural class considering long term trends in the snow climate (Dresden 1978-2007)</figcaption></figure>
<table><caption>**Table 8:** Best estimates for the design snow load in kN/m² for Dresden considering long-term trends in the snow climate for different starting years of exposure (L = 50 years, structural class 3)</caption><thead><tr><th>year</th><th>1986</th><th>1996</th><th>2006</th><th>2016</th><th>2026</th></tr></thead><tbody><tr><th>S<sub>des</sub></th><td>1.461</td><td>1.337</td><td>1.216</td><td>1.096</td><td>0.983</td></tr></tbody></table>

Compared to the best estimate considering a stationary snow climate, more favourable design snow loads are obtained for exposures starting in 2016 and later. It is worth mentioning that for exposures starting earlier than 2016 the design with the best estimate assuming a stationary snow climate would be on the unsafe side. In other words: The decision to exploit favourable long-term trends leads inevitably to the conclusion that existing structures are under-designed. Considering only the favourable effects and neglecting the consequences of a refined model of the snow climate for the evaluation of the existing structures is double thinking.

## 8. Summary and conclusions

The design of a structure in regard to snow actions can be separated in the following steps:

* Specification of the design target

With reference to the projected lifetime, an appropriate value of the exceedance probability of the design snow load has to be chosen. It is reasonable to distinguish structural classes, thus considering the different consequences of a failure. Since the estimation of the corresponding design value is based on confined ensembles, additionally a target value for the one-sided confidence interval has to be specified. The basic document EN 1990 recommends a confidence target of 75%.

* Specification of the design value for the snow loads on ground

The probability distribution of the underlying process can be identified based on annual extremes or on extremes above a specific threshold. An appropriate theoretical model for the annual extremes is given with the Generalized Extreme Value Distribution; the peaks-over-threshold approach uses the Generalized Pareto Distribution. The latter approach requires as additional random variable the number of events per twelve month cycle; this random parameter can be modelled by the Poisson distribution. Then, both approaches theoretically lead to identical traces of annual non-exceedance probabilities. Since usually the observation data are limited, the true distribution can hardly be identified. A fair compromise may be obtained by assuming an extreme value distribution type I for the annual non-exceedance probabilities. The uncertainties in the estimated design values depend on the size of the ensemble and the estimation

method. Generally, an adjusting factor is required to meet the target confidence. A further biasing effect may be obtained if the ensemble contains an outlier. A consistent treatment of outliers requires a high detection rate of true outliers and a small error rate for erroneously identified outliers. The strategy developed in the scope of the European research project is based on a least-square fit in the probability paper of the Gumbel distribution with plotting position of 1/(N+1) for N being the ensemble size. This method has a fairly poor detection rate. The majority of outliers will not be identified. A more consistent approach uses BLUE for the identification of outliers, leading to a considerably increased detection rate. Simply removing an identified outlier from the ensemble may lead to considerable underestimations of the design value. Specific adjusting factors are required to balance this short-coming. It is strongly recommended to revise the European national snow load maps in this regard.

* Specification of the design value of snow loads on the roof

The snow loads accumulated on a roof are influenced by the wind. It is reasonable to consider the shape factor as a random variable. Additionally, the distribution over the roof can be assumed as random. The design value for these parameters can be obtained from a convolution of the contributing probability densities. Based on the JCSS probabilistic model code, a refined proposal has been worked out. The shape factor explicitly depends on the wind speed during snow fall events. Respective statistics reveal that for many regions in Germany the wind speeds tend to zero for increasing snow loads on the ground. Then, the shape factor for flat roofs tends to unity. It is important to note that the corresponding shape factor specified in the Eurocode is 0.8 and therefore underestimates the design snow loads by 20%.

* Specification of a combined wind-snow action scenario

The Eurocode like most standards uses for the estimation of action effects under combined actions the strategy of specifying accompanying loads, which are specified based on combination factors. The specified values in Germany for snow and wind are rather based on tradition than on recent scientific analysis. For many regions in Germany, strong storms will not occur when there is snow on the ground. Then, the combination factor for accompanying snow will be zero. If additionally the wind speeds tend to small values for increasing snow load values, the corresponding wind loads become negligible small, i.e. the combination factor for accompanying wind is zero as well. In other regions, storm and snow may occur simultaneously. Then, two scenarios have to be covered: snow without storm and snow together with storm. The joint probability density for the latter scenario can be modelled in a first approach as uncorrelated extreme value distributions type I. Then, the accompanying load amplitude is given with the mode of the respective probability distribution. The
final evaluation of the combined wind and action effects has to be based on structural responses. Basically, the method of accompanying loads tends to underestimate the structural responses. For any refined approach it is therefore strongly recommended to consider the structural behaviour.

* The treatment of a possible climate change or of natural long-term trends in the wind-snow climate will be one of the challenges in future codification work. A consistent approach accepts that the true parameters of the probability distributions of the underlying physical processes cannot be identified since the ensembles of observed data are much too small. The same conclusion is true for the basic shape of trends. A reasonable approach therefore describes the trends by a simple linear law. For larger negative trends, the amplitudes in the future may become smaller than zero, thus leading to physical inconsistencies. For these cases, the snow load simply can be set to zero, and the corresponding non-exceedance probability becomes unity. A refined shape of negative trends is based on a linear model for the logarithms of the respective variable, thus avoiding negative values. Based on simulations, a range of possible true climates and trends can be identified which randomly lead to the observed values. In a next simulation step, the corresponding joint-probability densities for the parameter quintuple (a<sub>N</sub>, λ, a<sub>sc</sub>, sc, k) can be obtained. Each possible combination of the parameters leads to a specific design value. Sorting these design values in ascending order together with the corresponding probability and integrating the joint-probability densities finally leads to a trace of random design values. The best estimate of the design value is obtained considering the target confidence. Even with the bold approach of accepting zero snow loads in the future, for the chosen example of Dresden the best-estimates of the design snow load are larger than those values obtained with the naïve approach of assuming that the observed values are the true values and long-term trends do not exist. It is important to note that negative trends may be observed randomly even if the true trend is positive. While for dominant positive trends the design values increase in the future, a dominant negative trend leads to a decrease of future design values. On the other hand, existing structures may turn out to be under-designed.

## Acknowledgement

The author thanks the organizers of ICSE 2016 for the kind invitation to deliver this keynote paper. Part of the presented work has been sponsored by the German Research Foundation under different contract numbers for KA675 and by the Federal Ministry of Research and Education as part of the joint-research project RegioExAKT. The German weather service DWD (Deutscher Wetterdienst) has provided the meteorological data. These supports are gratefully acknowledged.

## References

1. ISO 22111 (2007): Bases of design of structures – General requirements

2. ISO 4354 (2009): Wind actions on structures

3. W.H. Press, S.A. Teukolsky (1992): Numerical Recipes in Fortran 77: The Art of Scientific Computing, Cambridge University Press,

4. M. Matsumoto, T. Nishimura (1998): Mersenne twister. A 623-dimensionally equidistributed uniform pseudo-random number generator, ACM Transactions on Modeling and Computer Simulation, 3-30

5. M. Matsumoto (2010) : Mersenne Twister homepage, A very fast random number generator of period $2^{19937}-1$, www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT

6. G. Marsaglia (2003): Random number generators, Journal of Applied Statistical Methods, 2, 2-13

7. BAST – Bundesanstalt für Straßenwesen (2009) Effects of heavy goods vehicle traffic on the bridges of the federal roads (in German)

8. E.J. Gumbel (1958): Statistics of Extremes, Columbia University Press, New York · London

9. I.I. Gringorten (1963): A plotting rule for extreme probability paper, Journal of Geophysical Research, 68, 813-814

10. M. Kasperski (2012): A consistent approach for estimating the design value of the snow load on the ground from confined ensembles, in Proc. of 7th Int. Conf. on Snow Engineering, Fukui, Japan

11. B. Czwikla, M. Kasperski (2016): Specification of the design value of the ground snow load considering measurements of the snow height – part 1: single stations, ICSE 2016, Nantes

12. L. Sanpaolesi et al. (1998): Phase 1 Final Report to the European Commission, Scientific Support Activity in the Field of Structural Stability of Civil Engineering Works:Snow Loads. Dept of Structural Engineering, University of Pisa

13. N. Balakrishnan, P.S. Chan (1992) Order statistics from extreme value distribution, II: Best linear unbiased estimates and some other uses, Commun. Statist.-Simula., 21(4), 1219-1246

14. A. Childs, N. Balakrishnan (2002) Series approximations for moments of order statistics using MAPLE, Computational Statistics & Data Analysis 38, 331-347

15. JCSS, Probabilistic Model Code Part 2: Load Models, 2.12 Snow Load, 99-CON-DYN/M0097 February 2002 1 JCSS-VROU-07-02-96

16. B. Czwikla, M. Kasperski (2012) Design snow loads for gable roofs, Proc. 7 ICCSM - International Congress of Croatian Society of Mechanics, Zadar, Croatia, May 22 - 25, 2012, Zadar

17. M. Kasperski (1992) Extreme wind load distributions for linear and nonlinear design. Engineering Structures 14, 27-34

18. M. Kasperski (2010) A new approach for the specification of the design wind speed considering long-term trends, Global Environmental Research, Vol. 13 No. 2, 161-168
# Climate Change: impact on snow loads on structures

Pietro Croce¹, Paolo Formichi¹,ᵃ, Filippo Landi¹ and Francesca Marsili¹

¹University of Pisa, Department of Civil and Industrial Engineering-Structural Division, Largo Lucio Lazzarino - 56126 Pisa, Italy

**Abstract.** A general procedure to evaluate future trends in snow loads on structures is illustrated aiming to study influences of climate change at European scale, to assess its impact on the design of new structures as well as on the reliability levels of existing ones, also in view of the next revision of the Eurocodes. Analysing high quality registered meteorological data of daily temperatures, rain and snow precipitations in nine Italian weather stations, conditional probability functions of occurrence of snow precipitation, accumulation and melting have been preliminarily determined as functions of daily air temperatures. By means of Monte Carlo simulations and based upon daily output of climate models (daily max. and min. temperatures and water precipitation) yearly maxima of snow loads for various time intervals of 40 years in the period 1980-2100 have been simulated, deriving, via the extreme value theory, the characteristic ground snow loads at the sites. Then, the proposed procedure has been implemented in a more general methodology for snow map updating, in such a way that the influence of gridded data of precipitation, predicted by global climate models, on extreme values of snow loads is duly assessed. Preliminary results demonstrate that the outlined procedure is very promising and allows to estimate the evolution of characteristic ground snow loads and to define updated ground snow load maps for different climate models and scenarios.

## 1 Introduction

As structural design is often governed by climatic actions, alterations of climatic actions caused by climate change could significantly impact design of new structures as well as the reliability of existing ones. For this reason some research groups, based in different areas of the world, are investigating the problem, especially aiming to study the possibility to transfer the outputs of climate models which refer to geographical cells of large dimensions on a suitably smaller scale, allowing to derive characteristic values of climatic actions from the above mentioned outputs. In particular, effects of climate change on snow loads have been recently investigated in Germany [1], Norway [2] and Canada [3], focusing on their influence on the reliability of built environment.

A general procedure for the estimation of climate change influence on ground snow loads is presented; such a procedure, combining the outputs of global climate models with measured data, allows to refine the snow load predictions arriving to a geographical resolution adequate to take into account also micro-climate effects, depending on the local orography.

The procedure, which is general enough to permit sound estimation of future trends in snow loading on structures, aims more specifically to address the impact of climate change on the assessment of new and existing structures, also in view of the evolution of the second generation of Eurocodes, and Eurocode EN 1991 - Part 1-3 - Snow loads [4] in particular, as requested by Mandate M/515 of European Commission [5] to CEN [6].

## 2 Snow load on structures: State of Art

### 2.1 Definition of ground snow load

The definition of snow loads on structures in the current version of the structural Eurocodes is largely based upon the results of the European Snow Load Research Project

(ESLRP) [6] and, in particular, on the European Ground Snow Load Map elaborated within that research. This map was the first snow load map derived at European scale with a strong scientific basis [7] and its elaboration started from the analysis of the basic snow data collected across 18 European countries, which at that time were members of CEN. By applying extreme value statistics to these observations, the characteristic ground snow loads were derived corresponding to a given probability of exceedance (the 2% upper percentile on yearly base, corresponding to a return period of about 50 years), taking into account both the influence of exceptional snowfalls or no snowy winters. It must be also stressed that in [7], exceptional snow loads were considered for the first time at European scale, on sound scientific bases, giving also appropriate guidance for their identification.

### 2.2 European snow load maps

Beside the characteristic ground snow loads at weather stations, the research allowed to identify ten major European climatic regions, and the corresponding snow maps, elaborated through GIS software: generally, in each climatic region the ground snow load depends on the altitude according to suitably derived law, except in Norway and Iceland, where snow load is practically independent of altitude. Based on the above snow map, each National Standard Body of CEN member states, produced their own national map of ground snow load, published in the National Annex to EN 1991-1-3.

## 3 Problem statement

### 3.1 Climate change impacts

The evidence of climate change is unequivocal [8]: the average global temperature, currently around 0,8°C above pre-industrial level, continues to rise, even more

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

<footnote>ᵃ Corresponding author: p.formichi@ing.unipi.it</footnote>
evidently in Europe [9]. Consequences of climate change, like alteration of natural processes, modifications of precipitation patterns, melting of glaciers, rise of sea levels and so on, are increasingly evident. For this reason, whatever the warming scenarios and the level of success of mitigation policies, in the coming decades the impact of climate change needs to be considered, taking into account possible delay in achieving targets in reduction of greenhouse gas emissions as well as adaptation measures to deal with economic, environmental and social consequences of climate change [10].

A global assessment of climate change science is presented by the recent IPCC (International Panel on Climate change) Report [11]. The assessment considers new evidence of past, present and projected future climate change foresees by several independent scientific studies, based on observations of the climate system, theoretical studies of climate processes and simulations using climate models. Predictions of climate change are based on a set of four Representative Concentration Pathways (RCPs) scenarios, where radiative forcing target level, estimated for the year 2100, are fixed in comparison with pre-industrial values [12]. Globally, these scenarios generally lead to further warming and changes in the water cycle, but trends can be different or even opposite on local scale.

The EU response to climate change is an adaption strategy to enhance the capacity to withstand it and the readiness to respond to its impacts, particularly in most vulnerable key sectors like infrastructures and buildings, characterised by a long life span and high costs [10]. In this respect, a central role is played by technical standards and by their evolution during the lifetime cycle of the infrastructures and buildings [13]; implications of climate change for new and existing structures is then a key aspect in the development of second generation of Eurocodes [6]. In this context, the present study aims to determine the effects of Climate Change on snow loads on structures, implementing suitable numerical procedures to derive characteristic snow loads from available global and regional climate models.

## 3.2 Snow load in a changing climate

One of the trite remarks about global warming is that, as obvious consequence of it, a reduction of snow is generally expected, but it is not said. In fact, although in consequence the global warming frequency of snow events should reduce, the intensity of extreme snow events may increase, since the capacity of the atmosphere to hold moisture (and to form snow particles in case of sudden temperature's drop) increases with temperature. This may lead to the increase of both snow density and occurrence of extreme snowfalls in regions where temperatures still may happen to be below freezing level during precipitation events [1].

The above considerations are confirmed by some spectacular and catastrophic collapses of lightweight roof structures caused by snow, which occurred in the last years all around the Europe [14] and [15] like the roof collapse of the Katowice Exhibition Hall (2006), which

caused 65 victims and over 170 injured people (Figure 1). A more wide illustration of some relevant failures can be found in [16], while the reliability level of existing structures designed with the current versions of Eurocodes, especially for those structures more sensitive to an increase of snow loads, like lightweight structures of wide-span roofs, is discussed in [16] and [17].

<figure><img src="image_121.png" /><figcaption>Figure 1. Roof collapse of Katowice Exhibition Hall</figcaption></figure>

To assess the evolution of the ground snow load and its impact on new and existing structures, future trends, in terms of intensity and frequency of precipitations in cold areas, should be compared with those adopted to develop current versions of snow load maps for structural design.

A suitable numerical methodology is proposed in order to estimate characteristic ground snow loads from available results of global and regional climate models (GCMs and RCMs) in the mid to long term time, according different scenarios RCPs.

As just said, the interpretation of the outputs of climate models is complicated by their spatial resolution, typically too coarse to study specific sites, which does not allow a direct evaluation of snow load trends. For this reason, the proposed procedure has been calibrated and validated against high quality measured data available for time series of 50 or more years and particular attention has been devoted to identify suitable techniques to check the homogeneity of the available data populations.

Finally, it must be stressed that intrinsic uncertainty affecting climate projections in different climate models reflects on their outcomes, so that conclusions about increase/decrease of actions need careful interpretations.

# 4 Procedure for the estimation of future trends of ground snow loads

The proposed procedure consists of four steps:

- analysis of observed data series in order to derive conditional probability functions, linking snowfall and snow melting conditions with precipitation data and air temperature;

- development of a predictive model to evaluate snow loads from available meteorological data;

- calibration and validation of the model predictions against observed data series;

- implementation of the model on projected data series from global climate change models.
## 4.1 Analysis of observed data series

Relevant meteorological data recorded in the past, daily temperatures and precipitation (water equivalent and snow cover depth), have been collected and analysed with the purpose to seek the conditions of maximum and minimum daily temperature, $T_{max}$ e $T_{min}$, at which snow cover depth increases, in case of precipitation of height $h_r$, or decreases, melting or increasing densities, and rainfalls are likely to occur.

In particular, seven different situations have been identified comparing day n and day n-1:

1. total melting of snow cover present at day n-1;

2. partial melting of the snow cover present at day n-1;

3. constant snow cover depth without precipitation;

4. rainfall in absence of snow cover at day n-1;

5. rainfall and total melting of the snow cover present at day n-1;

6. precipitation on snow cover with snow depth decreasing;

7. snowfall with increasing of snow cover;

and, following the flowchart shown in figure 2, the daily measured data can be allocated according to them.

<figure><img src="image_122.png" /><figcaption>Figure 2. Flowchart for allocation of daily measured data</figcaption></figure>

For each relevant situation, the frequency histograms have been then derived according to the conditions

$$ Z(\bar{T}_{min}, \bar{T}_{max}) = \text{number of cases of each situation for which} \\ \left( \bar{T}_{max} - \frac{\Delta T}{2} \right) < T_{max,n} < \left( \bar{T}_{max} + \frac{\Delta T}{2} \right) \\ \left( \bar{T}_{min} - \frac{\Delta T}{2} \right) < T_{min,n} < \left( \bar{T}_{min} + \frac{\Delta T}{2} \right) \\ -20^{\circ}\text{C} < \bar{T}_{min} < 40^{\circ}\text{C}; -20^{\circ}\text{C} < \bar{T}_{max} < 40^{\circ}\text{C}; \Delta T = 1^{\circ}\text{C} \quad (1) $$

The histogram of frequencies for each of the seven situations are then converted into continuous surfaces $f_j(T_{max}, T_{min})$ by means of up to three two-dimensional Gaussian functions

$$ f_j = \sum_i a_i e^{-\frac{1}{2} \left( \frac{\left[ \left( T_{max} \cos \beta_i + T_{min} \sin \beta_i \right) - \left( \mu_{max,i} \cos \beta_i + \mu_{min,i} \sin \beta_i \right) \right]^2}{\sigma_{max,i}} + \frac{\left[ \left( T_{min} \cos \beta_i - T_{max} \sin \beta_i \right) - \left( \mu_{min,i} \cos \beta_i - \mu_{max,i} \sin \beta_i \right) \right]^2}{\sigma_{min,i}} \right)} \quad (2) $$

being i the number of the Gaussian functions, which is assumed equal to the number of local maxima of the histogram.

In eq. (2), each Gaussian function is defined by 6 parameters:

- $a_i$ amplitude

- $\beta_i$ angle of rotation with respect to x-axis ($T_{max}$)

- $\mu_{max,i}$ median value with respect to x-axis ($T_{max}$)

- $\mu_{min,i}$ median value with respect to y-axis ($T_{min}$)

- $\sigma_{max,i}$ standard deviation with respect to x-axis ($T_{max}$)

- $\sigma_{min,i}$ standard deviation with respect to y-axis ($T_{min}$).

The parameters $\beta_i$, $\mu_{max,i}$, $\mu_{min,i}$ are derived from the peak values of the histograms, while the other parameters, $\sigma_{max,i}$, $\sigma_{min,i}$, $a_i$ are estimated by means of least squares method and imposing the condition of unitary volume underlying the surface ($\sum_i a_i = 1$).

Once derived the probability distribution function of the seven conditions $f_j$ for each relevant condition, $j=1,...,7$, it is possible to define conditional probability functions of snowfall and snow melting for given values of the daily temperatures, $\bar{T}_{max}$ and $\bar{T}_{min}$. In particular, it results that

- the conditional probability function of snowfall in presence of precipitation is

$$ P(\text{snowfall}(\bar{T}_{\text{max}} \bar{T}_{\text{min}}) | \text{prec}) = \frac{n_7 \cdot f_7}{n_7 \cdot f_7 + n_4 \cdot f_4} \quad (3) $$

where $n_7$ is the number of cases of effective snowfall (condition 7) and $f_7$ is the consistent probability function previously defined; $n_4-5$ is the number of cases of rainfall (conditions 4 plus 5) and $f_4-5$ is the consistent probability function (an example of such a function for the Italian weather station of Spigno Monferrato is reported in Figure 3);

<figure><img src="image_123.png" /><figcaption>Figure 3. Conditional probability function of snowfall in presence of precipitation – (Spigno Monferrato – I)</figcaption></figure>

- the conditional probability function of melting of snow cover in presence of snow cover is

$$ P(\text{snowmelt}(\bar{T}_{\text{max}} \bar{T}_{\text{min}}) | \text{snowcover}) = \frac{n_{1-2} \cdot f_{1-2}}{n_{1-2} \cdot f_{1-2} + n_3 \cdot f_3} \quad (4) $$

where $n_{1-2}$ is the number of cases of snow melting (conditions 1 plus 2) and $f_{1-2}$ is the consistent probability function; $n_3$ is the number of cases of constant snow cover (condition 3) and $f_3$ is the consistent probability function previously defined and, obviously, $\bar{T}_{max} > 0^{\circ}\text{C}$ (an example, for the Italian weather stations is reported in Figure 4)
<figure><img src="image_124.png" /><figcaption>Figure 4. Conditional probability function of snow melting in presence of snow cover – Italy</figcaption></figure>

## 4.2 Predictive model

A predictive model to evaluate ground snow loads has been then developed based on the previously determined conditional pdfs. In Figure 5, it is presented the flowchart of the algorithm which has been implemented through a suitable Monte Carlo simulation.

The input data of the algorithm are three relevant meteorological daily data: the maximum and minimum air temperatures, $T_{max,n}$, $T_{min,n}$, and the precipitation in mm of water, $h_{r,n}$, at the n-th day.

The probability of snowfall with increasing snow cover depth is estimated by checking the following conditions:

$$T_{min,n} < 5^{\circ}\text{C} \land h_{r,n} \neq 0 \land P(\text{snowfall}(\bar{T}_{max,n}, \bar{T}_{min,n}) | \text{prec}) > R, (5)$$

where $0 \le R \le 1$ is derived from the randomly generated number, modified with the hypercube latin sampling technique.

When conditions (5) are satisfied, the increase of the ground snow load $\Delta q_n$ at the n-th day is estimated in terms of water equivalent measured by the rain gauge:

$$\Delta q_n = 0,01 \cdot h_{r,n} \text{ kN/m}^2 \text{ with } h_{r,n} \text{ in [mm]} . \quad (6)$$

At the (n+1)-th day, several alternative and mutually exclusive events can happen, with an associated probability of occurrence:

- if snowfall conditions (5) are satisfied again, the ground snow load $q_{n+1}$ is calculated as

$$q_{n+1} = q + 0,01 \cdot h_{r,n+1} \text{ kN/m}^2 \quad (7)$$

- if snow melting conditions are satisfied, i.e.

$$P(\text{snowmelt}(\bar{T}_{max,n+1}, \bar{T}_{min,n+1}) | \text{snowcover}) > R \quad (8)$$

snow melts, partially or totally, being R defined before: the melting is assumed to be proportional to the value of $\bar{T}_{max,n+1}$, so that the updated ground snow load becomes

$$q_{n+1} = q_n \cdot \left(1 - \left(\frac{T_{max,n+1}}{20}\right)\right) \text{ where } 0^{\circ}\text{C} \le \bar{T}_{max,n+1} \le 20^{\circ}\text{C} \quad (9)$$

but, if $\bar{T}_{max,n+1} > 20^{\circ}\text{C}$ or $q_{n+1} < 0,025 \text{ kN/m}^2$ total melting is considered.

- if rainfall precipitation occurs in case of snow cover, the new ground snow load is given by

$$q_{n+1} = q_n + h_{r,n+1} \cdot 0,01 \text{ kN/m}^2 \le 1,25 \cdot q_n \quad (10)$$

but, if $q_n < 0,01 \cdot h_{r,n+1}$ total melting is considered.

<figure><img src="image_125.png" /><figcaption>Figure 5. Flow chart of the algorithm for the estimation of yearly maximum ground snow load</figcaption></figure>
The Monte Carlo simulation has been repeated for all days of the year, in order to estimate the yearly maximum ground snow load, $q_{max}$. For each year, the process has been iterated 10'000 times: the median value of the distribution of each series has been taken as the best estimate of the maximum ground snow load for the examined year. Concerning the number of iterations, calibration exercises proved that $10^4$, $10^5$ or $10^6$ iterations give practically coincident results.

From the set of the N yearly simulated maxima, the characteristic value $q_k$, associated to the probability of 2% to be exceeded in one year, has been finally derived, via extreme values analysis, as sketched in Figure 6.

<figure><img src="image_126.png" /><figcaption>Figure 6. Characteristic value of ground snow load</figcaption></figure>

## 4.3 Calibration of the model

The methodology illustrated above has been initially tested against the observed data series of nine Italian weather stations (Figure 7), where high quality daily data are available for long time series (more than 50 years).

These weather stations have been selected between the 125 Italian stations, whose data were collected in the Database of the University of Pisa during the European Snow Load Research Project (ESLRP) [7].

The results compared with those obtained by ESLRP (Table 1) show that the outlined procedure is very promising and allows not only to derive characteristic ground snow loads, but also to assess the effects on solid precipitation measurements of snow transport phenomena caused by wind (Figure 8) [18] through an appropriate parameter, $K_{corr}$.

The parameter $K_{corr}$, which is characteristic of each climatic area, takes into account the systematic errors of snow precipitation measurements at the rain gauge [19] [20] which are mainly due to wind field deformation above the gauge orifice (Figure 9), and depend on the wind speed, the intensity of precipitation and the type of the gauge [19].

## 4.3 Analysis of climate projections

Once performed the calibration of the procedure against registered data series, which allowed to estimate the

correction factor $K_{corr}$, the analysis has been extended to projected data series provided by different climate models (GCMs and RCMs) and different scenarios (RCPs) up to 2100. In particular, they have been used available high resolution data EUR11 (grid resolution of 12.5 km) developed within the EUROCORDEX initiative [21] [22] for the control period 1981-2005 (Historical Experiment) where “run” is forced by observed atmospheric composition changes[23] and for the future period 2006-2100 (RCPs Experiment) where “run” is forced by scenarios RCPs[23] illustrated in Table 2

The analysis has been carried out for subsequent time windows, each spanning over forty years; the obtained estimates of characteristic ground snow loads ($q_k$) are illustrated in Figure 10 for some of the nine tested Italian weather stations.

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

<table><caption>Figure 7. Italian weather station used in the analysis</caption><thead><tr><th>Country</th><th>Station</th><th>LON [°]</th><th>LAT [°]</th><th>Altitude [m]</th><th>q<sub>ESLRP</sub> [kN/m²]</th></tr></thead><tbody><tr><td>IT</td><td>Ascoli Piceno</td><td>13,56</td><td>42,90</td><td>136</td><td>1,12</td></tr><tr><td>IT</td><td>Bologna</td><td>11,36</td><td>44,50</td><td>51</td><td>1,67</td></tr><tr><td>IT</td><td>Castelnuovo Garfagnana</td><td>10,43</td><td>44,12</td><td>276</td><td>1,14</td></tr><tr><td>IT</td><td>Lodi</td><td>9,51</td><td>45,32</td><td>80</td><td>1,19</td></tr><tr><td>IT</td><td>Massafra</td><td>17,13</td><td>40,58</td><td>116</td><td>0,54</td></tr><tr><td>IT</td><td>Pesaro</td><td>12,91</td><td>43,92</td><td>11</td><td>1,07</td></tr><tr><td>IT</td><td>San Giuseppe Jato</td><td>13,19</td><td>37,97</td><td>450</td><td>0,94</td></tr><tr><td>IT</td><td>Spigno Monferrato</td><td>8,36</td><td>44,55</td><td>476</td><td>2,40</td></tr><tr><td>IT</td><td>San Severo</td><td>15,39</td><td>41,68</td><td>87</td><td>0,74</td></tr></tbody></table>

<table><caption>Table 1. Results of the algorithm for ESLRP data.</caption><thead><tr><th rowspan="2">Country</th><th rowspan="2">Station</th><th colspan="2">Characteristic Ground<br>Snow Load [kN/m²]</th><th rowspan="2">K<sub>corr</sub> = q<sub>k</sub>/q<sub>ESLRP</sub></th></tr><tr><th>Tested<br>Procedure<br>q<sub>k</sub></th><th>ESLRP [7]<br>q<sub>ESLRP</sub></th></tr></thead><tbody><tr><td>IT</td><td>Ascoli Piceno</td><td>0,68</td><td>1,04</td><td><strong>1,53</strong></td></tr><tr><td>IT</td><td>Bologna</td><td>0,99</td><td>1,66</td><td><strong>1,68</strong></td></tr><tr><td>IT</td><td>Castelnuovo<br>Garfagnana</td><td>1,14</td><td>1,20</td><td><strong>1,05</strong></td></tr><tr><td>IT</td><td>Lodi</td><td>0,76</td><td>1,20</td><td><strong>1,58</strong></td></tr><tr><td>IT</td><td>Massafra</td><td>0,16</td><td>0,54</td><td><strong>3,38</strong></td></tr><tr><td>IT</td><td>Pesaro</td><td>0,67</td><td>1,07</td><td><strong>1,60</strong></td></tr><tr><td>IT</td><td>San Giuseppe Jato</td><td>0,50</td><td>1,21</td><td><strong>2,42</strong></td></tr><tr><td>IT</td><td>Spigno Monferrato</td><td>0,88</td><td>2,55</td><td><strong>2,90</strong></td></tr><tr><td>IT</td><td>San Severo</td><td>0,46</td><td>0,80</td><td><strong>1,74</strong></td></tr></tbody></table>
<figure><img src="image_128.png" /><figcaption>Figure 8. Snow transport phenomena [18]</figcaption></figure>

<figure><img src="image_129.png" /><figcaption>Figure 9. Wind field deformation at the rain gauge [14]</figcaption></figure>

<table><caption>Table 2. Climate models and scenarios</caption><thead><tr><th>RCM</th><th>Resolution</th><th>Driving GCM</th><th>Driving Experiment</th><th>Period</th></tr></thead><tbody><tr><td rowspan="3">DMI-HIRHAM5</td><td rowspan="3">EUR11<br>(0,11"=12.5km)</td><td rowspan="3">EC-EARTH</td><td>Historical</td><td>1981-2005</td></tr><tr><td>RCP4.5</td><td>2006-2100</td></tr><tr><td>RCP8.5</td><td>2006-2100</td></tr><tr><td rowspan="3">CNRM-ALADIN5.3</td><td rowspan="3">EUR11<br>(0,11"=12.5km)</td><td rowspan="3">CNRM-CM5</td><td>Historical</td><td>1981-2005</td></tr><tr><td>RCP4.5</td><td>2006-2100</td></tr><tr><td>RCP8.5</td><td>2006-2100</td></tr><tr><td rowspan="3">KNMI-RACMO22E</td><td rowspan="3">EUR11<br>(0,11"=12.5km)</td><td rowspan="3">HadGEM2-ES</td><td>Historical</td><td>1981-2005</td></tr><tr><td>RCP4.5</td><td>2006-2090</td></tr><tr><td>RCP8.5</td><td>2006-2100</td></tr></tbody></table>

Finally, to assess the variability of climate models and scenarios, for the nine investigated stations the characteristic values of ground snow loads associated to different scenarios ($q_{KRCP4.5}$ and $q_{KRCP8.5}$) have been compared, adopting different climate models.

Considering 40 years intervals shifted by 10 years, the ratio between $q_{KRCP8.5}$ and $q_{KRCP4.5}$ has been determined for the period 1981-2100, obtaining the trends shown in Figure 11.

The results show a general decrease of characteristic snow loads for the investigated stations, but the expected decreasing trend of the ratio $q_{KRCP8.5}/q_{KRCP4.5}$, namely a more marked decrease of snow loads linked to increasing emissions, is not confirmed for all of them.

<figure><img src="image_130.png" /><figcaption>Figure 10. Future trends of $q_k$ [kN/m²] for different climate models and scenarios</figcaption></figure>
<figure><img src="image_131.png" /><figcaption>Figure 11. Trend of the ratio between qkRCP8.5 and qkRCP4.5 for different climate models and different stations</figcaption></figure>

## 5 Preliminary results on snow map comparison

The proposed procedure seems very appropriate for creation of snow maps taking into account the effects of climate change, since it allows to estimate characteristic ground snow loads on the basis of daily data for the minimum and maximum temperature and precipitation, which are typically available as outputs of climate change projections for all possible scenarios.

Then, extending the previous analysis to a larger number of weather stations it is possible to obtain enough local information to update snow maps. The updated versions of the maps, critically compared with the current versions, will provide so a comprehensive overview of the impact of climate change, according to various RCPs scenarios [4].

However, before arriving to obtain consistent snow maps, further calibration is needed.

Indeed, climate projections are provided, in the higher resolution models, for cell of about 12,5×12,5 km (EUR11 - 0,11° on a rotated grid) and are generally agreed to represent area averaged data rather than point

process data, especially for precipitation data [24]. Therefore, the analysis of gridded precipitation data can lead to an underestimation of extremes [24] [25].

### 5.1 Comparison of observed gridded data and observed point data

In order to assess a calibration method for the analysis of gridded datasets, able to correct the underestimation of extremes, the above mentioned issue has been especially studied for Italy, comparing characteristic ground snow loads obtained using grid cell data (E-OBS dataset) with those obtained by the European Snow Load Research Project (ESLRP) for the Mediterranean region [7]. In Figure 12, the location of the Italian weather stations considered in ESLRP is reported on the grid cell map of E-OBS.

The results confirm that characteristic values computed by the analysis of gridded data are typically lower than values obtained by point-data and that the difference

$$ \Delta q_k = q_{kESLRP} - q_{kEOBS} \quad (11) $$

generally increases with the site altitude (Figure 13).
Starting from these results, a regression analysis has been carried out and a relationship between the difference of characteristic values $\Delta q_k$ and altitude above the sea level A has been estimated as

$$ \Delta q_k = a \cdot \left[ 1 + \left( \frac{\Delta a}{b} \right)^2 \right] \frac{kN}{m^2} \quad \text{with A in [m]} \quad (11) $$

where the parameters a and b have been calculated for each one of the four climatic zones defined by EN1991-1-3 for Mediterranean region in Italy obtaining the curves illustrated in Figure 14, 15 and 16. As in zones 3 and 4 the amount of available data is limited, the zones themselves have been merged.

<figure><img src="image_132.png" /><figcaption>Figure 12. EOBS - Grid Cells and ESLRP weather stations for Italy</figcaption></figure>

<figure><img src="image_133.png" /><figcaption>Figure 13. $\Delta q_k$ - Altitude (Italy - Mediterranean Region)</figcaption></figure>

Once these relationships are available, it is possible to draw, according to the results obtained by the analysis of gridded data, the snow load map on the ground, corrected in such a way that the data are referred to the sea level. This correction is very important because it allows to evaluate also results obtained for characteristic ground

snow loads by the investigation of projected data with respect to the current code's provisions. (Figure 17).

Then, the next step will be the analysis of projected data, provided by different climate models, in order to derive maps of future trends of characteristic ground snow loads.

## 5.2 Actual trend of variation of ground snow load

Before the study of projected data, EOBS gridded data for Italian region have been analysed for subsequent time periods of thirty years in order to estimate both the actual trend of variation of characteristic ground snow loads and its future evolution.

The results are illustrated in Figure 18 and show different trends for Italian regions.

<figure><img src="image_134.png" /><figcaption>Figure 14. $\Delta q_k$ - Altitude (Italy - Mediterranean region - zone 1)</figcaption></figure>

<figure><img src="image_135.png" /><figcaption>Figure 15. $\Delta q_k$ - Altitude (Italy Mediterranean region-zone 2)</figcaption></figure>

<figure><img src="image_136.png" /><figcaption>Figure 16. $\Delta q_k$ - Altitude (Italy Mediterranean region-zones 3-4)</figcaption></figure>
<figure><img src="image_137.png" /><figcaption>**Figure 17.** Comparison of Snow Load Map of EN1991-1-3 and Snow Load Map obtained using EOBS data corrected with altitude relationship</figcaption></figure>

<figure><img src="image_138.png" /><figcaption>**Figure 18.** Snow load maps obtained using different periods of EOBS data</figcaption></figure>

# 6 Conclusions

Implication of climate change on snow loads could have very high impact on design of structures as well as in assessment of existing ones. The problem is very relevant as it is demonstrated by a high number of structural collapses of lightweight structures due to snow.

For this reason, the setup of a suitable procedure to predict the snow load on ground according to various climate change scenarios could represents a very strong improvement in definition and adaption of snow load maps. This issue is particularly relevant in Europe, also in view of the development of future generation of Eurocodes.

Aiming to tackle this relevant issue, a general procedure to evaluate future trends of ground snow loads
for structural design, taking into account the influence of
climate change is proposed, based on Monte Carlo
method.

The procedure, starting from the analysis of historical
meteorological observations, allows to estimate ground
snow loads and their characteristic values, on the basis of
daily outputs of climate models, which are typically
available (daily temperature extremes Tmax, Tmin and
height of precipitation hr).

The preliminary implementation of the outlined
procedure is very promising and shows that there is a
concrete possibility to arrive to define a general
methodology for the estimate of characteristic snow loads
on ground, according to prediction derived from
climatological models, also on a local scale. Since the
presented methodology seems to be able to appreciate the
evolution of the snow load during the time, it should also
allow proper refinement and updating of the snow load
maps, according the real trends.

Extending the analysis to predicted data series
covering a large enough number of weather stations
across Europe, it will be possible to compare the future
trends of snow loads, for each available emission
scenario, with the current European snow load maps,
assessing their evolution and their impact on structural
design, in line with the requests by the European
Commission to CEN in Mandate M515 within the
European adaptation to climate change policy, in view of
the development of the second generation of Eurocodes.

References

1. U. Strasser, Natural Hazard and Earth System Science, **8**, p.1-8 (2008)

2. H.T.T. Tayet, H.O. Hygen and T. Kvande., EMS 2013-230, **10** (2013)

3. CAN/CSA - S502-14 - *Managing changing snow load risk for building in Canada's North* (2014)

4. UNI EN 1991-1-3:2004 – Eurocode 1: Actions on structures – Part 1-3: General Actions – Snow Loads

5. M/515 EN - *Mandate for amending existing Eurocodes and extending the scope of Structural Eurocodes* (2012)

6. CEN/TC250 - *Response to Mandate M/515 - Towards a second generation of Eurocodes* (2013)

7. DGII-D3 - Scientific support activity in the field of structural stability of civil engineering works. Snow loads - FINAL REPORT (1996)

8. <http://climate.nasa.gov/evidence/>

9. European Environment Agency, *Climate Change, impacts and vulnerability in Europe* - EEA report n°12 (2012)

10. European Commission, *An EU strategy on adaptation to climate change* -COM 213 Final report (2013)

11. WG1AR5-IPCC, *Climate change 2013 – The physical Science Basis* -

12. D. P. Van Vuuren, J. Edmonds, M. Kainuma et al., *Climatic Change*, **109**, 5-31 (2011)

13. European Commission, *Adapting infrastructure to climate change* -SWD 137 Final report (2013)

14. J. Geis, K. Strobel and A. Liel, Journal of Performance of Construction Facilites, **26(4)**, 377-388(2012)

15. E. Früwald, E. Serrano, T. Toratti, et al., *Report TVBK-3053*, Lund University, Lund (2007)

16. M. Holicky, Engineering Sciences, **58**, 51-57 (2007)

17. M. Holicky and M. Sykora, Forensic Engineering, 2009:444-453

18. M. Mellor, Cold Regions Science and Engineering Part III, Section A3c (1965)

19. B.E. Goodison, P.Y.T. Louie and D. Yang, WMO/TD N872 (1998)

20. A. Lendvai, R. Ranzi, G. Peretti et al., AINEVA-Neve e Valanghe, **84**, 12-21 (2014)

21. D. Jacob, J. Petersen, B. Eggert et al., Reg Environ Change, **14**, 563-578 (2014)

22. S. Kotlarsky, K. Keuler, O.B. Christensen et al., Geoscientific Model Development, **7**, 1297-1333.

23. K.E. Taylor, R.J. Stouffer and G.A. Meehl, Bulletin of American Meteorological Society, **April 2012**, 485-498.

24. M.R. Haylock, N. Hosfra, A.M.G. Klein Tank et al., Journal of Geophysical Research, **113**,D20119 (2008)

25. E. C. Mannashardt-Shamseldin, R. L. Smith, S. R. Sain et al., The annals of Applied Statistics, **4**, 484-502 (2010)
# European snow load map – past and present

Jerzy Antoni Żurański¹ and Grzegorz Kimbar¹a

¹Instytut Techniki Budowlanej, ul. Filtrowa 1, 00-611 Warszawa, Poland

**Abstract.** For few decades there has been a natural trend for harmonization of structural codes in the regional and global scales. An example of this trend is the elaboration of Eurocodes. One of the main items in the Eurocode for snow actions is a common map of the ground snow load in Europe. The paper deals with old values of the ground snow load in the Central Europe as well as the present state and use of maps given in Annex C to the Eurocode 1. It is stated that the present European maps of snow loads have been used in limited range. An attempt is made to explain the reason of such situation and to present possible further works.

## 1 Introduction

The extensive work was done to elaborate European maps of snow loads when preparing the Eurocode for snow loads in the frame of the European Committee of Standardisation (CEN) [6]. On the base of data from 18 countries eight snow load regions were established, crossing the borders of countries. Results of the work in the form of maps and dependence of the snow load on the ground on the altitude above sea level constitute the Annex C to the Eurocode for snow loads [1]. Fifteenth or even more years later it would be useful to review how they have been used in practice. Another purpose of the paper is to summarize the gained experience and to plan further works. They are necessary because during the time which has gone the new countries have joined the European Committee for Standardisation.

## 2 Old values of snow load in Europe

For a long time snow load, which has to be taken into structural calculations, have been established on the base of everyday, winter weather observations. An example of such an approach can be found in the book by Thullie [7]: "In our countryside snow layer on roofs is usually not higher than 0.6 m. Taking this value and assuming that snow is approximately 8 times lighter than water is, we receive [...] a round value of snow load 80 kg/m²". It was rounded from 75 kg/m² (0.75 kN/m²) and it is worth to mention that this value has been used in Central Europe for many years, even until II WW. It should be also pointed out that it was the value of snow load on flat roofs. Similar approach was reported by Holzer in his historical article [3].

Probably two turning points of snow load definition in Europe can be seen: the first – just after the II WW and the second one – 50 years later when the Eurocode for snow load was elaborated. These two turning points show three periods in the history of snow load maps in Europe: before IIWW, after it and since the elaboration of Eurocodes.

Above mentioned value, recommended in XIX century, was applied in several countries of Central Europe to the end of II WW. Next a distinct difference was seen: in the Central – East European countries dominated by the former Soviet Union snow load values were established according to the Soviet code. The snow load values in that code were smaller than older ones in subordinated countries. For example in former German Democratic Republic the value of 0.5 kN/m² [8] was much smaller than 0.75 kN/m² in Germany, still in force in Federal Republic of Germany in 1976 [1]. In Polish standards snow load for the first zone, the majority of territory, was successively diminished from 0.8 kN/m² established before II WW to 0.5 kN/m² in 1960s [11]. By the end of 1970s this value was identified as having only 5 year return period. In other countries were maritime climate occurs, e.g. in France, lower values were usually applied. Two categories were distinguished: normal and exceptional snow loads [6].

In early codes and recommendations, as mentioned above, there was no difference between snow load on flat roofs and on the ground. Probably the first such a difference was introduced in the ISO standard [5] based on the Canadian experience. In 1970s the probabilistic analysis was undertaken, which was influential for snow load values in European countries by the end of XX century.

## 3 Snow load maps in Eurocode

For snow load map to the Eurocode a new approach was applied, proposed by Gränzer [3].[4] As a result the Annex C to the Eurocode 1 was elaborated [2] [8]. On the base of data from 18 countries eight snow load regions were established, crossing the borders of countries. Values of the snow load on the ground were reduced to the sea level and given also as a function of the snow zone number. Results of the work in the form of maps and parabolic dependences of the snow load on the

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

<footnote>a Corresponding author: g.kimbar@itb.pl</footnote>
ground on the altitude above sea level constitute the
Annex C to the Eurocode for snow loads [2].

<figure><img src="image_140.png" /><figcaption>**Figure 1.** Base snow load values in zones according to combined national anexes in Alpine region.</figcaption></figure>

<figure><img src="image_141.png" /><figcaption>**Figure 2.** Base snow load values in zones according to Appendix C of [1] in Alpine region.</figcaption></figure>

## 4 Present uses of snow load maps

Eighteen National Annexes to the Eurocode for snow loads and two separate maps are now available from the following countries: Belgium, Czech Republic, Denmark, Finland, France, Germany, Holland, Iceland, Italy, Lithuania, Luxembourg, Poland, Romania, Slovakia, Spain, Sweden, Switzerland and UK (together with a map of Ireland). Among them Czech Republic, Iceland, Lithuania, Poland, Romania and Slovakia were not taken into account when preparing maps of snow regions to the Eurocode.

Maps given in Annex C to the Eurocode for snow action have been conceived as a helpful material for national code-writers. However, the use of these maps is rather limited (Fig. 1 and 2). The snow load zones in National Annexes differ from these given in Eurocode. Only in two countries zone numbers are explicitly introduced into the formula for characteristic value of snow load, in the UK and Ireland.

In three countries (Germany, Italy and Switzerland) the parabolic formula for the dependence of snow loads

on altitude is used, according to recommendations in
Annex C. In other countries linear dependence is used.

Almost in all countries the minimum value of altitude is given. Below this altitude the value of snow load is constant. It seems that detailed dependence of snow load on low altitudes is not preferable. These minimum values range from 100 m to 400 m, even to higher altitudes in some places (as in Romania). Differences in approaches to the snow load value become clearly visible at administrative borders. Different base values and altitude formulae can result in step difference on borders even by a factor of two (example on Figure 3).

Hence, it may be stated that the present European maps of snow loads has been accepted in limited range.

The question arises – why? It may be said that probably the main reason is the devotion to the traditional approaches in particular countries. Another question arises: how to improve this situation and to harmonise national maps with neighbours?

<figure><img src="image_142.png" /><figcaption>**Figure 3.** Characteristic ground snow load on the Germany-France border according to national anexes of [1].</figcaption></figure>

## 5 Harmonisation of European snow loads

The first step may be a short comparison of snow loads at
borderline between several countries in Central West and
Central East Europe. Results are presented in the Table 1

As it can be seen there are discrepancies between
snow loads at borderlines in particular countries. Only in
one case there is the same value, 0.65 kN/m²: in French
zone C1 and German zone 1, however, for different
altitudes, 200 m in France and 400 m in Germany. If
altitude is 400 m, the same in both countries, then snow
loads are different, 0.85 kN/m² in France and 0.65 kN/m²
in Germany (Fig. 3). It would be desirable that
harmonization would be possible on the base of bilateral
agreements between neighbouring countries. As en
example of harmonisation of snow load may serve a map
of snow loads in Poland and its eastern neibourgs (Figure
4). However, it should be told that snow load zones of
Polish eastern neighbours were harmonised with others
on their own accord. Finally, it has not be the
harmonisation with the same values of snow loads. It
would be sufficient if the value would differ not much.
<table><caption>**Table 1.** Snow loads along borderlines in Central West and Central East Europe</caption><thead><tr><th>Country</th><th>A<sub>min</sub>, m</th><th>s<sub>k</sub></th><th>s<sub>k</sub></th><th>A<sub>min</sub></th><th>Country</th></tr></thead><tbody><tr><td>Belgium</td><td>100</td><td>0.5</td><td>0.45</td><td>200</td><td>France</td></tr><tr><td>Belgium</td><td>100</td><td>0.5</td><td>0.65</td><td>200</td><td>France</td></tr><tr><td>Belgium</td><td>200</td><td>0.62</td><td>0.45</td><td>200</td><td>France</td></tr><tr><td>Belgium</td><td>200</td><td>0.62</td><td>0.65</td><td>200</td><td>France</td></tr><tr><td>Belgium</td><td>100</td><td>0.50</td><td>0.85</td><td>285</td><td>Germany</td></tr><tr><td>Belgium</td><td>285</td><td>0.72</td><td>0.85</td><td>285</td><td>Germany</td></tr><tr><td>Belgium</td><td>as above</td><td>as above</td><td>0.70</td><td>Not given</td><td>Holland</td></tr><tr><td>France A1</td><td>200</td><td>0.45</td><td>0.85</td><td>285</td><td>Germany z.2</td></tr><tr><td>France B1</td><td>200</td><td>0.55</td><td>0.85</td><td>285</td><td>Germany z.2</td></tr><tr><td>France C1</td><td>200</td><td>0.65</td><td>0.85</td><td>285</td><td>Germany z.2</td></tr><tr><td>France C1</td><td><b>200</b></td><td><b>0.65</b></td><td><b>0.65</b></td><td><b>400</b></td><td>Germany z.1</td></tr><tr><td>France C1</td><td><b>400</b></td><td><b>0.85</b></td><td><b>0.65</b></td><td><b>400</b></td><td>Germany z.1</td></tr><tr><td>Denmark</td><td>All</td><td>0.90</td><td>0.85</td><td>285</td><td>Germany</td></tr><tr><td>Germany z.1</td><td>400</td><td>0.65</td><td>0.70</td><td>Not given</td><td>Holland</td></tr><tr><td>Germany z. 2</td><td>285</td><td>0.85</td><td>0.7</td><td>300</td><td>Poland z.1</td></tr><tr><td>Germany z. 2</td><td>285</td><td>0.85</td><td>0.9</td><td>300</td><td>Poland z.2</td></tr><tr><td>Germany z.3</td><td>255</td><td>1.10</td><td>0.7</td><td>300</td><td>Poland z.1</td></tr></tbody></table>

<figure><img src="image_143.png" /><figcaption>**Figure 4.** Map of the characteristic values of snow load on the ground in Poland and its neighbours</figcaption></figure>

## 6 Climate change and snow load

It may be stated that all predicted values of snow load in Europe have been calculated on the base of data measured during last several decades. These elaborations have been done with silent assumption about constant and the same climatic conditions in the past and in the future. Now a question arises how to take onto account the possible climatic change. A simplest answer would be: to look at the annual run of data. An example is shown in Figure 5. If the data from 1950/51 to 1994/95 are taken into account the decreasing of snow load is obvious.

<figure><img src="image_144.png" /><figcaption>**Figure 5.** Annual run of the yaerly (winter) maxima of the ground snow load at meteorological station Rzeszów Jasionka (Poland).</figcaption></figure>
Based on this period of observation one may come to
the conclusion that in 2019 there will be no snow in
Rzeszów. However, the maximum value from the period
1950 – 2010 was measured in 2005.

This example is unique from Poland. At other
meteorological stations there are different trends of
annual run of snow, increasing or decreasing of the
annual maxima snow load but not so defined as in
Rzeszów in 1950 – 1995. The main conclusion which
may be drawn is that longer period of observation is
necessary for a forecast of the influence of climate
change on the structural snow loads.

7 Concluding remarks

Initial simplistic European approach to snow load value
was substantially enhanced in the course of XX century
with a vast body of field measurements and use of
statistical methods. This approach found its most recent
conclusion in formulation of Annex C to the Eurocode 1.
The aim of this document was to provide national
regulatory bodies with guidance in the process of
preparing local snow load maps, and by this mean to
promote harmonisation in the CEN member countries.
Nevertheless, the adoption of this document happened to
be sparse. This situation yields inconvenient and
irrational inconsistencies, especially at borders as well as
complicates usage of the Eurocode 1 on an international
level. Improvement of this state of affairs is at hand,
because apparently there are no fundamental obstacles
other than local customs persisting in the engineering
community. The question remains if the advantage of
harmonisation is worth breaking those already established
customs.

References

1. EN 1991-1-3:2003 Eurocode 1 – Actions on structures. – Part 1-3: General actions – Snow load

2. M. Gränzer, Zur Festlegung der rechnerischen Schneelasten. „Bauingenieur", 58 (1983) 1-5

3. M. Gränzer, New ground snow loads for a European standard. Proc. of the Fourth Int. Conf. on Snow Engineering, Balkema (2000)

4. S. M. Holzer, Kleine Geschite der Schnee- und Windlast im 19.Jahrhundert, Bautechnik 83, no 11, (2006)

5. ISO 4355:1980 Basis for design of structures – Determination of snow loads on roofs (1980)

6. Régles definissants les effets de la neige et du vent sur les constructions et annexes (Régles NV 65 révisées 1967 et annexes), Paris (1968)

7. L. Sanpaolesi, Scientific support activity in the field of structural stability of civil engineering works.

Snow loads. Final Report. University of Pisa,
September (1999)

8. TGL 32274/05 Lastannahmen für Bauwerke. Schneelasten. DDR, Berlin 1976

9. M. Thullie, Handbook of building static for engineers, architects and students (in Polish), Lwow (1886)

10. J.A. Żurański, Snow Loads on Roofs in East European Standards and Codes of Practice. - First International Conference on Snow Engineering, Santa Barbara 1988, CRREL Special Report 89-6, 1989

11. J.A. Żurański, A.Sobolewski, Snow loads in Poland (in Polish). Instytut Techniki Budowlanej, Warszawa 2009
# Some ongoing researches to improve codified structural design under snow loads in China

Feng Fan<fnref n="1" />, <fnref n="2" />, a, Huamei Mo<fnref n="1" />, Qingwen Zhang<fnref n="1" />, <fnref n="2" />, Guolong Zhang<fnref n="1" /> and Mengmeng Liu<fnref n="1" />

**Abstract.** Roofs in cold regions with heavy snowfalls could experience unbalanced or non-uniform snow loads caused by snow drifting or sliding, this could cause significant damage or collapse to the roof. The unbalanced snow load patterns recommended by the Chinese design code are originally based on that for Union of Soviet Socialist Republics (USSR) and have been calibrated several times, but documents to support such codifications are unclear. Besides, the provisions are mainly specified for regular-sized roofs but not long-span roofs. With the growing of applications of long-span roofs in recent years in China, there are practical demands for investigations of unbalanced snow loads on them since such investigations are rarely available. Therefore, an investigation of the unbalanced snow loads on roofs, especially for long-span roofs, is warranted. This paper summarizes some of the research activities on this subject at Harbin Institute of Technology (HIT), including statistics analysis of meteorological data, at-site measurements of snow depth and snow density on flat roofs, wind tunnel tests of model buildings, as well as numerical simulations of snowdrifting. Some of the works that are to be conducted in the near future are also presented and discussed.

## 1 Introduction

Snow load, which is usually unbalanced under the action of wind, is one of the dominant live loads that should be considered carefully in the structural design for buildings in cold regions with severe winter climate, especially for long-span roofs, which are sensitive to unbalanced loads. Heavy snow loads have caused failures of roof structures for many buildings; a very pathetic example of such failures is the collapse of the Katowice Fair Building in Poland, which caused 65 fatalities and 130 injuries [1].

The unbalanced snow load patterns for regular roofs are specified by the Chinese Load Code for the Design of Building Structures (GB 50009-2012) [2]. However, the provisions in the earliest version of this code (published in 1954) were based on that for Union of Soviet Socialist Republics (USSR). Although the specifications have been calibrated for 5 times since then (in 1958, 1974, 1987, 2001 and 2012, respectively), documents supporting such calibrations are rarely available.

Besides, the specifications for unbalanced snow load patterns in the code are mainly for roofs with regular sizes but not for long-span roofs. As one can imagine, the wind field around the long-span roofs may differ from that for regular-sized roofs, which will result in different patterns of unbalanced snow loads. On the other hand, with the growing of economy and societal needs, many structures with long-span roofs such as airport terminals, railway stations, stadiums and exhibition halls have been constructed in China in recent years. The safety of these structures are of public concern since the buildings' costs are usually high and, more importantly, the number of people accommodated by these buildings is relatively large, many fatalities could be caused in the case of roof collapses.

Consequently, it's of great interests for one to investigate the distribution of snow loads on roofs, especially for long-span roofs, to improve the codified

structural design under snow loads in China. Previous studies in recent years have conducted some related investigations on this subject [3-5], but their investigations were mainly focused on specific buildings, which may not be applicable other buildings. Therefore, a more general study is expected.

Techniques for investigations of roof snow loads include at-site measurements, wind tunnel tests (or water flume tests) and numerical simulations. For sloped roofs, however, at-site measurements are usually unrealistic because of the danger of operating on the wet, slippery roof in a freezing and windy environment, and that remote measurements are not always available. When real snow particles are not available, wind tunnel tests can be considered only for flat roofs as well, because model particles, such as silica sands, can hardly stay on sloped roofs, which is not the case for the real non-disturbed snow. Tests carried out in climatic wind tunnel where temperature is controlled, such as those described by [6, 7], is not always available. In this context, numerical simulation techniques become prior to the other two approaches in the investigation, although the validity of the numerical model should firstly be verified based on at-site measurements and/or wind tunnel tests before its application.

In this paper, some of the research activities on the investigations of unbalanced snow loads on roofs at Harbin Institute of Technology (HIT, Harbin, China) are summarized, these include statistics analysis of meteorological data, at-site measurements of snow depth and snow density on flat roofs, wind tunnel tests of snowdrift around model buildings, as well as numerical simulations of snowdrifting. Some of the works that are to be conducted in the near future is also presented and discussed.

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

<footnote n="1">¹School of Civil Engineering, Harbin Institute of Technology, Harbin, 150090, China</footnote>

<footnote n="2">²Key Lab of Structures Dynamics Behavior and Control of the Ministry of Education (Harbin Inst. of Tech.), Harbin, 150090, China</footnote>

<footnote>a Corresponding author: fanf@hit.edu.cn</footnote>
## 2 Ongoing research programs focused on snow loads at HIT

### 2.1 Statistics analysis of meteorological data

Ground snow load is the basis for calculation of roof snow loads in design codes. It's indicated that insufficient code provisions are one of the main observed causes of structural damage [8]. Consequently, accurate estimation of ground snow load is of important significance for the design of roof snow loads.

Basic ground snow load in the Chinese design code [2] is defined as the 50-year return period value of annual maximum ground snow load. As stated by the code, the annual maximum ground snow load is modelled as a Gumbel variate in the analysis of extreme ground snow load in China; and that annual maximum snow depth with a regional average snowpack bulk density is used to estimate the extreme ground snow load when measurements of ground snow load are not available [2]. However, no document is available to detail the estimation process of the recommended basic ground snow load in the code. To possibly improve the specification of basic ground snow load in China, meteorological data including snow depths and ground snow loads are collected from local meteorological agencies. Also assembled are surface wind speeds, precipitation and air temperature. Information for above mentioned meteorological data is summarized in table 1.

<table><caption>Table 1. Information of collected meteorological data</caption><thead><tr><td>Variable</td><td>No. of stations</td><td>Time periods</td><td>Measurement frequency</td></tr></thead><tbody><tr><td>Ground snow depth</td><td>734</td><td>1951-2010</td><td>Daily</td></tr><tr><td>Ground snow load</td><td>733</td><td>1999-2008</td><td>5-daily</td></tr><tr><td>Wind velocity</td><td>839</td><td>1951-2013</td><td>Daily</td></tr><tr><td>Precipitation</td><td>839</td><td>1951-2013</td><td>Daily</td></tr><tr><td>Air temperature</td><td>839</td><td>1951-2013</td><td>Daily</td></tr></tbody></table>

With these data available, several topics could be discussed. Studies that are being conducted by the authors and being considered to be conducted in the near future include the probabilistic distribution model of annual maximum ground snow depth, estimation of regional average snowpack bulk density, and relation between snow density and environmental variables such as air temperature and precipitation.

### Probability model for annual maximum snow depth

As mentioned above, the annual maximum ground snow load (or snow depth) is modelled as a Gumbel variate in the Chinese design code [2]. This is consistent with the practice in Canada [9, 10], Europe [11, 12], and in Japan [13, 14] (the Gumbel distribution is fitted only to the upper 1/3 of all the data plotted on Gumbel probability paper in Japan). However, Log-normal distribution is adopted to estimate the extreme ground snow water equivalent (SWE) in the US [15, 16]. The inconsistency

of adopted probability model for annual maximum snow load (or SWE, or snow depth) is likely due to the different climate conditions in different regions. An analysis of annual maximum ground snow loads for a northeastern province in China showed that, when compared to Gumbel distribution, Log-normal distribution is preferred for most of the stations in that province [17]; this is inconsistent with the recommendation in the Chinese design code. Considering this inconsistency and the availability of up to 60 years of daily ground snow depth data for more than 700 stations nationwide, it is of great interest for one to examine the preferred probability model for annual maximum snow depth in China.

Preliminary investigation shows that among 330 considered stations, Log-normal distribution is preferred by 186 of them, and that Gumbel is preferred by the remaining 144 stations. This could be clear evidence that if a single probability model should be used in the analysis of extreme ground snow load, Log-normal distribution seems to be the better choice. Further investigation is needed to verify this finding and to evaluate the influence of switching the probability model from Gumbel to Log-normal on the estimated 50-year return period value of annual maximum ground snow load.

### Snowpack bulk density

Snowpack bulk density is a random variable with significant variability [15], it's a complex function of snow depth, snow temperature, snow deposition history, and the initial density of the individual snow layers [18]. Attempts had been made by [18] to relate the snowpack bulk density with snow depth, snow age, and snow classes, where snow classes has taken into account the effects of climate (temperature and wind), the relation is written as [18]:

$$ \rho_b=( \rho_{\max} - \rho_0)[1-\exp(-k_1 \cdot h - k_2 \cdot \text{DOY})]+ \rho_0 \quad (1) $$

where $\rho_{\max}$, $\rho_0$, $k_1$, and $k_2$ vary with snow class, DOY runs from -92 (1 Oct.) to +181 (30 Jun.) with no 0 value. With the data shown in table 1 available, a model similar to Eq. (1) but applicable to China could be developed in the near future, this is one of the objectives of the ongoing programs at HIT.

In the Chinese design code, it's stated that a fixed regional average snow density could be used to convert annual maximum snow depth into ground snow load in the case where direct measurement of ground snow load is not available [2]. The National Building Code of Canada (NBCC) also uses average seasonal snowpack density to convert snow depth into ground snow load [9], while a nonlinear relation between 50-year return period value of ground snow load and 50-year return period value of ground snow depth is adopted in the US [16, 19].

The recommended regional average snow density by the Chinese design code is shown in figure 1. It can be seen from figure 1 that the country is divided into 6 zones for the assignment of regional average snow density. In this manner, there is sharp increase or decrease for the
snow density on the boundaries between the 6 zones, resulting in larger differences in the estimated extreme ground snow loads for locations that are close enough to each other but belong to different zones. Consequently, a more detailed and smooth zoning map for regional average snow density is desired to reduce such effect. Efforts are being made by the authors on this subject.

<figure><img src="image_146.png" /><figcaption>**Figure 1.** Regional average snow density recommended by the Chinese design code</figcaption></figure>

Alternatively, a nonlinear relation between the 50-year return period value of ground snow load and 50-year return period value of ground snow depth, such as the one adopted by ASCE load standard[16], could be considered to aid the conversion between the two extreme variables. However, 10 years of ground snow load data is insufficient to achieve a reliable estimate for the 50-year return period value of ground snow load, thus, development of such nonlinear relation could be considered only when more data is made available for ground snow loads in the future.

## 2.2 At-site measurements for unbalanced snow load profile

Unbalanced roof snow load, caused by snow drifting or snow sliding, is the main contribution to snow-induced structural damage. For this reason, special attention are paid to the unbalanced snow load profile by building codes such as those in China [2], Canada [10], US [16] and in Europe [12].

As indicated by the 1954 version of the Chinese design code, the snow provisions were mainly based on that for USSR. Although it's emphasized that the provisions were made by considering the observations in China, no detailed reports are available to introduce the technical supports for such considerations. For example, it's mentioned by [20] that concurrent measurements of snow depth on roofs and on ground were conducted in Heilongjiang, Liaoning, Inner Mongolia, Xinjiang and Hunan provinces in 1979, but no reports related to this program were found by the authors.

In contrast, there are sufficient case histories in the US and Canada to aid the calibration of snow provisions in building codes there. For example, 28 case histories were developed in [21] to investigate unbalanced snow load on gable roofs in the US, and that snow provisions for stepped roofs in the ASCE load standard were based on more than 350 case histories [22, 23].

Snow provisions in National Building Code of Canada were also based on sufficient data measured on roofs [24-26].

In this context, to collect data for the distribution of roof snow loads, at-site measurements of snow depth and snow density have been carried out on two flat roofs and on open ground in the campus of HIT since the snow season of 2010-2011 (see figure 2). Figure 3 shows an example of measured snow depth profile for one of the cross sections selected in the measurements. It can be seen from this example that snow tend to deposit in the location where elevation changes on the roof, the deposited snow can reach the top of the parapet while the snow depth in the middle part of the section is very limited. This observation is common among the measurements up to now. Quantitative analysis for the relation between the deposited snow depth, ground snow depth, parapet height and wind speed could be expected when sufficient measurements become available, this will help one to understand the unbalanced roof snow loads on stepped roofs, and behaviour of snow in the wind flow.

<figure><img src="image_147.png" /><figcaption>**Figure 2.** Locations where at-site measurements of snow depth and snow density were conducted</figcaption></figure>

<figure><img src="image_148.png" /><figcaption>**Figure 3.** Example of measured snow depth profile</figcaption></figure>

While snow depth on flat roofs or on stepped roofs can be measured directly by rulers by hand, those on gable roofs and arched roofs (in general, sloped roofs) cannot be done in the same way because it's unlikely for one to stay on such roofs safely to conduct the measurement. Consequently, sloped roofs are not included in the measurement program mentioned above. Considering the fact that long-span roofs are usually designed to be arched (e.g., cylindrical or spherical), and that snow loads on such roofs is usually of particular interest, it's desired by the authors to include sloped roofs into the measurement program in the future. A possible way to reach this purpose is to employ photogrammetric techniques using Unmanned Aerial Vehicle (UAV). An UAV have been equipped to HIT (see figure 4) and efforts are being made to use this equipment to measure snow depth. However, issues such as influence of low air temperature on the UAV's operation and strong reflection of the snow surface to the camera are still to be solved
before such measurement can be conducted. Once this equipment is made feasible for the measurements, snow depths on sloped roofs, especially on long-span roofs, will be conducted immediately and will last for years.

<figure><img src="image_149.png" /><figcaption>**Figure 4.** UAV equipped for the purpose of photogrammetric measurement of roof snow depth</figcaption></figure>

Sufficient measurements of roof snow depth for various types of roofs are expected at the end of this program. These data will be used as the basis to examine the validity of snow provisions in the Chinese design code and subsequently possible suggestions for the calibration of the provisions.

## 2.3 Wind tunnel tests

Wind tunnel modelling is one of the most effective ways to investigate the role of wind on snow accumulation. Investigations of snowdrift around/on buildings have been carried out in several previous studies. For example, Iversen [27] conducted a series of experiments on scale models of snow fences in the Iowa State University environmental wind tunnel using glass sphere as model particles. Isyumov and Mikitiuk [28] examined the drifting and redistribution of snow on large-area two-level flat roof by carrying out wind tunnel model test by using bran as model particles, while artificial snow was used by Delpech and Gandemer [6] to simulate snow drifting around buildings for Antarctic Concordia research station. Nevertheless, similarities between prototype and model are still big challenges for wind tunnel tests of snowdrifting today. A reliable modelling depends on similitude criteria, including geometric, dynamic and kinetic criteria [6]. Yet, it's generally accepted that fully correct, reduced scale modelling of snowdrifting is not possible in practice [29], consequently, compromises need to be made in the modelling. It's indicated by Anno [30] that it would be allowable to abandon Froude scaling if the shape of a snowdrift depends on the shape of wakes and the wind speed at the reference height, but this is not the case if the shape of snowdrift depends on the saltation trajectory. Kind [29] recommended to use high density particles in air in the modelling where time scaling and quantitative accuracy are more important than convenience, this is based on the observation that saltation dynamics of high density particles in air are dynamically similar to those of real snow in air.

While snowdrift around buildings or snow fences are relatively widely investigated by wind tunnel tests (e.g.,

[6, 27, 30]), investigation of snowdrift on building roofs is rarely available, especially for sloped roofs, this is due to the restriction of model particles' properties. Model particles such as sand, sawdust, glass spheres do not reproduce all inter-particle forces [6], this would result in the difficulty to model the snow on sloped roofs, even when the angle of repose for the model particle is equal to that for real snow. Considering this limitation, wind tunnel tests are not conducted directly for snowdrift on building roofs, but for snowdrift around standard model: cubic building model. The test results of snowdrift around cubic model will then be used as the benchmark to verify the validity of numerical simulation model that will be presented in the next section, the investigation of snowdrift on roofs will primarily rely on the numerical simulations thereafter.

For the wind tunnel test of snowdrift around cubic building model, a prototype with sizes of 1m×1m×1m was installed on an open ground in the campus of HIT and the snow depth around it was measured after each snowfall event. After conducting the measurement, the snow was cleaned to make sure there is no snow remained around the cube before next snowfall event. Figure 5 shows the distribution of snow around the cube after one of the snowfall events in the 2014-2015 season. The lateral erosion is clearly showed in this figure. However, there is no distinct deposition of snow in the windward or leeward side of the cube, this observation is different from that made by Oikawa et al. [31] in Sapporo, Japan and by Beyers et al. [32] in Antarctica; this is probably caused by the differences of snow properties and wind conditions in different regions. The measurement will keep going to collect as many cases as possible, and those cases with relatively steady wind direction will be selected as the prototype for the wind tunnel tests in the future.

<figure><img src="image_150.png" /><figcaption>**Figure 5.** Snow accumulation around the cube prototype</figcaption></figure>

Some preliminary investigations have been conducted by the authors in a wind tunnel of open circuit type with a cross section of 50cm × 50cm. Figure 6 shows the setup for the tests. In these tests, approaching wind is set uniform and the floor of the wind tunnel is covered by a uniform layer of model particles, no feeding or supply of particles is made during the experiments. The scaling of the test's parameters is carried out mainly by considering the criteria given by Eq. (2) and Eq. (20) in [33] and model particles are selected as bicarb soda powder.
Results showed that there are discrepancies between the tested model and its prototype. This is mainly because 1) the wind for the prototype is turbulent and the direction is varying during the snowfall, this cannot be modelled in the wind tunnel tests; and 2) snow particles are in falling movement in the prototype by are settled in the tests. Further improvements are needed to address these issues to reduce the discrepancies.

<figure><img src="image_151.png" /><figcaption>**Figure 6.** Schematic view of the wind tunnel test's setup</figcaption></figure>

Considering the fact that the wind condition is not controlled in the prototype experiment, one may desire to control the wind to reduce the uncertainties involved in the experiment. For this purpose, a 2×2 fan matrix is being constructed at HIT to pursue a stable wind direction and a desired wind speed for a specific duration in the prototype experiment (see figure 7 for illustration). Facilities such as temporal walls will be used to make sure the snow was uniformly distributed around the cube model during each snowfall event; the temporal walls will then be removed and the fans started. By doing in this way, the authors are expecting to make the prototype much more similarly modelled in the wind tunnel tests than it is in the current stage.

<figure><img src="image_152.png" /><figcaption>**Figure 7.** Schematic view of the proposed wind-controlled prototype experiment</figcaption></figure>

## 2.4 Numerical simulations

As mentioned above, at-site measurements and wind tunnel tests for snowdrift on sloped roofs is unrealistic due to some technical limitations such as the human safety and the satisfaction of similitude criteria. Consequently, the Computational Fluid Dynamics (CFD) technique may be an appropriate tool for the investigation of snow load on such roofs.

Application of CFD technique on simulation of snowdrift has been investigated for decades and was review by Tominaga et al. [34]. A typical model for simulating snowdrift using CFD technique is the one presented by Beyers et al. [32]. In this model, it's assumed that the air-snow system is highly dilute and that

only one-way coupling between the fluid and discrete phase is considered, i.e., the snow particles follow the continuous phase but do not influence it. The air-snow system is treated as turbulent, incompressible, Newtonian two-phase fluid flow and Eulerian approach was adopted. The inlet snow flux for the saltation layer was modelled as the empirical relations obtained by measured data by Pomeroy and Gray [35] and the one for suspension layer is based on Pomeroy and Male [36].

It's found that in most of the previous studies, the investigation is focused on snowdrift around buildings, for example, around snow fences [37, 38], cube models [32, 34], and building complex [39, 40]; those related to snowdrift on building roofs is quite rare. In this program, a model similar to the one presented by Beyers et al. [32] will be developed and to be validated by the wind tunnel tests mentioned earlier in the previous section, snowdrift on building roofs, especially on long-span roofs, will then be simulated by using the validated numerical model. These simulation results, supplemented by the at-site measurements, will help us to interpret the snow provisions in the design code and to possibly improve the provisions.

Besides, once the validity of numerical simulation model is verified, it could be a very useful tool for the design of buildings in snowy regions, especially for those with complicated shape and that shape coefficient is not given by the code. For one aspect, based on the simulation results, enhancements can be made to the structural elements where snowdrift tends to accumulate; second, the shape of the roof can be improved to avoid deposition of snow in the locations where snow accumulation is not expected.

## 3 Closing remarks

The snow load provisions in the Chinese design code are originally based on that for USSR and have been calibrated several times, but documents supporting such codifications are rarely available. Besides, with the growing of applications of long-span roofs in recent years in China, there are practical demands for investigations of unbalanced snow loads on these roofs since such investigations are rarely available. Therefore, an investigation of the unbalanced snow loads on roofs, especially for long-span roofs, is warranted. In this paper, some of the research activities on this subject at HIT, including statistics analysis of meteorological data, at-site measurements of snow depth and snow density on flat roofs, wind tunnel tests of model cube, as well as numerical simulations of snowdrift on building roofs are summarized. The programs are supposed to be finished at the end of 2018 and a detail report related to the interpretation to the snow load provisions in the Chinese design code and suggestions to possible calibration of the provisions is expected at that time. For this objective, several key techniques problems are to be solved or improved in the near future, these include:

1) Development of a more detailed and smooth zoning map for regional average snow density;
2) Photogrammetric measurement of roof snow depth using UAV installed with high-resolution camera;

3) Development of acceptable wind tunnel tests for cube model that are in good agreement with its prototype; and 4) Development and validation of CFD numerical model for simulating snowdrift around cube model.

## References

1. A. Biegus, K. Rykaluk, Eng. Fail. Anal., **16**, 1643 (2009)

2. GB 50009-2012, *Load code for the design of building structures* (China Architecture & Building Press, Beijing, 2012) (in Chinese)

3. X.Y. Zhou, M. Gu, X.F. Li, J. Build. Struct., **29**, 7 (2008) (in Chinese)

4. X.Y. Sun, C.B. Hong, Y. Wu, J. Vib. Shock, **33**, 36 (2014) (in Chinese)

5. Y. Li, X.F. Yuan, J. Build. Struct., **35**, 130 (2014) (in Chinese)

6. Ph. Delpech, P. Palier, J. Gandemer. J. Wind Eng. Ind. Aerod., **74-76**, 567 (1998)

7. L. Sanpaolesi, Commission of the European Communities DGIII/D-3 Final report phase II. (1999)

8. M. Holicky, M. Sykora. *Proc. Fifth Congress on Forensic Engineering*, ASCE. 444 (2009)

9. M. J. Newark et al. Can. J. Civil. Eng., **16**, 267 (1989)

10. NRC. *National Building Code of Canada*. Ottawa (2010)

11. L. Sanpaolesi, Commission of the European Communities DGIII/D-3 Final report phase I. (1998)

12. BSI. *Eurocode 1 -Actions on structures - part 1-3: General actions - Snow loads*. London (2003)

13. M. Izumi, H. Mihashi, T. Takahashi, *First International Conference on Snow Engineering*, Santa Barbara, California (1988)

14. AIJ. *AIJ recommendations for loads on Buildings*. Tokyo (2004)

15. B. Ellingwood, R. Redfield, J. Struct. Eng., **109**, 950 (1983)

16. ASCE. *Minimum Design Loads for Buildings and Other Structures*. Reston (2010)

17. H. M. Mo, F. Fan, H. P. Hong, Nat. Hazards, **77**, 543 (2015)

18. M. Sturm et al., J. Hydrometeorol., **11**, 1380 (2010)

19. R. L. Sack. J. Struct. Eng. **142** (1) (2015)

20. B.N. Hou, C.A. Wei, Metallurgical Buildings., (2), 61 (1982) (in Chinese)

21. M. O'Rourke, M. Auren, J. Struct. Eng., **123**, 1645 (1997)

22. M. O'Rourke, W. Tobiasson, E. Wood, J. Struct. Eng., **112**, 2080 (1986)

23. M. O'Rourke, C. DeAngelis, J. Struct. Eng., **128**, 1330 (2002)

24. D. A. Taylor, Can. J. Civil Eng., **6**, 85 (1979)

25. D. A. Taylor, Can. J. Civil Eng., **8**, 63 (1981)

26. D. A. Taylor, Can. J. Civil Eng., **19**, 59 (1992)

27. J. D. Iversen, J. Wind Eng. Ind. Aerod., **8**, 231 (1981)

28. N. Isyumov, M. Mikitiuk, J. Wind Eng. Ind. Aerod., **36**, 893 (1990)

29. R. J. Kind, Cold Reg. Sci. Technol., **12**, 217 (1986)

30. Y. Anno, Cold Reg. Sci. Technol., **8**, 241 (1984)

31. S. Oikawa, T. Tomabechi, T. Ishihara, J. Snow Eng. Japan, **15**, 283 (1999)

32. J. H. M. Beyers, P. A. Sunsbo, T. M. Harms, J. Wind Eng. Ind. Aerod., **92**, 725 (2004)

33. R. J. Kind, Atmos. Environ., **10**, 219 (1976)

34. Y. Tominaga, T. Okaze, A. Mochida, Build. Environ., **46**, 899 (2011)

35. J. W. Pomeroy, D. M. Gray, Water Resour. Res., **26**, 1583 (1990)

36. J. W. Pomeroy, D. H. Male, J. Hydrol., **136**, 275 (1992)

37. T. Uematsu et al., Cold Reg. Sci. Technol., **20**, 65 (1991)

38. M. Naaim, F. Naaim-Bouvet, H. Martinez, Ann. Glaciol., **26**, 191 (1998)

39. Y. Tominaga, A. Mochida, J. Wind Eng. Ind. Aerod., **81**, 273 (1999)

40. M. Beyers, B. Waechter, J. Wind Eng. Ind. Aerod., **96**, 1603 (2008)
# Specification of the design value of the ground snow load considering measurements of the snow height – part 1: single stations

Benjamin Czwikla, Michael Kasperskiª

Ruhr-Universität Bochum, Research Team EKIB. 44780 Bochum, Germany

**Abstract.** The straight-forward approach to the specification of the design value of the ground snow loads is based on measurements of the water equivalent. In Germany, this parameter however is only available for a few stations and the most recent decades. Based on a snow density model for each individual station, the additional years of observed snow heights can be used to increase the statistical stability of the estimated design value of the snow load on the ground.

## 1 Introduction

The straight-forward approach to the specification of the design value of the ground snow loads is based on measurements of the water equivalent. Since the acquisition of this observation variable is rather time-consuming not all meteorological stations provide this information. For stations which do provide the water equivalent, the observations usually cover only the recent decades, i.e. the majority of observation data provides only the height of snow cover. Using these data for the specification of snow loads requires a consistent model for the snow density ρs which translates the snow heights hs to snow loads s:

$$s = \rho_s \cdot h_s \quad (1)$$

For exploiting the observation years of snow heights of a single station, obviously it is sufficient to develop an individual snow density model for each station. For Germany, this step increases the number of observations on average by a factor of about 2, thus improving the statistical stability of the estimated design values. The individual models require simultaneous measurements of snow heights and water equivalents.

The number of German stations only measuring the snow heights is by a factor of about 14 larger than the number of stations measuring the water equivalent. Exploiting this huge amount of data requires a regional model of the snow density. This paper discusses the available models for the snow density, specifies required features and presents a new approach which is based on extreme value statistics for both the water equivalent and the snow heights.

## 2 Models for the snow density

The basic demand for any model probably deals with the number of parameters to be identified, i.e. the best approach follows the demand to have as few parameters as possible and as many parameters as required. Basically, snow density in the context of translating the height of a snow cover to snow loads means the average density of the snow pack. Then it seems to be reasonable that the model should be monotonic in regard to the

height of snow cover and probably should have a limiting value.

The simplest model, used e.g. in [1] is based on a linear approach and uses only two parameters:

$$\rho_s = a \cdot h_s + b \quad (2)$$

In Russia [2], the influence of the height is considered to be smaller by introducing the square root:

$$\rho_s = a \cdot \sqrt{h_s} + b \quad (3)$$

A more general approach introduces as third variable an unknown exponent for the height [1]:

$$\rho_s = a \cdot h_s^c + b \quad (4)$$

These models have in common that a limiting value has to be specified for the snow density.

In [3], a three-parameter model has been developed which describes the snow density between two limiting values ρ0 and ρ∞:

$$\rho_s = \rho_\infty \cdot \frac{h_{\text{ref}}}{h_s} \cdot \ln \left[ 1 + \frac{\rho_0}{\rho_\infty} \cdot \left( \exp \left( \frac{h_s}{h_{\text{ref}}} \right) - 1 \right) \right] \quad (5)$$

Finally, in [4], a four-parameter model is proposed which uses as additional parameter the age of the snow cover in days DOY:

$$\rho_s = (\rho_\infty - \rho_0) \cdot \left[ 1 - \left( \exp(-k_1 \cdot h_s - k_2 \cdot \text{DOY}) \right) \right] + \rho_0 \quad (6)$$

Other approaches aim in translating the snow height into snow loads only for a specific annual non-exceedance probability [2], e.g. in USA for a return period of 50 years or in Japan for a return period of 100 years.

## 3 New approach

The new approach extends the above idea of translating the snow height into a corresponding snow load for the complete trace of non-exceedance probabilities. This approach leads to an average probabilistic model of the snow density depending on the snow height.

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

<footnote>ª Corresponding author: michael.kasperski@rub.de</footnote>
In a first step, it is assumed that the traces of the annual non-exceedance probability follow the extreme value distribution type I (Gumbel distribution) which is given as:

$$F(x) = \exp\left[-\exp\left[-\left(\gamma + \frac{\pi}{\sqrt{6}} \cdot \left(\frac{x-m}{\sigma}\right)\right)\right]\right] \quad (7)$$

with m – mean value, σ - standard deviation, γ - Euler constant ≈ 0.5772

In the probability paper, both traces appear as straight lines; dividing the two linear equations leads to the following simple approach:

$$\rho_s = \rho_\infty - \frac{\Delta\rho \cdot h_{\text{ref}}}{h_s} \quad (8)$$

The limiting density $\rho_\infty$ depends on the standard deviations of the snow load and the snow heights:

$$\rho_{\infty} = \frac{\sigma_{s}}{\sigma_{h} \cdot g} \quad (9)$$

with g – gravitational acceleration ≈ 9.81 m/s²

The term $\Delta\rho$ is obtained as the weighted difference of the two mean values. For consistency in the dimensions, it is required to introduce a reference height:

$$\Delta\rho = \frac{\rho_{\infty} \cdot m_{h} - m_{s}}{h_{\text{ref}}} \quad (10)$$

In this approach, the reference height is an arbitrary value and can be assumed as $h_{ref} = 1$ m.

The individual traces of extreme snow heights and snow loads are based on the maximum observations during independent events with a continuous snow cover, i.e. the observations do not necessarily occur on the same day. The approach can either be based on yearly extremes or on individual events above a specific threshold value. For the latter ensemble, the non-exceedance probability per event is translated to the annual non-exceedance probability by using the following relation:

$$p_{\text{annual}}(x \le x_{\text{ref}}) = \sum_{K=0}^{\infty} p(K) \cdot p_{\text{event}}(x \le x_{\text{ref}})^K \quad (11)$$

with $p(K)$ – probability of having 0, 1, 2, 3, 4... independent snow events in a 12 month period
$p_{event}$ – non-exceedance probability per event

Basically, the annual probabilities in equation (11) follow the generalized extreme value distribution if the number of events per year follows a Poisson distribution and the non-exceedance probability per event corresponds to the Generalized Pareto distribution. If the annual probabilities follow the Gumbel distribution, then, the consistent model demands that the individual events follow an exponential distribution which is given as:

$$F(x) = 1 - \exp\left(-\frac{(x - x_{\text{th}})}{sc}\right) \quad (12)$$

with sc – scale parameter, $x_{th}$ – threshold value

The Poisson distribution is given as:

$$p(K) = \frac{\lambda^K}{K!} \exp(-\lambda) \quad (13)$$

with $\lambda$ – average number of events per 12 month cycle

The translation scheme of the probabilistic event-based model to the probabilistic annual model looks as follows:

$$\begin{aligned} m &= x_{\text{th}} + (\gamma + \ln \lambda) \cdot sc \\ \sigma &= \frac{\pi}{\sqrt{6}} \cdot sc \end{aligned} \quad (14)$$

For the identification of the ‘true’ model for the snow density, the fitting procedure should aim in unbiased estimates of the four describing parameters which are the mean values and the standard deviations of the extremes of the snow heights and the water equivalents, respectively. This obviously can be achieved by applying BLUE, i.e. the Best Linear Unbiased Estimator. Usually, for the Gumbel distribution, an alternative notation is used:

$$F(x) = \exp\left[-\exp\left[-\left(\frac{x-\mu}{\beta}\right)\right]\right] \quad (15)$$

The two parameters $\mu$ and $\beta$ are obtained as weighted sums of the observations $x_i$ which have to be sorted in ascending order:

$$\begin{aligned} \mu &= \sum_{i=1}^{N} x_i \cdot A_i \\ \beta &= \sum_{i=1}^{N} x_i \cdot B_i \end{aligned} \quad (16)$$

with N – ensemble size

The mean value and the standard deviation are obtained as follows:

$$\begin{aligned} m &= \mu + \beta \cdot \gamma \\ \sigma &= \frac{\pi}{\sqrt{6}} \cdot \beta \end{aligned} \quad (17)$$

Usually, the required weighting coefficients A and B are provided for a set of annual extremes, i.e. the non-exceedance probabilities are assumed to be equidistant. Respective tables for ensembles sizes up to N = 30 can be found in [5]. For larger ensembles sizes, the weighting coefficients can be obtained based on a Taylor series expansion [6].
Equation (11) leads to non-equidistant probability values. The required weighting coefficients can be achieved based on [6] by simply exchanging the individual p-values $p_i = i/(N+1)$ by the corresponding value from equation (11).

Alternatively, a least square fit in the Gumbel probability paper provides unbiased estimates if the appropriate plotting position as recommended by Gringorton [7] is used, i.e. the probability $p_i$ becomes:

$$p_i = \frac{i-0.44}{N+0.12} \qquad (18)$$

The target probability distribution occurs as a straight line. The demand that the sum of the squared differences between observation and theoretical value should be minimal leads to a simple system of two linear equations:

$$\left[ \begin{matrix} \sum_{i=1}^{N} x_i^2 & \sum_{i=1}^{N} x_i \\ \sum_{i=1}^{N} x_i & N \end{matrix} \right] \cdot \left[ \begin{matrix} c_1 \\ c_2 \end{matrix} \right] = \left[ \begin{matrix} \sum_{i=1}^{N} x_i \cdot y_i \\ \sum_{i=1}^{N} y_i \end{matrix} \right] \qquad (19)$$

The coordinate y in the probability paper is given as:

$$y_i = -\ln(-\ln(p_{i, \text{annual}})) \qquad (20)$$

Finally, the two unknown parameters m and σ are obtained from $c_1$ and $c_2$ as follows:

$$m = \frac{\gamma - c_2}{c_1}, \quad \sigma = \frac{\pi}{\sqrt{6}} \cdot \frac{1}{c_1} \qquad (21)$$

Once the model for the snow density has been obtained, the next step is to translate the additional observations of $h_s$ to corresponding snow loads s. Then, the design value can be estimated based on an increased ensemble.

Tentative target values of the annual exceedance probability of the design snow load on the ground are specified in [8] for different consequence classes. The design value $s_{des}$ then is obtained as follows:

$$s_{des} = m_s + s_{des, red} \cdot \sigma_s \qquad (22)$$

The reduced design value depends on the target exceedance probability:

$$s_{des, red} = \left( -\ln(-\ln(p_{\text{target}})) - \gamma \right) \cdot \frac{\sqrt{6}}{\pi} \qquad (23)$$

It is important to note that the increased ensemble size is still far away from infinity, i.e. any estimation of the design value is influenced by statistical uncertainties. In accordance to the demands in the EN 1990 – Basis of Design, the final extreme value statistics aim in a one-sided confidence interval of 75%. Since an estimation

based on BLUE leads on average to the true value, the corresponding one-sided confidence is only in the range of 50%. Therefore, an additional adjusting factor is required to increase the confidence in the final design values to the desired level:

$$s_{des} = m_s + s_{des, red} \cdot \sigma_s \cdot c_{ad} (p_{\text{target}}, N) \qquad (24)$$

The adjusting factor depends on the ensemble size and the target exceedance probability. It can be obtained for any estimation method based on simulations.

## 4 Derivation of the adjusting factor

For confined ensembles, the estimations of the mean value and the standard deviation of the extremes contain some uncertainties. Basically, the smaller the ensemble, the larger will be the random scatter in the estimated values. The randomness can be analysed based on simulations. The following results are based on 1,000,000 independent runs, assuming that the basic variable follows the Gumbel distribution with zero mean and a standard deviation of 1. Figure 1 shows the estimated statistical parameters for different ensemble sizes for a set of annual extremes.

<figure><img src="image_154.png" /><figcaption>**Figure 1:** Randomness in the estimated mean value and standard deviation for an ensemble of annual extremes following the Gumbel distribution for different ensemble sizes</figcaption></figure>

The design value is estimated introducing a target exceedance probability of 1/1000 per year. Figure 2 shows the probability distribution of the estimated design value for different ensemble sizes of annual extremes. The true target value is obtained from equation (19) as
$s_{des, red} = 4.936$. Additionally, the 90%-confidence interval is shown.

Since a reasonable two-sided confidence interval requires a very large ensemble size, which usually is not available, the basic strategy in design is to avoid under-estimations in case of actions and over-estimations in case of the resistance. Then, a one-sided confidence interval is used. Applying the target confidence recommended in the Eurocode requires shifting the traces of the non-exceedance probabilities to the right side with the basic aim that all traces go through the anchor point which is given with the true design value and the chosen error target. The latter is obtained as 1-conf, where conf is the confidence target. This shifting of the traces increases inevitably the probability of overestimating the design value which leads to uneconomic design. The higher the confidence demand, the higher is the probability of an uneconomic design. The 75%-target value is therefore a fair compromise to balance the demands to safety and economy.

<figure><img src="image_155.png" /><figcaption>**Figure 2:** Randomness in the estimated design value with $p_{target} = 0.999$ for an ensemble of annual extremes following the Gumbel distribution for different ensemble sizes and two-sided 90%-confidence interval</figcaption></figure>

The required adjusting factor can be obtained iteratively based on simulations. Figure 3 shows the respective results for ensemble sizes from 8 to 128. It is important to note that in this approach all annual extremes have to be larger than zero, i.e. this approach is recommended for snow climates which in each year show a continuous snow cover over the whole winter season.

For snow climates where in some individual years there is no snow at all and / or for snow climates which might have more than one independent period of

continuous snow cover per winter season, the approach based on equation (11) has to be used. Then, beside the actual number of observation years the average number $\lambda$ of independent periods with continuous snow cover has to be considered for the specification of the adjusting factor. The respective results are shown in figure 4.

<figure><img src="image_156.png" /><figcaption>**Figure 3:** Adjusting factor for the estimation of the design value for annual extremes following the Gumbel distribution, $p(s > s_{des}) = 0.001 / \text{year}$</figcaption></figure>

<figure><img src="image_157.png" /><figcaption>**Figure 4:** Adjusting factor for the estimation of the design value for extremes above a threshold value following the exponential distribution and for the number of events following the Poisson distribution for different average number of independent events $\lambda$, $p(s > s_{des}) = 0.001 / \text{year}$</figcaption></figure>

# 4 Example of application

One of the places in Germany with really long weather records is the station at Aachen. Observations for the height of snow cover start already in 1892. However, there are some gaps. The years 1904 and 1915-1929 are missing. Since 1950, water equivalents are sampled if the snow cover is at least 1 cm. The snow climate at Aachen is characterized by eventual mild winters without considerable snow. In the period from 1950 to 2006 there are 17 winters without snow. Altogether, in 56 winter seasons 83 independent events with a continuous snow cover are observed providing both the snow load and the snow height. Figure 5 shows the traces of the non-exceedance probability and the fitted theoretical model for the snow loads. In figure 6, the corresponding trace for the snow heights is shown.
The final model for the snow density for Aachen is then identified as:

$$ \rho_s = 161.238 - \frac{3.129 \cdot h_{ref}}{h_s} [kg/m^3] \quad (25) $$

<figure><img src="image_158.png" /><figcaption>Figure 5: Probability plot for the snow loads at Aachen (1950-2006)</figcaption></figure>

<figure><img src="image_159.png" /><figcaption>Figure 6: Probability plot for the snow heights at Aachen (1950-2006)</figcaption></figure>

Figure 7 compares the identified model of the snow density at Aachen to the observed values for the 83 independent events. It is important to note, that the snow density model has to extend to values which are in the range of the design value of the snow heights.

Based on the identified snow-density model, the observation years from 1892 to 1949 lead to 163 further observations, i.e., the ensemble size is almost three times larger. However, many early observations correspond to snow heights of only a few centimetres. Therefore, it is reasonable to introduce an individual threshold value which for Aachen can be assumed as $s_{lim} = 0.05$ kN/m². Then, the average number of independent events becomes 1.37 per twelve month cycle. The total number of events is increased to 129. Figure 8 shows the extended trace of annual non-exceedance probability for the observation period from 1892 to 2006.

<figure><img src="image_160.png" /><figcaption>Figure 7: Probabilistic model of the snow density (Aachen 1950 – 2006)</figcaption></figure>

<figure><img src="image_161.png" /><figcaption>Figure 8: Probability plot for the estimated snow loads at Aachen (1892-2006)</figcaption></figure>
The extended ensemble obviously contains an outlier which occurred in the winter of 1908. This value is almost twice as large as the second largest observed snow load in 1987. For the outlier, the corresponding observed snow height is 75 cm, which is 1.75 times larger than the second largest snow height of 43 cm which was observed in 1969.

The least-square fit therefore neglects this outlier. The final estimation of the design value requires an appropriate adjusting factor. For the initial ensemble, the adjusting factor is 1.131, for the increased ensemble it slightly reduces to 1.102. The final design value in accordance to equation (22) then becomes 1.10 kN/m². This value is smaller than the actual DIN-EN value which is 1.5 * 0.85 = 1.275.

For most stations, the records of the water equivalent start in the 1950s. However, a few stations show considerable shorter periods either due to a later starting year or due to gaps in the observation, e.g. Dresden with the starting year 1978 or Jena with missing data from 1980 to 1990. Therefore, the statistical uncertainties in the identified parameters may vary from station to station.

Furthermore, not all traces show similar well-behaved shapes as in the above example for Aachen. Figures 9 and 10 show the respective traces of the observed snow loads and snow heights for Schwerin from 1947 to 2006.

<figure><img src="image_162.png" /><figcaption>**Figure 9:** Probability plot for the snow loads at Schwerin (observation period 1947-2006)</figcaption></figure>

While the trace of the snow loads suggests a type II extreme value distribution, the trace of the snow heights seems to follow a type III extreme value distribution. The question arises, if the largest observed value in the snow loads really is an outlier.

Nevertheless, for the identification of the snow density model the largest observed snow load is ignored. The identified model for the snow density at Schwerin is shown in figure 11. The supposed outlier occurs with an extremely large density of more than 450 kg/m³ for a snow height of 30 cm.

<figure><img src="image_163.png" /><figcaption>**Figure 10:** Probability plot for the snow heights at Schwerin (observation period 1947-2006)</figcaption></figure>

<figure><img src="image_164.png" /><figcaption>**Figure 11:** Probabilistic model of the snow density for Schwerin (observation period 1947 – 2006)</figcaption></figure>

# 4 Versus a snow density model for Germany

In a first step, the data of 225 stations are analysed. Figure 12 shows the positions of the stations on a map of Germany. The minimum requirement for applying the new approach is 10 years of simultaneous observations of snow load and snow height.

Figure 13 shows the interdependency of the parameters of the snow density model. For the majority
of stations, the model predicts an increase of the snow density with snow height. Only about 4% of the stations show the opposite behaviour, i.e. for small snow heights the snow density is large, and the snow density decreases for increasing snow height. The results so far suggest that there is a non-linear decreasing trend for $\Delta\rho$ if $\rho_\infty$ increases.

<figure><img src="image_165.png" /><figcaption>**Figure 12:** Distribution of analysed German stations</figcaption></figure>

<figure><img src="image_166.png" /><figcaption>**Figure 13:** Interdependency of the parameters of the snow density model for 225 stations in Germany</figcaption></figure>

The results so far suggest that the snow densities do not depend on the altitude for stations below 250 m. This finding is in agreement with the basic setting of base values of the design snow load in the German snow load standard DIN-EN 1991-1-3/NA:2010-12. For Germany, three basic zones are specified. For zone 1 and 2, additionally regions 1a and 2a are distinguished. The minimum characteristic values are 0.65, 0.85 and 1.1 kN/m² for zone 1, 2 and 3. The appropriate values for the zones 1a and 2a are obtained by multiplying the

corresponding values of zone 1 or 2 with 1.2. For Germany, the partial factor for the snow load is 1.5

In figure 14 the identified design snow loads are plotted versus the altitude for altitudes below 250 m, for sake of simplicity however only for stations in zone 1 and 2. It is important to note that these estimates do not contain the adjusting factor to meet the target confidence, i.e. on average the results correspond to the true design values.

The results suggest that for zone 1 the minimum design value can be further reduced. For zone 2, many stations show a considerable larger design snow load than is specified in the German snow load code although the adjusting factor has not yet been considered. This comparison underlines the need of a consistent analysis of the German snow climate to meet the demands of safety and economy.

<figure><img src="image_167.png" /><figcaption>**Figure 14:** Improved design values for German stations with altitudes below 250 m in comparison to the specifications of DIN-EN 1991.1</figcaption></figure>

# 5. Summary and conclusions

The straight-forward approach to the specification of the design value of the ground snow loads is based on measurements of the water equivalent. These observations are available in Germany for about 550 stations, however only for the recent decades. Observations for the snow height on the other hand start for some stations already at the end of the 19th century. On average, the observation period of the snow height is twice the observation period of snow loads. Therefore, the statistical stability of estimated design values of the snow load can be increased by exploiting the additional data for the snow heights. This requires a model for the snow density

The new approach identifies the snow density as a probabilistic optimum, i.e. the approach demands that the traces of the annual non-exceedance probabilities of the observed snow heights and snow loads fall together if the optimum snow density model is applied. Assuming that both traces follow the extreme value distribution type I leads to a two-parametric model for the snow density.

It is important to note that the identification of the snow density model requires an unbiased estimator. For annual extremes, this unbiased estimator is provided with
the Best Linear Unbiased Estimator for the Gumbel
distribution. If the model is based on independent events
above a threshold value, new sets of BLUE-coefficients
are required. Alternatively, the estimation can be based
on Gringortons plotting position which also leads to
unbiased estimates.

The final estimation of the design value has to meet a
once-sided confidence which in the Eurocode is specified
as 75%. Since the increased ensemble, which combines
truly observed snow loads and estimated snow loads
based on observed snow heights, is still far away from
having a size close to infinity, adjusting factors are
required to push the confidence of the unbiased estimates
to the target value. These adjusting factors can be derived
based on simulations. The paper presents corresponding
figures for annual extremes and independent events.

Acknowledgement

Part of this work has been sponsored by the German
Research Foundation under the contract number
KA675/14-3. This support is gratefully acknowledged.

**References**

1. T. Jonas, C. Marty, J. Magnusson: Estimating the snow water equivalent from snow depth measurements in the Swiss Alps, Journal of Hydrology 378 (2009) 161–167

2. ISO 4355 (2013) Determination of snow loads on roofs

3. M. Gränzer, *Angaben von Schneelasten, geographisch nach Zonen gegliedert für den Eurocode "Lasten"*. Landesstelle für Baustatik, Baden-Württemberg. Landesstelle für Baustatik, Baden-Württemberg, Tübingen 1989 (in German)

4. J.L. McCreight, E.E. Small: Modeling bulk density and snow water equivalent using daily snow depth observations, The Cryosphere, 8, 521–536, 2014

5. N. Balakrishnan, P.S. Chan : Order statistics from extreme value distributions, II : Best linear unbiased estimates and some other uses, Communications in Statistics – Simulation and Computation, 1219-1246, 1992

6. A. Child, Balakrishnan, N. (2002) Series approximations for moments of order statistics using MAPLE, Computational Statistics & Data Analysis 38, 331-347

7. I.I. Gringorton : A plotting rule for extreme probability paper, Journal of Geophysical Research, 68, 1963, 813-814

8. ISO 22111 (2007): Bases of design of structures – General requirements
# Research on reliability of roof structures designed by Chinese codes

Shengguan Qiang, Xuanyi Zhouª, Ming Gu

*State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China*

**Abstract.** In this paper, the annual maximum ground snow loads of three representative cities in China are simulated by a snowmelt model based on available historical meteorological data, and six probabilistic models are chosen for describing them and are evaluated by K-S test, AIC and Q-Q plots. Additionally, basing on the six different probabilistic models for annual maximum ground snow loads, the reliability indices of steel roof members in the three cities are assessed. The results indicate that the fit of a probabilistic model to annual maximum ground snow loads, especially to the high quantiles, has great influence on reliability assessment of roof structures.

## 1 Introduction

Load and resistance factor design (LRFD) method has been widely used to building design. In LRFD method, structural safety is expected to achieve a targeted level by using different partial factors for resistance and different loads. However, the reliability of roof structures whose design is governed by snow loads is actually of much uncertainty. For instance, Takahashi and Ellingwood (2005) investigated reliability of roofs subjected to extreme snow load in Japan, they presented that steel roof members designed by the *Japanese Building Code* (JBC) appear to have a lower reliability against snow load compared to roof members designed by codes in the United States, Canada, or Western Europe, because that the partial factor for snow load and the snow density in JBC are not conservative [1]. Kozak and Liel (2015) studied the reliability of two types of steel roofs under snow loads designed according to ASCE 7-10, they found that the reliability of studied roof types is highly dependent on the location and seasonal snowfall patterns at the building location and the inclusion of serviceability deflection limits in design [2]. However, the influence on reliability assessment of the uncertainty of the probabilistic model for the annual maximum ground snow loads of a city has not been considered in their researches. In fact, there are maybe multiple underlying probabilistic models for annual maximum ground snow loads of a city, and different results of reliability assessment will be got if different probabilistic models are used to characterize annual maximum ground snow loads for assessing structural reliability.

In this paper, the annual maximum ground snow loads of three representative cities in China, which include Beijing, Nanjing and Hohhot, are simulated by a snowmelt model presented by Zhou et al. (2013) [3]. Then, three underlying probability distributions and three types of parameter estimation methods are combined into six different probabilistic models to fit the annual maximum ground snow loads, and K-S test, Akaike information criterion (AIC) and Q-Q plots are adopted as tools to evaluate the quality of these models. Finally, the reliability indices of steel roof members designed by Chinese load code GB50009-2012 [4] in the three cities are investigated on the basis of the six probabilistic models for snow loads, and the influence on reliability

assessment of different probabilistic models for snow loads is discussed.

## 2 Samples of annual maximum ground snow loads

Most weather stations in China give information on ground snow accumulation simply in terms of the depth of snow without simultaneous measurements of snow density [4]. Hence, a mean snowpack density has to be used to estimate ground snow load in Chinese load code GB 50009-2012. Thus, uncertainty of snowpack density which is an important source of the variability of ground snow loads has been ignored.

Zhou et al. (2013) presented a snowmelt model based on the mass-and-energy balance models, which can be used to simulate ground snow loads with available historical meteorological data, including daily average wind speed, temperature, relative humidity, and amount of snowfall or rainfall. With the help of the snowmelt model, the annual maximum ground snow loads of Beijing, Nanjing and Hohhot during 60 winters from November 1, 1951 to March 31, 2011 are simulated. Fig.1 shows the sample of simulated snow loads of Beijing, and Tab.1 lists the mean value and coefficient of variation (COV) of the annual maximum ground snow loads of the three cities. The coefficient of variation (COV) is defined by the ratio of standard deviation to mean value.

<figure><img src="image_168.png" /><figcaption>**Fig. 1.** Sample of annual maximum ground snow loads in Beijing</figcaption></figure>

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

<footnote>ª Corresponding author: zhouxytj@tongji.edu.cn</footnote>
<table><caption>Tab.1. Mean value and coefficient of variation</caption><thead><tr><th></th><th>Beijing</th><th>Nanjing</th><th>Hohhot</th></tr></thead><tbody><tr><td>mean (kPa)</td><td>0.14</td><td>0.23</td><td>0.30</td></tr><tr><td>COV</td><td>0.79</td><td>1.02</td><td>0.30</td></tr></tbody></table>

### 3. Structural design

In accordance with GB 50009-2012, LRFD method should be used to design a building structure. For the ultimate limit states, considering the combination of dead and snow loads on a roof, the design value of load effect $S_d$ is determined by formula (1), meanwhile structural resistance must meet formula (2).

$$S_d = \max\{1.2S_{Dk} + 1.4\gamma_L S_{Srk}, 1.35S_{Dk} + 1.4\gamma_L \varphi_c S_{Srk}\} \quad (1)$$

$$\gamma_0 S_d \leq R_k / \gamma_R \quad (2)$$

where 1.2 and 1.35 are partial factors for dead load when design is governed by variable load and dead load, respectively. 1.4 is partial factor for variable loads. $\varphi_c$ is the coefficient of combination value of variable loads which is 0.7 for snow load. $S_{Dk}$ and $S_{Srk}$ are the effects of characteristic dead load and characteristic snow load, respectively. $\gamma_L$ is the adjustment coefficient which should be considered if the design working life of a building is different from design reference period which is 50 years defined by GB 50009-2012. $\gamma_0$ is the importance factor of structures which is related to the category of structural safety. $R_k$ is the characteristic value of structural resistance, and $\gamma_R$ is partial factor for resistance which varies with the type and material of structural components and a value of 1.1 is used herein.

The characteristic value of snow load on a building roof is defined by the product of characteristic ground snow load and ground-to-roof conversion factor as follows,

$$S_{rk} = \mu_k S_{gk} \quad (3)$$

in which $\mu_k$ is the characteristic value of ground-to-roof conversion factor, which is also termed as the shape coefficient. $S_{gk}$ is the characteristic value of ground snow load $S_g$ and is specified as the 0.98 quantile of the probability distribution of $S_g$. Only uniform distribution of snow load is taken into consideration for the design of a flat roof in GB 5009-2012, and a conservative value 1.0 of $\mu_k$ is specified for the design of flat roofs. The values of $S_{gk}$ for different locations are tabulated in GB50009-2012.

Since the characteristic value of dead load depends on various factors, a ratio of characteristic snow-to-dead load, which is assumed at the range of [0.25, 2.0] in this paper, is used to determine it.

### 3 Assessment of structural reliability

Structural reliability index depends on the variability of structural resistance and loads, and corresponds to the

design working life of a building which is 50 years specified by *Unified standard for reliability design of engineering structures* (GB 50153-2008) for common buildings [5]. According to the severity of possible consequences caused by structure failure, buildings are classified into three safety categories in GB 50153-2008. The safety category of the studied roof structures in this paper is assumed to be II, therefore the importance factor $\gamma_0$ is 1.0 and the target reliability index is 3.2 which means the probability of structural failure in 50 years is not more than $6.9 \times 10^{-4}$ if the failure of structural members is ductile. The reliability indices of steel roof members in the three cities which include Beijing, Nanjing and Hohhot are assessed through Monte-Carlo simulation in this paper, and the influences on reliability assessment of different probabilistic models for ground snow loads is also investigated.

#### 3.1 Limit state function

Reliability of roof structures can be defined by the limit state function as follows:

$$g(R, D, S) = R - (C_D D + C_S S) \quad (4)$$

where $C_D$ and $C_S$ are the conversion coefficients of load to effect for dead load and snow load, respectively. Reliability assessment would not be influenced by assumption $C_D = C_S = 1$ when snow load and dead load are all uniformly distributed on the roof surface, which means that the spatial variability of dead and snow load are ignored. R, D, and S are dealt as random variables, and Monte-Carlo simulation is adopted to calculate structural reliability indices. Probability of structural failure equals to the probability of $g(R, D, S) < 0$.

#### 3.2 Structural resistance

The suggestion of Ellingwood (1982) is taken herein that resistance of a steel roof member can be defined by a lognormal distribution, with a mean 1.07 times its characteristic value and coefficient of variation 0.15 [6]. Characteristic value of resistance is determined by formula (2) and (3) according to Chinese load code GB 50009-2012.

#### 3.3 Dead load

Dead load of a roof structure mainly refers to its self-weight, which hardly changes with working time in general. Thereby, dead load is also treated as a random variable that does not vary with time. Normal distribution was suggested by Ellingwood (1982) to describe it [6], the mean value is 1.05 times its characteristic value, and coefficient of variation is 0.10.

#### 3.4 Snow load

Snow load on a roof is generally defined by the product of the ground-to-roof conversion factor $\mu$ and the ground snow load $S_g$. Conversion factor $\mu$ can be
modelled as a lognormal random variable referring to O'Rourke (1983) [7]. This and O'Rourke (2015) statistically analysed the database containing simultaneous measurements of roof and ground snow loads in southern Norway [8]. They obtained the shape coefficient μ for low-slope roofs (0-5°) with mean 0.73 and standard deviation 0.12. It is assumed herein that the ground-to-roof conversion factor μ for flat roofs in the three cities of China has the same parameters.

JCSS proposes that Ferry Borges-Castanheta process can be used to describe variable loads. Thus, the maximum value of snow load in the design working life T of a building is given by,

$$F_T(x) = [F(x)]^T \quad (5)$$

where F(x) is the cumulative distribution function of the annual maximum ground snow load, F_T(x) is the cumulative distribution function of the maximum ground snow load in the design working life of a building.

GB 50009-2012 proposes that Gumbel distribution should be used to characterize the annual maximum ground snow depth, and method of moments (MOM) should be used to estimate the parameters of Gumbel distribution. Nonetheless, Gumbel distribution, generalized extreme value (GEV) distribution, and lognormal distribution are considered for the simulated annual maximum ground snow loads. Method of moments (MOM), maximum-likelihood estimation (MLE), and least square method (LSM) are employed to estimate the parameters of Gumbel distribution. MOM and MLE are applied to lognormal distribution. MLE is adopted for GEV distribution. The Kolmogorov-Smirnov (K-S) test is adopted to examine the applicability of above mentioned probabilistic models, and Akaike information criterion (AIC) defined by formula (6) is used to measure the relative quality of above probabilistic models for annual maximum ground snow loads.

$$\text{AIC} = 2k + n \ln(\text{RSS}) \quad (6)$$

where n denotes the sample size of simulated snow loads and k denotes the number of parameters of distribution. RSS is the residual sum of squares. The preferred model is the one with the minimum AIC value.

The results of Kolmogorov-Smirnov test indicate that all of the above mentioned distributions and corresponding parameter estimation methods are applicable to the annual maximum ground snow loads. The relative quality of different models according to AIC are compared in Tab. 2. The ‘1st’ means the minimum AIC value, and the ‘6th’ implies the worst model.

<table><caption>Tab.2 AIC of different models</caption><thead><tr><th>Models</th><th>Beijing</th><th>Nanjing</th><th>Hohhot</th></tr></thead><tbody><tr><td>Gumbel(MOM)</td><td>5th</td><td>4th</td><td>2nd</td></tr><tr><td>Gumbel(MLE)</td><td>2nd</td><td>6th</td><td>5th</td></tr><tr><td>Gumbel(LSM)</td><td>1st</td><td>5th</td><td>1st</td></tr><tr><td>Lognormal(MOM)</td><td>6th</td><td>3rd</td><td>4th</td></tr><tr><td>Lognormal(MLE)</td><td>4th</td><td>2nd</td><td>3rd</td></tr><tr><td>GEV(MLE)</td><td>3rd</td><td>1st</td><td>6th</td></tr></tbody></table>

Furthermore, Q-Q plots are used to graphically compare the simulated data sets to the underlying theoretical models, as showed in Fig.2 (a)~1(c). The x

coordinate of a point on the plot corresponds to one of the quantiles of the simulated snow load data set, and the y coordinate of the point is the same quantile of an underlying theoretical distribution. If a distribution characterizes the data set well, the points in the Q-Q plot will approximately lie on the line y=x. It can be found that the model with the minimum AIC value fit the overall data best on the Q-Q plots of the three cities, but many outliers appear at low quantiles and high quantiles. Notably, the goodness of fit to high quantiles is more important for reliability assessment, compared to low quantiles, because that the maximum value of ground snow load in the design working life which is usually 50 years is directly related to structural failure. However, the high quantiles are too scattered to fit well. In order to clarify the influence of the models for annual maximum ground snow loads on structural reliability assessment, all the underlying distributions for annual maximum ground snow loads mentioned above are considered.

<figure><img src="image_170.png" /><figcaption>(a) Beijing</figcaption></figure>

<figure><img src="image_171.png" /><figcaption>(b) Nanjing</figcaption></figure>

<figure><img src="image_172.png" /><figcaption>(c) Hohhot</figcaption></figure>

<caption>Fig.2 Q-Q plots for comparison of different probability models to annual maximum ground snow loads</caption>
Basing on the distribution of annual maximum values, Monte-Carlo method is used to simulate the maximum ground snow loads in the design working life of a building, which is 50 years in general. A procedure of generating the maximum ground snow loads in 50 years on the basis of the distributions of annual maximum ground snow loads is explained as follows:

(1) Firstly, 1,000,000 sets of annual maximum ground snow load samples are generated by Monte-Carlo simulation, each set contains 50 data.

(2) Then, the maximum values in every set of annual maximum ground snow loads are picked out to act as the maximum ground snow loads in 50 years. Thus, the sample of maximum ground snow loads in the design working life of a building is assembled.

# 4 Results and discussion

As showed in Tab. 1 and discussed in section 3.4, the annual maximum ground snow loads in Beijing, Nanjing, and Hohhot have different variability. K-S test indicate that all of the 6 probabilistic models considered in this paper can be used for the annual maximum ground snow loads of the three cities, whereas different distributions and parameter estimation methods lead to great difference in estimating high quantiles. In order to explore the influences on structural reliability assessment of this difference, all of the 6 models are adopted as the basis of the distribution of the maximum ground snow loads in 50 years, and Monte-Carlo simulation is employed to calculate structural reliability indices. The results are displayed as Fig. 3.

<figure><img src="image_173.png" /><figcaption>(a) Beijing</figcaption></figure>

<figure><img src="image_174.png" /><figcaption>(b) Nanjing</figcaption></figure>

<figure><img src="image_175.png" /><figcaption>**Fig. 3.** Reliability indices for steel roof members designed by Chinese codes in different cities.</figcaption></figure>

Fig.3 (a) shows that different distributions and parameter estimation methods lead to great difference in results of structural reliability assessment for Beijing, which is consistent with the influence on high quantile estimation. Combining Fig.3 (a) with Fig.2 (a), one can find that the model which is a combination of Gumbel distribution and MLE gets the minimum estimation of high quantiles in Q-Q plots of Beijing, but lead to the maximum estimation of structural reliability indices. Conversely, the model that is composed of lognormal distribution and MLE gets the maximum estimation of high quantiles in Q-Q plots of Beijing, but lead to the minimum structural reliability indices. It also can be found from Fig.2 (a) that the model consisted by GEV distribution and MLE fit well to the high quantiles of Beijing's annual maximum ground snow loads, therefor the reliability assessment results basing on it is the most believable and other results are overestimated or underestimated.

From Fig.3 (b), one can find the similar results of Nanjing compared to Beijing that if a model underestimate the high quantiles, it will lead to an overestimated structural reliability assessment. Fig.2 (b) shows that the six models all underestimate the high quantiles of Nanjing, therefor the reliability indices of steel roof members in Nanjing are overestimated.
Fig.3 (c) shows that the reliability indices of steel roof members in Hohhot approximately lies at the range [1.7, 2.3], which are influenced less by the models for annual maximum ground snow loads compared to Beijing and Nanjing. The reason is that the 6 models have almost the same goodness of fit to the high quantiles of Hohhot.

Comparing the results of the three cities, one can find that the difference of high quantile estimation by the six models tends to become more and more obvious as the COV increases, and the same trend of the influence on reliability assessment can be found.

Furthermore, it can be found that the influence on structural reliability assessment of probabilistic models for snow loads become larger slightly with the increase of the ratio of characteristic snow-to-dead load. The reason for this is that the total load effect will have greater
variability if the ratio is larger on the condition that the variability of snow load is greater than that of dead load.

## 5 Conclusions

GB 50009-2012 proposes Gumbel distribution for annual maximum ground snow loads, and method of moments (MOM) for estimating parameters. Nonetheless, six different models which are the combination of three distributions and three parameter estimation methods are considered for annual maximum ground snow loads in this paper to investigate structural reliability indices of steel roof members in three representative cities of China. The results indicate that an appropriate probabilistic model for annual maximum ground snow loads is very important to structural reliability assessment. The main conclusions of this paper are the following:

(1) If a model overestimate the high quantiles of annual maximum ground snow loads, it will lead to conservative assessment of structural reliability. Whereas structural reliability will be overestimated if the model used for annual maximum ground snow loads underestimate the high quantiles.

(2) The larger the COV of annual maximum ground snow loads is, the more difficult it is to fit the high quantiles well, and more uncertainties in structural reliability assessment will be found.

(3) The influence on structural reliability assessment of probabilistic models for snow loads tends to be a little greater with the increase of the ratio of characteristic snow-to-dead load.

## References

1. T. Takahashi, and B. R. Ellingwood, Reliability-based assessment of roofs in Japan subjected to extreme snows: incorporation of site-specific data, Eng. Struct., 27:89-95 (2005)

2. D. L. Kozak, and A. B. Liel, Reliability of steel roof structures under snow loads, Struct. Saf., 54: 46-56 (2015)

3. X. Y. Zhou, Y. Q. Zhang, M. Gu, et al, Simulation method of sliding snow load on roofs and its application in some representative regions of China, Nat. Hazards, 67: 295-320 (2013)

4. Ministry of Construction of the People's Republic of China, General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China, Load code for the design of building structures (GB 50009-2012), China Building Industry Press, Beijing (2012)

5. Ministry of Construction of the People's Republic of China, General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China, Unified standard for reliability design of engineering structures (GB 50153-2008), China Building Industry Press, Beijing (2008)

6. B. Ellingwood, T. V. Galambos, Probability-based criteria for structural design, Struct. Saf., 1(1): 15-26 (1982)

7. M. O'Rourke, U. Stiefel, Roof snow loads for structural design, J. Struct. Eng., 109(7): 1527-1537 (1983)

8. T. K. Thiis, M. O'Rourke, Model for Snow Loading on Gable Roofs, J. Struct. Eng., 04015051 (2015)
Session 3

# Avalanche / Snow physics