# The importance of field measurements and observation systems in snow engineering: examples in the fields of avalanche dynamics and drifting snow

Florence Naaim-Bouvet

Univ. Grenoble Alpes, IRSTEA, UR ETNA, Saint Martin d'Hères, France

**Abstract.** At a time when public finances are low and the cost of setting up, operating, maintaining and managing observation systems high, it is essential to understand the importance of field data in snow engineering. This paper will explain and demonstrate this importance using practical examples in the fields of avalanche dynamics and drifting snow. For snow avalanche dynamics, a huge effort has been made to carry out and manage field measurements and observations and to disseminate data at the European or even global level. This has led to significant progress in related snow engineering applications. For drifting snow, however, much still needs to be done. Ongoing work and propositions aimed at enhancing the data base and improving the management of field observations will be presented.

## 1 Introduction

Observation systems capable of providing data for snow engineering purposes are costly to set up, operate, maintain and manage. At a time when public funding is increasingly hard to come by, it is important to understand the essential nature of field data in snow engineering. This paper presents a number of examples in the field of snow avalanche dynamics and drifting snow that clearly illustrate this point. They are based on my professional experience acquired in these fields as an engineer and researcher.

In snow science and more generally in geophysics, observation lies at the heart of research (Figure 1). The first process in the scientific method [1] involves the observation of a phenomenon, event or “problem.” Observation leads to a question that needs to be answered to satisfy human curiosity or solve a given problem. To answer the question, a hypothesis is formed. Once the hypothesis has been established, an experiment must be designed to prove or disprove the hypothesis.

In Earth sciences, including snow science, application of the scientific method can sometimes be complicated by the fact that the Earth is not a controlled laboratory. Many variables are involved. Scientists must design controlled and repeatable experiments to develop and validate theories. Note that there is a permanent exchange and mutual enrichment between theory, observation, numerical experiments and laboratory experiments.

Working at very different spatial scales represents an additional complexity. For example, for drifting snow, the scale of interest extends from that of the particle (if we focus on the splash function) to the continent (if we are interested in surface mass balance in Antarctica).

Snow engineers have various tools at their disposal to solve problems and design various works. These tools are mainly based on research results. As in snow science, observation is important and provides the necessary background (Figure 2). The choice of final design values requires statistical methods based on observations. Moreover, the analysis of experts relies on personal observations compiled in their minds over the years to form "living" data bases. Physical and numerical simulations are often used to solve a given problem but

field experiments first provide the benchmarks for the development and calibration of such models.

Snow research contributes to snow engineering and snow engineering in turn provides feedback, questions to investigate and practical problems to solve.

<figure><img src="image_176.png" /><figcaption>Figure 1. The role of observation in snow science</figcaption></figure>

<figure><img src="image_177.png" /><figcaption>Figure 2. The role of observation in snow engineering</figcaption></figure>

In an attempt to be more practical than theoretical, the importance of field measurements and observation systems in snow engineering will be illustrated in this paper with specific examples. For snow avalanche dynamics, a huge effort has been made to carry out and manage field measurements and observations and to disseminate data at European and even global levels.

<figure><img src="image_178.png" /></figure>
Significant results have been transferred to snow engineering. For drifting snow, much work remains to be done and many questions have not yet been answered.

The paper is organized as follows: Section 2 discusses the importance of field measurements from avalanche test sites, presents data from avalanche surveys and deals with the combined use of avalanche surveys and numerical models. Section 3 focuses on the lack of data for drifting snow and briefly introduces ongoing attempts to improve and better manage observations in this field.

## 2 Avalanche dynamics

Snow avalanches represent a serious problem in our society. In the United States and Canada, most avalanche accidents involve recreational users, but in other part of the world houses and businesses are also at risk (01-02/1992: 261 deaths in Turkey – 05/1998: 250 deaths in India – 01/1993: 239 deaths in Russia – 03/1997: 100 deaths in Afghanistan – 01/1954: 200 deaths in Austria – 1950/1951: 98 deaths in Switzerland – 02/1970: 39 deaths in France – 02/1999: 17 deaths in Switzerland and 12 deaths in France).

As we can see, there are many areas in Europe that have a long and tragic history of avalanche accidents. This is probably the reason that research on this topic is very active in Europe, as can be seen from an analysis of search results in the Scopus data base presented for different countries (Figure 3).

<figure><img src="image_179.png" /><figcaption>**Figure 3.** Scopus data base search results presented for different countries (search keywords: “Avalanche” + “Flow” / date: May 2016)</figcaption></figure>

### 2.1 Field measurements on avalanche test sites

All European research centres dealing with avalanche dynamics are aware of the importance of field measurements and perform full-scale avalanche experiments even though they are costly, time-consuming and dangerous. Such experiments provide data that is essential to the understanding and modelling of avalanche motion, even if it can take several years to obtain reliable results. Many different configurations of avalanche types and paths need to be studied. A number of European countries have jointly planned and conducted such experiments through three successive European FP projects (SAME, CADZIE and SATSIE) [2]) and two bilateral INTERREG ALCOTRA projects (DYNAVAL

and MAP³) [3]. A dozen full-scale test sites have been involved, though at present only five are still operative (Table 1).

Rather than describe the sites and their instrumentation and configuration in detail, we will focus on their complementarities in terms of configuration and range of pressure, velocity, path length and drop height, as can be seen clearly in Table 1. The smallest site is in Italy, with a path length of 400 m and a maximal pressure of 0.5 tonnes/m². The biggest site is in France, with a path length of 7000 m and a maximal pressure of 100 tonnes/m². Note that all these sites are open to collaborations and can accommodate sensors provided by other sites or teams.

The Taconnaz avalanche path (Chamonix, France) (Figure 5) [4] is extremely long. An international road and several buildings are located in the area at risk, leading to the design of an avalanche protection system. Avalanches mainly consist of dry snow with ice blocks breaking from the glacier and cannot be artificially released.

The Vallée de la Sionne (Switzerland) avalanche path (Figure 4) [5] ends with an inverse slope on which a bunker houses staff during the experiments. The site is probably the best equipped with a unique collection of field-hardened sensors to measure avalanche velocity, impact pressure, flow depth and density.

The Ryggfonn (Norway) (Figure 6) [6] runout ends at a catching dam which was built to test the effectiveness of such a mitigation measure.

As opposed to others full-scale avalanche test sites in Europe, the runout at Lautaret (Figure 7) [7] is a natural open slope. Flow height, vertical profiles of pressure and velocity are measured along a mast which was destroyed in 1995 and subsequently rebuilt. Note that measurement equipment has been destroyed at least once on all the experimental sites.

<figure><img src="image_180.png" /><figcaption>**Figure 4.** Vallée de la Sionne avalanche test site (photo courtesy of SLF)</figcaption></figure>
<figure><img src="image_181.png" /><figcaption>**Figure 5.** Taconnaz avalanche test site (inset: ice blocks in snow avalanche deposit)</figcaption></figure>

<figure><img src="image_182.png" /><figcaption>**Figure 6.** Ryggfonn (avalanche test site (photo courtesy of NGI))</figcaption></figure>

<figure><img src="image_183.png" /><figcaption>**Figure 7.** Lautaret avalanche test site</figcaption></figure>

Seehore peak (Figure 8) [8] in Italy is the latest test site to have been set up and is used to study artificially released small and medium avalanches, while the majority of the existing sites study large avalanches.

<figure><img src="image_184.png" /><figcaption>**Figure 8.** Seehore peak avalanche test site (photo courtesy of E. Bovet)</figcaption></figure>

Much has been learned over the last 10-15 years from this well-organized network of avalanche test sites. Recorded data helps us understand the physical processes of a flowing avalanche or an avalanche impacting on catching dams or others obstructions. Such data is also used to calibrate avalanche-dynamics models. In this paper, we will focus on two main results now used in snow engineering.

### 2.1.1 Erosion and deposition processes

The analytical and numerical models used in practical calculations of snow avalanche dynamics contain a questionable simplification: they assume that avalanche mass is constant meaning that no snow cover entrainment or mass deposition takes place along the avalanche path.

In order to establish the mass balance of an avalanche, he distribution of the snow before and after avalanche release (Figure 9) has been measured first by photogrammetry and then by laser scanning [9]. These measurements have shown that the avalanche mass can be 6-7 times higher than the released mass [10]. This has led to the development of advanced erosion and deposition models used in numerical models. The validation of such advanced models requires in turn more detailed measurements of the processes within the avalanche and therefore more sophisticated field measurements.

<figure><img src="image_185.png" /><figcaption>**Figure 9.** Map of the elevation difference as measured by laser scanning for the avalanche released on the Col du Lautaret path on 13 February 2013. [7]</figcaption></figure>
<table><caption>Table 1. Key features of avalanche test sites (updated from [2])</caption><thead><tr><th></th><th>Seehore<br>Italy<br>2009<br>Small size</th><th>Lautaret<br>France<br>1973<br>Small to medium size</th><th>Ryggfonn<br>Norway<br>1980<br>Medium size</th><th>Vallée de La Sionne<br>Switzerland<br>1999<br>Large size</th><th>Taconnaz<br>France<br>2010<br>Large size</th></tr></thead><tbody><tr><th>Release</th><td>Artificial</td><td>Artificial</td><td>Artificial</td><td>Artificial/natural</td><td>Natural</td></tr><tr><th>Path length</th><td>400-500 m</td><td>800 m</td><td>2100 m</td><td>2700 m</td><td>Up to 7000 m</td></tr><tr><th>Drop height</th><td></td><td>450 m</td><td>975 m</td><td>1300 m</td><td>3000 m</td></tr><tr><th>Mean slope</th><td>38°</td><td>36°</td><td>31°</td><td>29°</td><td>25°</td></tr><tr><th>Deposit volume</th><td>200-400 m³</td><td>500-10,000 m³</td><td>10,000 m³<br>- 100,000 m³</td><td>300,000 m³</td><td>1,600,000 m³</td></tr><tr><th>Pressure</th><td>0.5 t/m²</td><td>3-5 t/m²</td><td>60 t/m²</td><td>Up to 100 t/m²</td><td>Up to 100 t/m²</td></tr><tr><th>Velocity</th><td>Up to 30-40 m/s</td><td>Up to 30-40 m/s</td><td>Up to 60 m/s</td><td>Up to 50-70 m/s</td><td>Up to 70 m/s</td></tr></tbody></table>

### 2.1.2 Impact pressure

The impact pressure $P$ can be expressed as a function of the dynamic pressure where $C_D$ is the effective drag coefficient (eq. 1).

$$P = \frac{C_D \rho U^2}{2} \quad [1]$$

where $\rho$ the flow density and $U$ the flow velocity upstream of the obstacle.

$C_D$ is generally taken as 1. It has been shown that the drag coefficient increases as the Reynolds number decreases (Figure 10) so that high pressures can still be exerted at low velocities in the runout zone, as observed in MontRoc (Figure 11) near Chamonix during the 1999 avalanche disaster in France. Avalanche hazard maps are used as regulation for land use planning. It shows the maximum between the centennial and maximal observed run-outs. The hazard levels are colour-coded red, blue and white according to the impact pressure. Using a constant value of 1 for the drag coefficient leads to an underestimation of the avalanche hazard in the locations where the velocities are low such as run-out zones.

<figure><img src="image_186.png" /><figcaption>Figure 10. $C_D$ as function of the Reynolds number [11]</figcaption></figure>
<figure><img src="image_187.png" /><figcaption>**Figure 11.** High pressure exerted at low velocities in the runout zone of Péclerey avalanche that hit MontRoc village (9 February 1999) – F. Rapin/IRSTEA.</figcaption></figure>

## 2.2 Observatories: permanent avalanche surveys

Observatories generally measure fewer variables but offer greater data availability over time compared to specific field measurement campaigns that are more often dedicated to particular research topics and practical questions.

In France, an extensive avalanche data base is available. The Enquête Permanente sur les Avalanches (E.P.A.) is a descriptive report including quantitative (runout altitudes, deposit volumes,...) and qualitative (flow regime, release cause,...) information for selected avalanche paths. This chronicle is characterized by its long duration. Furthermore, it was set up by a forest engineer and not by a scientist. The first avalanches were recorded in the Savoie region in 1899. 3936 paths are monitored each winter, involving 465 towns and 260 observers. In all, 85,000 events are available over the last century (www.avalanches.fr).

The data series in this unique data base are unusually long, even if some of them are incomplete and suffer from quality problems. This data base was initially devoted to the evaluation of forest damage but is now highly valuable for other applications. For instance, it provides an accurate view of the spatio-temporal fluctuations of avalanche activity over the last century in the context of climate change.

At the scale of the entire French Alps, a major turning point can be noted around 1978 [13]: between 1946 and 1978 we observe an increase of 0.1 avalanche occurrences per winter and per path and a roughly 55 m lower runout altitude. And between 1978 and 2010, the reverse occurs: we notice a decrease of 0.1 avalanche occurrences per winter and per path and the runout altitude rose 55 m.

The same trend is observed if we consider the runout altitude corresponding to a return period of 10 years. Considering the clearly non-stationary nature of this phenomenon, the issue that must be addressed is: How can we evaluate long term avalanche hazards?

There has also been a general decrease of around 12% in the proportion of powder snow avalanches since 1973.

This decrease appears to be consistent with increasing winter temperatures over the same period.

## 2.3 Combined use of permanent avalanche surveys and numerical models

### 2.3.1 Local predetermination using statistical-dynamical approaches

Even if we consider stationary conditions, determination of a 100-year event in the framework of protection design is not so obvious.

Legal regulations for land use planning are based on return periods that vary depending on the country. For example, it is a 100-year return period in France. The 100-year volume and runout distance can be directly calculated from the EPA data base. However, protection dam design requires estimates of impact pressure. The example of the Taconnaz avalanche provides greater insight into how the problem has been solved by combining an avalanche-dynamics model and the EPA data base.

The Taconnaz avalanche path (Figure 5) is located in the Chamonix valley close to Mont Blanc. It is the largest European avalanche path with roads and buildings located within the zone of risk. The Taconnaz path is 7 km long with a vertical drop of 2 km and a mean slope of 25°. A total of 75 events have been recorded in the data base. Each recorded event was back-calculated using a numerical model of dense flow to define parameters (the volume and the Froude number, F defined as the square root of the ratio between the kinetic and potential energies corresponding to the maximum discharge) at the entry of the runout zone upstream of the defence structures. Then, a statistical analysis of these parameters allowed characterization of 100-year return period events. The 100-year values were 3.58 for the Froude number and 1.6 Mm³ for the volume. Finally, physical (Figure 12) and numerical (Figure 13) models of dense avalanches interacting with defence structures were combined in order to design the most effective passive structure to contain the reference scenarios.

<figure><img src="image_188.png" /><figcaption>**Figure 12.** Small-scale physical model of the runout zone fed with granular flows from a channel of adjustable slope and equipped with a reservoir at the top to release granular material</figcaption></figure>
<figure><img src="image_189.png" /><figcaption>**Figure 13.** Final deposits determined by numerical simulation</figcaption></figure>

### 2.3.2 Correlation between friction coefficients and the physical properties of snow

Even rough data can be used to gain an in-depth understanding of avalanche flow dynamics.

Naaim et al. [15] have combined the EPA data base with avalanche-dynamics models and snowpack models to link avalanche friction coefficients and the physical properties of snow.

A total of 735 historical avalanches released on 26 paths and recorded in the EPA chronicle since 1958 have been back-analyzed using a depth-averaged model with a high-resolution digital model acquired by laser scanning. The initial snow depth and snow cohesion were computed using the detailed CROCUS snowpack model fed with data from a SAFRAN meteorological analysis. More than 940,000 simulations were performed.

The results show that the static friction coefficient increases quasi-linearly with temperature (Figure 14), meaning that the runout distance should decrease with temperature. This is consistent with previous results which suggested that the decrease in runout distance was due to atmospheric warming that modified the average characteristics of snow. Focusing on liquid water content (Figure 15), it can be shown that the static friction coefficient decreases strongly for high values, typically those higher than 35 kg m⁻³

<figure><img src="image_190.png" /><figcaption>**Figure 14.** Back-analyzed static friction coefficients $\mu_0$ as function of average snow temperature T [15]</figcaption></figure>

<figure><img src="image_191.png" /><figcaption>**Figure 15.** Back-analyzed static friction coefficients $\mu_0$ as versus liquid water content [15]</figcaption></figure>

This section has provided a brief overview of practical lessons learned from avalanche dynamics field measurements and observations, which are not dedicated only to the validation of numerical models, as some may think. A huge effort has clearly been made to carry out and manage field measurements and observations and to disseminate data. The progress achieved in the field of avalanche dynamics is the result of many years of data acquisition and human and financial investments. In the field of drifting snow, much remains to be done, as will be discussed in the next section. Neither numerical and physical models nor expert knowledge can answer all practical questions, as shown by numerous examples. For this reason, efforts must be continued to set up and manage observations in the field.
# 3 Drifting snow

## 3.1 Estimation of drifting snow fluxes

The estimation of drifting snow fluxes remains an open question. Consider for example the case of Gouter refuge [16].

The new Gouter refuge was constructed in 2012 at a height of 3835 m on the normal route to Mont-Blanc, 200 meters from the old refuge. Drifting snow was a problem at the old refuge (Figure 16). For the new refuge, the architects wanted to design a building with main and backup exits free of snow and a drinking water supply fed by melting snow (Figure 17). Small scale experiments were carried out in a wind tunnel (Figure 18) and led to several recommendations for the building.

- Shape: aerodynamic (ovoid)

- Orientation: major axis of the ellipse parallel to the prevailing wind direction

- Layout: entrance in the front of the building and the melting container in the back

- Overhang: as much as possible.

<figure><img src="image_192.png" /><figcaption>**Figure 16.** Snowdrift around the former buildings. The red arrow indicates the prevailing wind direction</figcaption></figure>

<figure><img src="image_193.png" /><figcaption>**Figure 17.** Melting container in the back of the building</figcaption></figure>

<figure><img src="image_194.png" /><figcaption>**Figure 18.** Small-scale model in the wind tunnel (H. Raguet)</figcaption></figure>

However one question remains unanswered. It concerns the amount of snow that can be expected to cover the melter for a given duration. This is an important question given that the melter container must supply water for 120 people. Unfortunately, it cannot be answered for several reasons. First of all, there is no consensus on similitude requirements for storm duration [17]. However, the main problem concerns data. In fact there are not enough data to quantify a reference snow storm for the project (annual snow storm, for example). Furthermore even if wind data were available, drifting snow quantities could not be estimated accurately. Empirical formulae available in the literature produce divergent results [18].

## 3.2 Projected trends in the context of climate change

The formation of snowdrifts on roads causes disturbances to communication and transportation systems.

This was the case during the winter of 2005 when the A75 motorway linking Belgium and Spain was closed for 6 days due to the formation of snowdrifts (Figure 19) at La Fageolle pass, one of the highest motorway passes in Europe [19].

<figure><img src="image_195.png" /><figcaption>**Figure 19.** A75 closed to traffic: snow-removal operations</figcaption></figure>

In France, drifting snow is generally controlled by artificial snow fences. However, forested strips are being
used increasingly because they are better adapted to multidirectional wind. Unfortunately, such a solution takes time to become effective, as illustrated by the photos in Figure 20 and Figure 21, taken 26 years apart.

<figure><img src="image_196.png" /><figcaption>**Figure 20.** Forested strips planted in 1991 to control snowdrifts at La Fageolle pass</figcaption></figure>

<figure><img src="image_197.png" /><figcaption>**Figure 21.** The same forested strips in 2007</figcaption></figure>

Decision makers voiced strong doubts when it was suggested to extend the existing forested strip to strengthen snowdrift control. They wanted assurance that drifting snow storms would still exist at this location 30 years from now, when the proposed solution would become effective.

Unfortunately, we lack the necessary data to assess the effects of climate change on blowing snow storms.

## 3.3 Field measurements and observations for drifting snow

Even if drifting snow represents a serious problem, it is less important in our society than avalanches. Research on drifting snow is not mainly conducted in the countries of the Alps (Figure 22) and the research network is less organized. Even if there is a need, the strategy differs from that adopted for avalanches.

<figure><img src="image_198.png" /><figcaption>**Figure 22.** Scopus data base search results presented for different countries (search keywords: "Drifting snow" + "Blowing snow" / date: May 2016)</figcaption></figure>

For France there are two experimental test sites dedicated to drifting snow within the framework of the CRYOBS-CLIM observatory: Lac Blanc Pass and Adélie Land.

CRYOBS-CLIM (the CRYosphere, an OBServatory for the CLIMate) is aimed at gathering monitoring results and long-term observations performed by French research institutes on mountain glaciers, polar ice-sheets, seasonal snow cover and mountain permafrost and to provide relevant tools and observational data to a wide community so that they can address major societal issues such as water resources, ice-related hazards, atmospheric processes, avalanches, sea level rise and global climate change. The long-term observation strategy is based on in-situ (automated and manual) and remote (photogrammetry, LiDAR, satellite) monitoring of snow- and ice-related variables. Observations will be disseminated to the entire scientific community through a data base (http://www.cryobsclim.osug.fr/) that will be released in 2017.

### 3.3.1 Lac Blanc Pass

Lac Blanc Pass (2720 m a.s.l.), in the French Alps, is a part of the CRYOBS-CLIM network. It has been operated by IRSTEA and Météo France since 1990 and is dedicated to the observation of a typical high altitude snowpack strongly affected by drifting and blowing snow.

The aim is to study the physical processes involved in a mountainous context, create a data base for numerical modelling and create a 25-year long climatological reference at 2720 m a.s.l.

Due to the surrounding topography, winds are 80% from the north-east or the south. Moreover snow transport is observed 10.5% of the time in winter [21]. Initially, the site was mainly equipped with conventional meteorological stations and a network of snow stakes (Figure 23). Then new sensors (Figure 24) and technologies (ultrasonic anemometers, drifting snow sensors including snow particle counters and laser scanners [22]) (Figure 25) were installed, making it possible to develop new knowledge dealing with the relationships between threshold velocities and
morphological features of snow grains and snow flux profiles including parameters such as fall velocity and Schmidt number, histograms of particle widths, aerodynamic roughness, sastrugi and gust factors.

The physical processes involved in drifting snow have thus been studied more deeply, making it possible to better take them into account in models. At the same time, the instrumentation has made it possible to estimate the output parameters (snow depth) with a higher spatial resolution.

Note that data from 28 Flowcapt acoustic snowdrift flux sensors, operated in mountainous environments (France, Spain, Italy, Switzerland, Chili), are available on line (http://www.iav-portal.com/).

FlowCapt acoustic sensors, designed to measure the fluxes of snow transported by wind, have been compared to optical snow particle counters, considered herein as the reference at the Lac Blanc Pass, where a bench test for the wind transport of snow was set up. The two existing generations of FlowCapt have been compared. Both seem to be good detectors for the transport of snow by wind, but refinements still need to be made to accurately estimate drifting snow fluxes [23]. The third generation is currently being tested.

<figure><img src="image_199.png" /><figcaption>**Figure 23.** Conventional meteorological station and a network of snow stakes</figcaption></figure>

<figure><img src="image_200.png" /><figcaption>**Figure 24.** Ultrasonic anemometer, wind anemometer profile, acoustic and optical snow transport sensors</figcaption></figure>

<figure><img src="image_201.png" /><figcaption>**Figure 25.** Terrestrial laser scanner and general view of the experimental site.</figcaption></figure>

### 3.3.2 Adélie Land

More recently a similar experimental test site has been set up in Antarctica (Figure 26). Measurements are less intensive due to logistical constraints [24]. The site is managed by LGGE (Laboratory of Glaciology and Environmental Geophysics, CNRS, Grenoble) with the logistical support of the French Polar Institute.

Net erosion of snow by the wind may contribute significantly to the surface mass balance of the Antarctica coastal zone. However, there are very few field observations to confirm this hypothesis and to evaluate numerical models developed for this purpose [25]. Adélie Land is characterized by strong and frequent katabatic winds which are the strongest in the world in terms of mean wind speed at sea level. They are often associated with wind transport of snow. A field campaign was therefore launched in 2009. The observation system consisting in three automatic weather stations was deployed with Flowcapt acoustic drifting snow sensors close to the coast. The stations were set up at different locations, ranging from 1 to 100 km inland. One of them has a 7-m mast with six levels of anemometers and thermohygrometers (Figure 27).

<figure><img src="image_202.png" /><figcaption>**Figure 26.** Work environment in Antarctica (photo courtesy of LGGE - Grenoble)</figcaption></figure>
<figure><img src="image_203.png" /><figcaption>**Figure 27.** Flowcapt sensor (foreground) with a 7-m mast with six levels of anemometers and thermohygrometers (photo courtesy of LGGE - Grenoble)</figcaption></figure>

Drifting snow is a complex phenomenon at different spatial scales and is representative of the mutual interaction between the atmosphere and snow cover. Accurate field measurements covering long periods of time are extremely difficult to acquire. As is the case with avalanches dynamics, addressing this problem therefore requires a combination of resources, teams, sites, specialties and expertise.

Such a measurement programme is underway at Lac Blanc Pass and Adélie Land, two French experimental sites managed by three different institutes (IRSTEA, Météo France and LGGE). It will make it possible to observe and model drifting snow and validate and refine new and existing models. Further progress could be achieved by sharing this information at a global level.

## Conclusions

The importance of field data in snow engineering has been demonstrated thorough the paper using specific examples dealing with avalanche dynamics and drifting snow.

It is true that operating and maintaining experimental sites and observatory systems requires tremendous effort.

Nevertheless, given the importance of such measurements in snow research and risk management, research institutes should be encouraged to secure funds and human resources to maintain the level of these data series.

The importance of sharing data to cover a large spatial and temporal scale is well recognized and we must enhance and take full advantage of the complementary nature of the different sites.

Data exchange is of course a sensitive topic. Nevertheless several initiatives to this end have already been taken at a global level concerning different aspects of the cryosphere including snow. The objectives include facilitating:

- the standardisation of measurements and data collection

- the analysis, archiving and publication of such data

- the interaction and collaboration between scientific and operational communities.

Some examples we can cite include GCW (Global Cryosphere Watch from World Meteorological Organization), IACS (International Association of Cryospheric Science) and GEWEX-INARCH (International Network for Alpine Research Catchment Hydrology). However, effective use of the different networks can be confusing. All aspects related to snow engineering may not be taken into account and the cryosphere includes components other than snow such as lake, river and sea ice, glaciers, ice caps and ice sheets, frozen ground and permafrost.

In conclusion, much work remains to be done and observation systems and field measurements relevant to snow engineering still need to be organized at a global level.

## References

1. C.V. McLelland, The Nature of Science and the Scientific Method, The Geological Society of America, (2006).

2. Barbolini, M., Issler, D., Avalanche Test Sites and Research Equipment in Europe: An Updated Overview. SATSIE Avalanche Studies and Model Validation in Europe, (2006).

3. Manuel pour l'étude de l'intéraction de l'écoulement d'une avalanche avec un obstacle, projet n°48 DYNAVAL, (2012).

4. H. Bellot, F. Naaim-Bouvet, M. Naaim, P. Caccamo, T. Faug, F. Ousset, Taconnaz avalanche path: pressure and velocity measurements on breaking mounds, *Proceedings of International Snow Science Workshop Grenoble – Chamonix Mont-Blanc*, 1378-1382, (2013).

5. B. Sovilla, J. McElwaine, W. Steinkogler, M. Hiller, F. Dufour, E. Surinach, C Perre Guilen, J-T Fisher, E. Thibert, D. Baroudi, The full-scale avalanche dynamics test site Vallée de la Sionne, *Proceedings of International Snow Science Workshop Grenoble – Chamonix Mont-Blanc*, 1372-1377, (2013).

6. P. Gauer, K. Kristensen, Four decades of observations from NGI's full-scale avalanche test site Ryggfonn – Summary of experimental results, Cold Regions Science and Technology, **125**, 161-176 (2016).

7. E. Thibert, H. Bellot, X. Ravanat, F. Ousset, G. Pulfer, M. Naaim, P. Hagenmuller, F. Naaim-Bouvet, T. Faug, K. Nishimura, Y. Ito, D. Baroudi, A. Prokop, P. Schon, A. Soruco, C. Vincent, A. Limam, R. Hero, The full-scale avalanche test-site at Lautaret Pass (French Alps), Cold Regions Science and Technology, **115**, 30-41, (2015).

8. M. Maggioni, M. Freppaz, E. Ceaglio, D. Godone, D. Viglietti, E. Zanini, M. Barbero, F. Barpi, M. Borri Brunetto, E. Bovet, B. Chiaia, V. De Biagi, B. Frigo, O. Pallara, A new experimental snow avalanche test site at Seehore peak in Aosta Valley (NW Italien Alps) – part I: Conception and logistics, Cold
Regions Science and technology, **85**, 175-182,
(2013).

9. A. Prokop, P. Schon, F. Singer, G. Pulfer, M. Naaim, E. Thibert, A. Soruco, Merging terrestrial laser scanning technology with photogrammetric and total station data for the determination of avalanche modeling parameters, Cold Regions Science and Technology, **110**, 223-230, (2015).

10. Sovilla, B., and Bartelt, P., 2002. Observations and modelling of snow avalanche entrainment, Natural Hazards and Earth System Sciences, **2**, 169–179, (2002)

11. M. Naaim, N. Eckert, G. Giraud, T. Faug, G. Chambon, F. Naaim-Bouvet, D. Richard, Impact du réchauffement climatique sur l’activité avalancheuse et multiplication des avalanches humides dans les Alpes Françaises, La Houille Blanche, (2016).

12 E. Bourova, E. Maldonado, J-B Leroy, R. Alouani, N. Eckert, M. Bonnefoy-Demongeot, M. Deschatres, A new web-based system to improve the monitoring of snow avalanche hazard in France, Nat. Hazards Earth Syst. Sci., **16**, 1205-1216, (2016).

13. N. Eckert, C.J Keylock, H. Castebrunet, A. Lavigne, M. Naaim, Temporal trends in avalanche activity in the French Alps and subregions : from occurrences and runout altitudes to unsteady, Journal of Glaciology, **59**(213), 93-114, (2013).

14. M. Naaim, T. Faug, F. Naaim-Bouvet, N. Eckert, Return period calculation and passive structure design at the Taconnaz avalanche path, France, Annals of Glaciology, **51**(54), 89-97, (2010).

15. M. Naaim, Y. Durand, N. Eckert, C. Chambon, Dense avalanche friction coefficients: influence of physical properties of snow, Journal of Glaciology, **59**(216), 771-782 (2013).

16. F. Naaim-Bouvet, J. Pietu, Building's design in a blowing snow context: case of mountain retreat « Gouter » (Mont-Blanc Massif), 12th congress Interpraevent, Grenoble, France, (2013).

17. F. Naaim-Bouvet, Comparison of requirements for modeling snowdrift in the case of outdoor and wind tunnel experiments, Surveys in Geophysics, **16**, Number 5-6, 711-727, (1995).

18. J.J Doorschot and M. Lehning, Equilibrium saltation: mass fluxes, aerodynamic entrainment, and dependence on grain properties, Boundary layer Meteorology, **104**, 111-130 (2002).

19. F. Naaim-Bouvet, S. Monier, V. Mauduit, M. Boulet, France's A75 Motorway: protection from snowdrifts at the Fageole summit, 14th International Winter Road Congress, Andorra, (2014).

20. F. Naaim-Bouvet, G. Guyomarc'h, H. Bellot, Y. Durand, M. Naaim, V. Vionnet, C. Genthon, K. Nishimura, Y. Ito and A. Prokop, Lac Blanc Pass: a natural wind-tunnel for studying drifting snow at 2700ma.s.l., International Snow Science Workshop Grenoble Chamonix-Mont-Blanc, France, October 7-11, 1332-1339, (2013).

21. V. Vionnet, G. Guyomarc'h, F. Naaim Bouvet, E. Martin, Y. Durand, C. Bel, H. Bellot, P. Puglièse, Occurrence of blowing snow events at an alpine site

over a 10-year period: observations and modeling, Adv. Water Res., **55**, 53-63, (2013).

22. P. Schön, A. Prokop, V. Vionnet, G. Guyomarc'h, F. Naaim-Bouvet, M. Heiser, Improving a terrain-based parameter for the assessment of snow depths with TLS data in the Col du lac Blanc area, Cold regions Science and Technology, **114**, 15-26, (2015).

23. A. Trouvilliez, F. Naaim-Bouvet, H. Bellot, C. Genthon, and H. Gallee, Evaluation of the Flowcapt acoustic sensor for the aeolian transport of snow, Journal of Atmospheric and Oceanic Technology, **32**(9), 1630–1641, (2015).

24. A. Trouvilliez, F. Naaim-Bouvet, C. Genthon, V. Favier, L. Piard, H. Bellot, C. Agosta, C. Palerme, C. Amory and H. Gallée, A novel experimental study of aeolian snow transport in Adelie Land (Antarctica), Cold regions Science and Technology, **108**, 125-138 (2014).

25. H. Gallee, A. Trouvilliez, C. Agosta, C. Genthon, V. Favier, F. Naaim-Bouvet, Transport of snow by the wind: a comparison between observations in Adélie Land, Antarctica, and simulations made with the Regional Climate Model MAR. Boundary-Layer Meteorology, **146**, Issue 1, 133-147 (2013).
# Effect of reporting rate on vulnerability with an example for snow avalanche risk to backcountry recreationists in Canada

Bruce Jamieson<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> and Alan S.T. Jones<fnref n="3" />

**Abstract.** The vulnerability component in risk assessments or risk-based design may be based on damage models, expert judgement or records of damage. When based on records of damage, we show the calculated vulnerability is sensitive to the reporting rate for events that did not cause the specified level of damage, and derive an adjustment to the vulnerability. Using recreationists caught in snow avalanches in Canada – for which non-fatal involvements are poorly reported – we apply the adjustment and discuss the resulting vulnerability.

## 1 Introduction

In assessments for snow and other hazards, risk includes three components: probability of the event (typically for each scenario of specified magnitude), exposure of elements of value, and vulnerability. This third component, vulnerability V, is the expected fraction of loss, $0 \le V \le 1$, given that the elements of value are impacted. Vulnerability can be based on damage models, records of damage, or expert judgement. When vulnerability V is based on records, the reporting rate for events that caused a specified level of damage, $N_D$, and events that did not, $N_N$, should be similar, otherwise the vulnerability will be biased towards the better recorded event.

<figure><img src="image_204.png" /><figcaption>**Figure 1.** Sensitivity graph for the effect of the reporting rate R on vulnerability.</figcaption></figure>

## 2 Adjusting vulnerability for the reporting rate

When based on records of damage, vulnerability V can be expressed

$$V = \frac{N_D}{N_{Impacted}} = \frac{N_D}{N_D + N_N} \quad (1)$$

While events that cause damage tend to be well recorded, in some situations events that do not cause damage are not as well recorded. If the reporting rate for events that do not cause the specified level of damage is R < 1, then the corrected vulnerability $V_R$ is

$$V_R = \frac{N_D}{N_D + N_N/R} \quad (2)$$

which assumes that all events that cause the specified level of damage, $N_D$, are recorded.

From Eq. 1,

$$N_N = N_D(1/V - 1) \quad (3)$$

Substituting Eq. 3 in Eq. 2 yields

$$V_R = \frac{V}{V+(1-V)/R} \quad (4)$$

Figure 1 shows the adjusted vulnerability $V_R$ for R = 0.1, 0.4 and 0.9. For most risk assessment applications, R ≥ 0.9 likely yields sufficient accuracy; however, lower values of the reporting rate could compromise the accuracy of vulnerability and hence risk calculations.

Equation 3 is not specific to risk due to snow hazards. It can be applied to any hazard for which vulnerability is based on records of damage. It may be relevant to risk to structures such as masts, snowpack support structures and buildings due to snow creep, snow glide or snow avalanches. However, we recognize that vulnerability may — alternatively — be based on damage models [1] or expert judgement. Also, many structures threatened by natural hazards are designed based on safety factors applied to expected loads [2] rather than on risk.

## 3 Example for the vulnerability of recreationists caught in potentially fatal snow avalanches

Most studies of vulnerability for people caught in snow avalanches assume that all non-fatal involvements are reported. However, in a recent Canadian study Jamieson and Jones [3] used the results of a survey to argue that only about 10 % of persons caught in, but not killed by, an avalanche during recreation in Canada are reported.

In Canada, avalanche size is reported based on a five level scale for destructive potential from 1 (relatively harmless) to 5 (largest known) [4] (Table 1). Half-sizes are often used. The size classes are sometimes prefixed by a D

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

<footnote n="1">¹ Snowline Associates Ltd., Calgary, Canada</footnote>

<footnote n="2">² Dept. of Civil Engineering, University of Calgary, Canada</footnote>

<footnote n="3">³ Dynamic Avalanche Consulting Ltd., Revelstoke, BC, Canada</footnote>

<footnote>a Corresponding author: bruce.jamieson@snowline.ca</footnote>
to distinguish the classes from other international classifications.

<table><caption>**Table 1.** Classes of avalanche size by destructive potential [4]</caption><thead><tr><th>Size<sup>a</sup></th><th>Destructive potential</th><th>Typical mass (t)</th><th>Typical path length (m)</th></tr></thead><tbody><tr><td>D1</td><td>Relatively harmless to people.</td><td>&lt;10</td><td>10</td></tr><tr><td>D2</td><td>Could bury, injure or kill a person.</td><td>10<sup>2</sup></td><td>100</td></tr><tr><td>D3</td><td>Could bury a car, destroy a small building, or break a few trees.</td><td>10<sup>3</sup></td><td>1000</td></tr><tr><td>D4</td><td>Could destroy a railway car, large truck, several buildings, or a forest with an area up to 4 ha.</td><td>10<sup>4</sup></td><td>2000</td></tr><tr><td>D5</td><td>Largest snow avalanches known; could destroy a village or a forest of 40 ha.</td><td>10<sup>5</sup></td><td>3000</td></tr></tbody><tfoot><tr><td colspan="4"><sup>a</sup> the D prefixing the number is a recent addition to distinguish the classes from another numerical avalanche size classification also used in North America.</td></tr></tfoot></table>

Based on Canadian data for recreationists caught in potentially fatal avalanches (Size D2 and larger), the unadjusted vulnerability V is 0.19 [5]. The vulnerability adjusted for R = 0.10 is V<sub>R</sub> = 0.03 (Eq. 4) as shown in Table 2.

Jamieson and Jones [3] identified three factors that contribute to this surprisingly low vulnerability:

* The recreationists are typically caught when on skis or a snowmobile high in an avalanche path (Figure 2).

* They have some skill to escape a moving avalanche.

* They wear avalanche transceivers, and are accompanied by other recreationists capable of companion rescue with transceivers, shovels and probes.

In the absence of these factors, vulnerability would be higher [3].

We recognize that the definition of “potentially fatal” based on the Canadian size classification [4] may contribute to the low vulnerability. This adjusted vulnerability does not apply to non-recreationists (Figure 3) and will not apply to countries with a different definition of a potentially fatal avalanche.

The vulnerability adjusted for non-reporting is difficult to validate. In a study of the efficacy of balloon packs, Haegeli et al. [6] reasonably assumed that most non-fatal involvements would be reported in fatal avalanches. However, the selection criteria included deployed balloon packs, excluding many of the smaller potentially fatal avalanches in which recreationists are often caught. Further, Haegeli et al. did not partition their results by avalanche size [4], which is required for comparison with the adjusted vulnerability in Table 2. We sought individual well documented avalanches with numerous people caught

and at least one fatality, and found one such avalanche. The Size D3 avalanche on 2013-03-13 about 18 km west of Revelstoke, BC caught 50 to 100 people in the runout zone, injured 32 and killed 2 [7], making it suitable for comparison with the adjusted vulnerability V<sub>R</sub>. In this single avalanche with many involvements, 0.02 to 0.04 of those caught were killed. While vulnerability from a single avalanche cannot be used to validate the vulnerability, the fatality rate is not inconsistent with vulnerability V<sub>R</sub> = 0.07 for a Size 3 avalanches based on R = 0.10 (Table 2).

<table><caption>**Table 2.** Avalanche vulnerability for recreationists in Canada by avalanche size, 1984-2011 [5]</caption><thead><tr><th rowspan="2">Avalanche size [4]</th><th rowspan="2">Relative frequency of recreationists caught (n = 1343)</th><th colspan="2">Probability of death if caught, V<sub>R</sub></th></tr><tr><th>R = 1</th><th>R = 0.10</th></tr></thead><tbody><tr><td>D2</td><td>0.55</td><td>0.07</td><td>0.007</td></tr><tr><td>D2.5</td><td>0.22</td><td>0.20</td><td>0.02</td></tr><tr><td>D3</td><td>0.17</td><td>0.43</td><td>0.07</td></tr><tr><td>D3.5</td><td>0.06</td><td>0.63</td><td>0.15</td></tr><tr><td>Frequency-weighted vulnerability</td><td></td><td>0.19</td><td>0.03</td></tr></tbody></table>

<figure><img src="image_206.png" /><figcaption>**Figure 2.** When caught, recreationists are typically caught near the top of an avalanche such as this whereas non-recreations are typically caught in the runout of the avalanche. B. Jamieson photo.</figcaption></figure>
<figure><img src="image_207.png" /><figcaption>**Figure 3.** This person standing at the toe of a stopped avalanche illustrates two common characteristics of a non-recreationist caught in an avalanche, e.g. caught low in path, and not on skis or snowmobile. B. Jamieson photo.</figcaption></figure>

In a 2012 study [8], Jamieson et al. assumed a reporting rate for non-fatal involvements of R = 0.05 to 0.10 and used the adjusted vulnerability to calculate the risk to backcountry skiers for various levels of avalanche danger. While this required simplifying assumptions for exposure and terrain, it illustrates an application of adjusting vulnerability for reporting rate in risk calculations.

## 4 Summary

For risk-based design or risk assessments, vulnerability can be based on records of damage, expert judgement, or damage models. When based on records of damage, incomplete recording of events that did not cause specified levels of damage can lead to overestimation of vulnerability. We are unsure how many vulnerability calculations, including those relevant to snow hazards, are based on records of damage. This paper derives an equation to adjust vulnerability for incomplete recording of specified levels of damage.

We demonstrate the vulnerability adjustment using recreationists caught in snow avalanches, for which underreporting of non-fatal involvements is common in Canada. The adjusted vulnerability is lower than most previously published values. This is attributed to underreporting of non-fatal involvements as well as factors relating to recreational exposure to avalanches, the skill of victims to sometimes escape a moving avalanche, and rescue capability by suitably equipped companions. The definition of “potentially fatal” varies between countries and may also affect the adjusted vulnerability. The vulnerability to snow avalanches of people outside buildings or on roads will be higher.

For avalanche vulnerability calculations, we propose that recreational activity be analyzed separately from non-recreational activity. Avalanche involvements during recreation or avalanche mitigation work are characterized by

1. the victim or another person in the group of recreationists triggers the avalanche

2. the victim is usually caught high in the path (Figure 2)

3. the victim is usually on skis, snowboard or snowmobile

4. the victim has some training and skill to escape the moving avalanche

5. the victim is usually wearing an avalanche transceiver and surviving members of the group have transceivers, avalanche probes, shovels and skill in companion rescue.

None of these five characteristics are common for non-recreational exposure to avalanches (Figure 3). While the contribution of the individual factors to vulnerability has received little attention, this would be worthwhile research.

## References

1. D. Bertrand, M. Naaim and M. Brun. Nat. Hazards Earth Syst. Sci., **10**, 1531-1545 (2010)

2. European Commission. *The Design of Avalanche Protection Dams. Practical and theoretical developments and results*. Ed. by T. Jóhannesson, D. Issler, P. Gauer and K. Lied, European Commission, EUR 23339 (2009)

3. B. Jamieson and A.S.T. Jones. The effect of underreporting of non-fatal involvements in snow avalanches on vulnerability, *12th International Conference on Applications of Statistics and Probability in Civil Engineering*, June 2015, Vancouver, BC, Canada (2015)

4. D.M. McClung and P.A. Schaerer. Snow avalanche size classification. Canadian Avalanche Committee. Avalanche Workshop, 3-5 Nov. 1980. Associate Committee on Geotechnical Research. Technical Memorandum 133, National Research Council of Canada, 12-27 (1981)

5. B. Jamieson and A.S.T. Jones. Vulnerability: Caught in an avalanche - then what are the odds? *Proceedings of the 2012 International Snow Science Workshop in Anchorage, Alaska*, 1-8 (2012)

6. P. Haegeli, M. Falk, E. Proctor, B. Zweifel, F. Jarry, S. Logan, K. Kronholm, M. Biskupic and H. Brugger. Resuscitation 85, 1197-1203 (2014)

7. J. Kelly. Avalanche.ca - The Journal of Canada's Avalanche Community **93**, 20-21 (2010)

8. B. Jamieson, J. Schweizer, J., and C. Shea. Simple calculations of avalanche risk for backcountry skiing, *International Snow Science Workshop*. Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Davos, Switzerland, pp. 336-340 (2009)
# Seismic responses of snowpack on a slope in a vibration experiment

Yusuke Harada<fnref n="1" />,ᵃ, Wataru Takahashi<fnref n="1" />, Satoshi Omiya<fnref n="1" />, Hiroki Matsushita<fnref n="2" />, Takahiro Chiba<fnref n="3" /> and Masaru Matsuzawa<fnref n="1" />

**Abstract.** A vibration experiment was conducted to understand the seismic response of the snowpack on a slope. The following were found. 1) The dynamic magnification factor of a snowpack on a slope increases linearly with increase in the height of the snowpack when the input seismic motion has a short period. 2) Regardless of the snow grain shape, the natural mode of a snowpack on a slope is the primary mode. 3) The dynamic magnification factor is greater for melt forms than for rounded grains, and is greater for wet snow than for dry snow.

## 1 Introduction

When a large-scale earthquake occurs, it is necessary to know the conditions of damage and to quickly start rescue and evacuation operations. During the snowy season, avalanches caused by earthquakes have the potential to obstruct information-gathering activities and rescue operations. Recently in Japan, roads were closed at many locations by earthquake-induced avalanches in the northern area of Nagano Prefecture in March 2011 [1] and in the northern area of Tochigi Prefecture in February 2013 [2]. It is very important to prepare for such situations under the assumption that earthquakes can cause avalanches. The local disaster prevention plans of some Japanese municipalities clearly recognize the potential for avalanches to occur during and after an earthquake in the snowy season. However, the snowpack and seismic motion conditions that assess the risk of earthquake-induced avalanche have not been clarified. To enhance disaster prevention and damage mitigation with respect to earthquake-induced avalanches, clarification of the occurrence mechanism of earthquake-induced avalanches and proposals for a risk assessment method for such avalanches are thought to be useful.

There have been many reports of avalanche damage following earthquakes in the U.S.A., Russia, the western Himalaya mountains, Switzerland, Japan and elsewhere [1-7]. To clarify the conditions under which earthquake-induced avalanches occur, Podolskiy et al. (2010) [3] investigated the relationship between the location of earthquake-induced avalanches and the scale and distance from the epicenter of the earthquake by using the previous cases. For the earthquake-induced avalanche that occurred during the March 12, 2011, earthquake (M. 6.7) in northern Nagano Prefecture, Kamiishi et al. (2012) [1] conducted a snow-profile observation in the starting zone of the avalanche, and based on the observations, examined the stability of the snowpack on slopes at the time of earthquakes. Matsushita et al. (2013) [2] conducted observations and examinations similar to those of Kamiishi et al. for an earthquake-induced avalanche that occurred after the February 25, 2013 earthquake in northern Tochigi Prefecture. Matushita et al. (2014) [8] examined cases of earthquake-induced avalanches for which the observation data of the

snowpack around the avalanche occurrences were reported [1, 2, 5, 6] and investigated the conditions under which earthquake-induced avalanches tend to occur based on the concept of the stability of the snowpack on the slope. As for the stability of the snowpack on a slope, Chernouss et al. (2002) [9], Matsuzawa et al. (2007) [10] and Pérez-Guillén et al. (2014) [7] conducted theoretical examinations in which seismic motions were considered. Chiba et al. (2015) [11] conducted a seismic response analysis in which the snowpack on the slope was modeled by a particle system shear model.

There have been reports on analyses based on vibration experiments using snow blocks that reproduced the snowpack on a slope during an earthquake. Podolskiy et al. (2010) [12] demonstrated that snowpack failure resulting from seismic motion is initiated when tensile failure occurs in a weak layer. Kamiishi et al. (2012) [13] examined the failure of a snowpack that included a weak layer and showed that it is possible for the horizontal vibration to affect the snowpack more than the vertical vibration does. Matsushita et al. (2013) [14] demonstrated that the acceleration in the direction in which the snow layers are subjected to tensile force was greater than the acceleration in the direction in which the snow layers are subjected to compressive force.

The above-stated assessment methods for the stability of the snowpack on the slope, which consider the seismic motion, assume that an earthquake-induced avalanche occurs when the ground surface acceleration generated by the seismic motion is the same as the acceleration of the snowpack on the slope. However, depending on the magnitude and frequency of the seismic motion, the snowpack may undergo greater acceleration than the ground surface does. In addition, the transfer of seismic acceleration in the snowpack may differ with differences in grain shape. It is thought to be necessary to clarify the basic relationship between the state of snowpack and the seismic motion. In this study, the authors conducted vibration experiments to investigate the characteristics of the seismic response of the snowpack on a slope. Based on the experimental findings, we proposed a method for estimating the risk of earthquake-induced avalanches in which snowpack stability is considered in static analysis. Our current study addresses only the lateral vibration, and

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

<footnote n="1">¹PWRI, Civil Engineering Research Institute for Cold Region, 1-3-1-34 Hiragishi, Toyohira-ku, Sapporo, 062-8602, Japan</footnote>

<footnote n="2">²PWRI, Snow Avalanche and Land Slide Research Center, 2-6-8 Nishiki, Myoko, Niigata, 944-0051, Japan</footnote>

<footnote n="3">³Hokkaido University of Science, 7-15-4-1 Maeda, Teine-ku, Sapporo, 006-8585, Japan</footnote>

<footnote>ª Corresponding author: harada-y@ceri.go.jp</footnote>
the avalanches are assumed to start as a result of shear fracture.

## 2 Experimental methods

The experiment was conducted at the Ishikari Test Field, Japan (N43°12′55″, E141°23′23″) from January 2014 to March 2015. The measured items were the acceleration of the snowpack specimen (hereinafter: "the specimen"), the acceleration of the shaking table, and certain physical properties of the specimen (snowpack height, and snow layer structure, grain size, density, hardness and temperature).

First, snow blocks were cut from natural snowpack to make the specimens. Each snow block was reshaped into a parallelogram specimen (angles: 60° and 120°) 40 cm wide by 55 cm long by 40 to 70 cm high. A slope model (50 cm wide by 55 cm long) with an angle of 30° was anchored to a shaking table (200 cm wide by 200 cm long: SPTDU-20K-85L-50T) installed inside a building. The specimen was placed on the slope such that the snow layer structure would parallel the slope (Figure 1). To prevent the underside of the snowpack from sliding, 5-cm nails were installed on the slope at 6-cm intervals in the downslope and cross-slope directions. In addition, the L-steel member on the acceleration (positive) side was fixed with a clear acrylic board and vises.

Next, 3 to 6 accelerometers (18 × 18 × 24 mm, 40g: ASW-5A) were inserted into the snowpack with a vertical spacing of 10 to 20 cm. For some experiments, we used specimens weighted with a 15-kg lattice weight (30 cm wide by 45 cm long by 5.3 cm high) fixed to the snowpack with 32 U-shaped fixtures each weighing 0.075 kg (8 cm wide by 15 cm long).

After the accelerometers were installed, the specimen was accelerated in the horizontal direction for 3 min (Figure 1). By means of sweep acceleration where the acceleration of the shaking table was fixed and the frequency was varied between 1 and 10Hz. The range of input acceleration to the shaking table was varied between 0.1G and 0.9G (1G = 9.81 m/s²). The accelerations of the specimen and the shaking table were measured every 0.04 sec.

## 3 Experimental results

### 3.1 Response characteristics of the specimens

The vibration experiment was conducted 194 times,

<figure><img src="image_209.png" /><figcaption>Figure 1. The shaking table and specimen.</figcaption></figure>

with combinations of specimen conditions and accelerations varied. In this experiment, the maximum value of the rate of acceleration of the snowpack on the shaking table (hereinafter: "the dynamic magnification factor" (DMF)) of snowpack on the slope was found to be about 10Hz (period: 0.10s) for all specimens. The DMF was higher for greater input accelerations and greater weights, and for higher locations in the snowpack varied from 1 to 10Hz. Figure 2 shows the relationship between frequency and DMF.

### 3.2 The DMF by grain shape and height

Using the results of this experiment, we compiled the following data on each specimen: snow grain shape, mean snow density (kg/m³), mean snow hardness (kN/m²), specimen weight (kg), snowpack depth (m), "weight converted snowpack depth" for the weighted specimen (the total of the specimen's weight and the weight is converted to snow depth), input acceleration (G) and the max. value of DMF (hereinafter: "max. DMF") measured by the accelerometers. We measured the max. DMF at various input accelerations and created 800 data profiles. In this study, we classified the snow type as "dry rounded grains", "wet rounded grains", "dry melt forms" or "wet melt forms", based on the 50% higher grain shape and dry or wet of the specimen (Fig.

<figure><img src="image_210.png" /><figcaption>Figure 2. DMF vs. frequency (example).</figcaption></figure>

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

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

<table><caption>Figure 3. The layered structure of the specimens, with dry or wet snow, and the date of the experiment.</caption><thead><tr><th>Class</th><th>symbol</th></tr></thead><tbody><tr><td>Precipitation Particles</td><td>+</td></tr><tr><td>Decomposing and Fragmented precipitation particles</td><td>/</td></tr><tr><td>Rounded Grains</td><td>●</td></tr><tr><td>Melt Forms</td><td>○</td></tr><tr><td>Faceted Crystals</td><td>□</td></tr><tr><td>Ice Formations</td><td>—</td></tr></tbody></table>
3). The accelerometers were installed at four height ranges: 10 - 20cm, 30 - 39cm, 40 - 49cm, and 50 - 60cm.

We performed single regression analyses to obtain correlation coefficients. The objective variable was the max. DMF for an input acceleration whose explanatory variable was one of following: mean snow density, mean snowpack hardness, specimen weight, snowpack depth, weight converted snowpack depth, if the specimen is a weighted one, or input acceleration (Figure 4). The analyses found a positive correlation between the heights of the accelerometers and the input acceleration (G), between the heights of the accelerometers and weight (kg), and between the heights of the accelerometers and the weight converted snowpack depth (m) in all types of specimens and for all levels of accelerometer height. For the correlations between accelerometer heights and snow depth (m), between accelerometer heights and mean density (kg/m³), and between accelerometer heights and mean hardness (kN/m²), some of the data points plot with positive correlations, others with negative correlations, and yet others with near zero correlations. We made multiple regression equations for each group of accelerometers grouped in accordance with the range of installed elevation from the underside of the snowpack (i.e., 10 - 20cm, 30-39cm, 40-49cm, and 50-60cm) for each snow type (hereinafter: "accelerometer elevation") (Table 1). Where, y is the max. DMF, $x_1$ is the input acceleration (G) and $x_2$ is the weight (kg). The weight converted snow depth was treated as being included in the weight, because the snow depth was determined based on the weight. Wet rounded grains was excluded from this analysis because too few cases with such snow quality were found.

When the DMF at the frequency of 10Hz is 100%, the average value of the increase ratio in the DMF was calculated by using the data for each grain shape at the 30 - 39cm accelerometer height and at 0.5Hz intervals (Table 2). For each grain shape, the acceleration at 1 Hz (period: 1.0 S) was 0%; in other words, the acceleration was about the same as that of the shaking table. However, when the frequency increased gradually and approached 10 Hz (period: 0.10 S), the acceleration gradually increased. The DMFs increased about 80 percentage points from 5 Hz (period: 0.20 S) to 10 Hz (period: 0.10 S), and about 25 percentage points between 9 Hz (period: 0.11 S) and 10 Hz. No marked differences in the tendency of increase based on the grain shape were found. The standard deviation in the rate of increase in DMF for melt forms was greater than the standard deviation in the rate of increase in DMF for rounded

<table><caption>**Table 1.** Relationship between the max. value of DMF (y), and input acceleration ($x_1$) and weight ($x_2$).</caption><thead><tr><th>Snow type</th><th>Accelerometer elevation</th><th>Mean elevation</th><th>N</th><th>Multiple regression equation</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="4">Dry rounded grains</td><td>10~20cm</td><td>15.7cm</td><td>140</td><td>y=0.116x<sub>1</sub>+0.003x<sub>2</sub>+0.944</td><td>0.452</td></tr><tr><td>30~39cm</td><td>30.3cm</td><td>97</td><td>y=0.190x<sub>1</sub>+0.005x<sub>2</sub>+0.919</td><td>0.449</td></tr><tr><td>40~49cm</td><td>40.7cm</td><td>91</td><td>y=0.257x<sub>1</sub>+0.007x<sub>2</sub>+0.868</td><td>0.487</td></tr><tr><td>50~60cm</td><td>54.4cm</td><td>77</td><td>y=0.399x<sub>1</sub>+0.007x<sub>2</sub>+0.844</td><td>0.407</td></tr><tr><td rowspan="4">Dry melt forms</td><td>10~20cm</td><td>15.0cm</td><td>51</td><td>y=0.261x<sub>1</sub>+0.005x<sub>2</sub>+0.808</td><td>0.513</td></tr><tr><td>30~39cm</td><td>30.0cm</td><td>44</td><td>y=0.513x<sub>1</sub>+0.008x<sub>2</sub>+0.676</td><td>0.778</td></tr><tr><td>40~49cm</td><td>40.8cm</td><td>49</td><td>y=0.546x<sub>1</sub>+0.009x<sub>2</sub>+0.639</td><td>0.523</td></tr><tr><td>50~60cm</td><td>50.4cm</td><td>39</td><td>y=0.677x<sub>1</sub>+0.011x<sub>2</sub>+0.524</td><td>0.539</td></tr><tr><td rowspan="4">Wet melt forms</td><td>10~20cm</td><td>15.0cm</td><td>55</td><td>y=0.357x<sub>1</sub>+0.006x<sub>2</sub>+0.761</td><td>0.582</td></tr><tr><td>30~39cm</td><td>30.0cm</td><td>48</td><td>y=0.577x<sub>1</sub>+0.009x<sub>2</sub>+0.630</td><td>0.554</td></tr><tr><td>40~49cm</td><td>40.0cm</td><td>48</td><td>y=0.725x<sub>1</sub>+0.011x<sub>2</sub>+0.543</td><td>0.607</td></tr><tr><td>50~60cm</td><td>50.0cm</td><td>41</td><td>y=1.043x<sub>1</sub>+0.012x<sub>2</sub>+0.442</td><td>0.629</td></tr></tbody></table>

Note: Wet rounded grains was excluded from this analysis due to the small number of such samples.

<figure><img src="image_213.png" /><figcaption>**Figure 4.** The average elevations of the accelerometers and the relationship between the maximum DMF and the correlation coefficient of each measured items. See Table 1 for the average elevation of accelerometer and N for each snow type.</figcaption></figure>
grains. The relationship between grain shape and increase in the DMF relative to the frequency increase was determined as shown in Table 3.

## 3.3 Trial calculation for DMF of snowpack on a slope

We calculated the max. DMF, y, by substituting an arbitrary input acceleration for $x_1$ and by substituting the mean snowpack weight for each snow type for $x_2$. The snowpack weight was calculated as the product of volume and density. The volume of snow for each snow

<table><caption>**Table 2.** Increase ratio for DMF when the DMF at 10 Hz is normalized to 100 %.</caption><thead><tr><th rowspan="3">Period (S)</th><th rowspan="3">Frequency (Hz)</th><th colspan="2">Dry rounded grains</th><th colspan="2">Dry melt forms</th><th colspan="2">Wet melt forms</th></tr><tr><th colspan="2">N = 97</th><th colspan="2">N = 44</th><th colspan="2">N = 48</th></tr><tr><th>Avg.</th><th>Std Dev.</th><th>Avg.</th><th>Std Dev.</th><th>Avg.</th><th>Std Dev.</th></tr></thead><tbody><tr><td>1.000</td><td>1</td><td>0.0%</td><td>0.0%</td><td>0.0%</td><td>0.0%</td><td>0.0%</td><td>0.0%</td></tr><tr><td>0.667</td><td>1.5</td><td>0.7%</td><td>0.7%</td><td>0.8%</td><td>0.5%</td><td>0.8%</td><td>0.6%</td></tr><tr><td>0.500</td><td>2</td><td>2.4%</td><td>1.6%</td><td>2.7%</td><td>1.4%</td><td>2.3%</td><td>1.2%</td></tr><tr><td>0.400</td><td>2.5</td><td>4.4%</td><td>2.0%</td><td>4.6%</td><td>1.8%</td><td>4.1%</td><td>1.6%</td></tr><tr><td>0.333</td><td>3</td><td>6.7%</td><td>2.5%</td><td>6.7%</td><td>2.3%</td><td>6.2%</td><td>2.1%</td></tr><tr><td>0.286</td><td>3.5</td><td>9.6%</td><td>2.9%</td><td>9.3%</td><td>2.7%</td><td>8.7%</td><td>2.7%</td></tr><tr><td>0.250</td><td>4</td><td>12.9%</td><td>3.2%</td><td>12.3%</td><td>3.3%</td><td>11.6%</td><td>3.2%</td></tr><tr><td>0.222</td><td>4.5</td><td>16.8%</td><td>3.3%</td><td>15.7%</td><td>3.8%</td><td>15.0%</td><td>3.9%</td></tr><tr><td>0.200</td><td>5</td><td>21.1%</td><td>3.6%</td><td>19.7%</td><td>4.5%</td><td>18.8%</td><td>4.6%</td></tr><tr><td>0.182</td><td>5.5</td><td>26.0%</td><td>3.9%</td><td>24.1%</td><td>5.1%</td><td>23.1%</td><td>5.4%</td></tr><tr><td>0.167</td><td>6</td><td>31.3%</td><td>4.3%</td><td>29.2%</td><td>5.8%</td><td>27.9%</td><td>6.1%</td></tr><tr><td>0.154</td><td>6.5</td><td>37.4%</td><td>4.5%</td><td>34.8%</td><td>6.4%</td><td>33.4%</td><td>6.8%</td></tr><tr><td>0.143</td><td>7</td><td>44.1%</td><td>4.8%</td><td>41.1%</td><td>6.8%</td><td>39.5%</td><td>7.4%</td></tr><tr><td>0.133</td><td>7.5</td><td>51.4%</td><td>4.9%</td><td>48.2%</td><td>7.2%</td><td>46.4%</td><td>7.8%</td></tr><tr><td>0.125</td><td>8</td><td>59.5%</td><td>4.8%</td><td>56.4%</td><td>7.3%</td><td>54.3%</td><td>8.1%</td></tr><tr><td>0.118</td><td>8.5</td><td>68.2%</td><td>4.3%</td><td>65.4%</td><td>7.2%</td><td>63.3%</td><td>7.9%</td></tr><tr><td>0.111</td><td>9</td><td>77.7%</td><td>3.8%</td><td>75.3%</td><td>6.1%</td><td>73.6%</td><td>6.9%</td></tr><tr><td>0.105</td><td>9.5</td><td>88.1%</td><td>2.7%</td><td>86.5%</td><td>4.1%</td><td>85.6%</td><td>4.7%</td></tr><tr><td>0.100</td><td>10</td><td>100.0%</td><td>0.0%</td><td>100.0%</td><td>0.0%</td><td>100.0%</td><td>0.0%</td></tr></tbody></table>

Note: The data for the accelerometer elevation of 30 - 39 cm
Avg.: Average, Std Dev.: Standard deviation

<table><caption>**Table 3.** Relationship between the increase ratio of DMF (y (%)) and the frequency (x), when the DMF at 10 Hz is 100 %.</caption><thead><tr><th>Snow type</th><th>N</th><th>Relational equation</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>Dry rounded grains</td><td>97</td><td>y=0.041x<sup>2.41</sup></td><td>0.989</td></tr><tr><td>Dry melt forms</td><td>44</td><td>y=0.045x<sup>2.34</sup></td><td>0.992</td></tr><tr><td>Wet melt forms</td><td>48</td><td>y=0.041x<sup>2.37</sup></td><td>0.994</td></tr></tbody></table>

type was obtained by multiplying the mean cross-sectional area of the specimen at the time of the experiment (0.176m²) by the mean value of converted snowpack depth, which was 0.689 m for the dry rounded grains, 0.658 m for the dry melt forms and 0.624 m for the wet melt forms. Density was set in accordance with the *Japanese Dictionary of Snow and Ice* as follows [15]: 325 kg/m³ for dry rounded grains, 400 kg/m³ for dry melt forms and 500 kg/m³ for wet melt forms. Consequently, the values calculated for normal snowpack weight were as follows: 39.4 kg for dry rounded grains, 46.3 kg for dry melt forms and 54.9 kg for wet melt forms.

From the normal snowpack weights and the multiple-regression equation for arbitrary acceleration (Table 1), the DMFs were obtained. Figure 5 plots the DMF for each accelerometer elevation. A linear approximation was made for the DMF plotted for each input acceleration. The results show that the DMF increases proportionally to the accelerometer elevation. The acceleration is higher for the accelerometers at higher elevations. Consequently, regardless of the type of snow, the natural mode of specimens that assumed snowpack on slope is the primary mode. At the acceleration of 0.2G, the distribution of DMFs for melt forms has roughly the same shape as that for rounded grains. At accelerations of 0.4G or greater, the DMF is greater for melt forms than for rounded grains, and the DMF is greater for wet snow than for dry snow.

## 3.4 Difference in DMF by grain shape

Fig. 6 shows the DMF by weight of specimen at the input acceleration of 0.5 G at the height of 30 - 39cm, and the relationship between the DMF and the mean snow hardness and between the DMF and the mean density. In this result it is found that, for the dry rounded grains, the DMFs are high for high mean density, hardness, and weight. In the dry melt forms and wet melt forms, the above relationship is not found. Many of the data points for the dry melt forms and wet melt forms are high in density and DMFs than those for rounded grains. We estimated the cause of the differences in the DMFs between dry and wet snow in the results discussed in Section 3.3. A snowpack behaves as an elastic body when

<figure><img src="image_214.png" /><figcaption>**Figure 5.** Relationship between accelerometer elevation and DMF based on experiment (10Hz).</figcaption></figure>
the magnitude and duration of forces acting on it are small [16]. In this experiment, the snowpack was given shear force by vibration in the horizontal direction. Elasticity (or stiffness) against such input acceleration can be considered to behave as follows. Yamanoi and Endo (2002) [17] expressed the shear strength of precipitation particles, decomposing and fragmented precipitation particles, rounded grains, and dry melt forms in terms of a power function of density, and showed that the sheer strength of dry melt forms was smaller than those of other types of snow, and that the sheer strength of wet snow was smaller than that of dry snow. Kinoshita (1960) [18] expressed the hardness of dry snow in terms of a power function of density. For the increase in snow density, the dry rounded grains increases in hardness and shear strength. The rigidity against horizontal force increases when hardness and shear strength increase. Therefore, it is suggested that the DMF decreases. However, the weight also increases with the increase in hardness and shear strength, which is the reason for the great DMFs for rounded grains shown Fig. 6. In contrast, the dry melt forms has lower shear strength than dry rounded grains, so the shear strain of the dry melt forms against horizontal input acceleration is greater than that of dry rounded grains. Weight also effects shear strength. Therefore, the greater is the input acceleration, the greater is difference in DMF between dry melt forms and dry rounded grains. This trend is much more obvious for wet snow than for dry snow (Figure 5).

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

Note: This is data for the accelerometer elevation of 30 - 39 cm and the input acceleration 0.5 G, without weighting.

<caption>Figure 6. Relationship between DMF and the mean density, and between DMF and the hardness.</caption>

# 4. Estimation of the likelihood of earthquake-induced avalanches

First, the snowpack stability on the slope (SI: stability index) was assessed. Then, the risk of avalanche occurrence was estimated by using the results of the vibration experiment in Chapter 3 and the snowpack data measured at the time of previous avalanches.

## 4.1. Stability of the snowpack on the slope

The stability of the snowpack on the slope in nature SI is expressed as a ratio of the shear strength Σs of the subject snow layer to the shear stress σn sinψ acting on the snow layer (Equation (1) and Fig. 7 (a)). The smaller is the SI value, the less stable is the snowpack on the slope. The SI value is used for assessing the stability of the snowpack on the slope and for indicating the occurrence of avalanches [19], [20].

$$ SI = \frac{\Sigma_s}{\sigma_n \sin \psi} = \frac{C + \sigma_n \cos \psi \tan \phi}{\sigma_n \sin \psi} \quad (1) $$

Where, Σs is the shear strength (N/m²) of the subject snow layer, such as a weak layer; and σn is the snowpack load (N/m²) per unit volume above the weak layer as determined from the thickness D (m), density ρ (kg/m³), and gravitational acceleration g (m/s²) of the layers above the weak layer (σn=gρD). ψ is the slope gradient (°). The right-hand side of Equation (1) is the shear strength Σs expressed in terms of the Mohr-Coulomb yield criterion. Maeno and Kuroda (1986) [16] maintained that when the density of snow is low or the vertical load σn is great, application of this criterion does not produce accurate results because of the compressibility of snow. However, Nakamura et al. (2010) [21] demonstrated that the criterion held well even for low-density snow. Therefore, the authors of this study assume that the Mohr-Coulomb yield criterion holds for all types of snow. C is the cohesive force (N/m²) of snow particles. In this study, a measured value of shear frame index SFI is used. tanφ is the internal friction coefficient of the snow particles.

For actual snowpack on a slope, the critical value for avalanche occurrence is not necessarily 1 at the stability expressed in Equation (1), because other factors such as lateral connection between the sides of the snowpack and the sides of the surrounding snow cover influence the stability of the snowpack on the slope SI [19]. To clarify the critical value for avalanche occurrence, examinations for critical stability index SI have been done using actual cases of avalanches. For example, in Canada, a value of about 1.5 was presented as the threshold value of stability indicator SI for avalanche occurrence [22]. A study on avalanche occurrences at roadside slopes in Hokkaido, northern part of Japan, in which a snowpack model was used to investigate the snowpack stability index SI at the time of avalanche occurrence [20], found that avalanches tend to occur when the snowpack stability index SI is 2.5 or lower and that the number of avalanches markedly increases when SI falls to 2.0 or lower. Based on the
above, a stability index of 2.0 or lower is considered to be the threshold for avalanche occurrence in this study, and an SI of 1.5 - 2.0 is classified as "small risk", 1.0 - 1.5 is classified as "middle risk", and 1.0 or lower is classified as "great risk".

## 4.2 Stability of the snowpack on the slope considering seismic motion

The stability of the snowpack on the slope at the time of an earthquake SI<sub>E</sub> is expressed by Equation (2), which is obtained by adding the horizontal seismic coefficient a (Fig. 7 (b)), as an external force, to Equation (1) [1], [2], [3], [10].

$$SI_E = \frac{C + \sigma_n (\cos \psi - a \sin \psi) \tan \varphi}{\sigma_n (\sin \psi + a \cos \psi)} \quad (2)$$

The horizontal seismic coefficient a of Equation (2) is a ratio of the horizontal acceleration of seismic motion (gal) to the gravitational acceleration (gal), which is a horizontal acceleration (G). The horizontal seismic coefficient a is used for assessing the stability of the embankment at the time of earthquake [24] and for examining the aseismic performance of road facilities [25]. Therefore, we used the horizontal seismic coefficient a in examining the conditions for avalanche occurrence at the time of earthquake.

When we assume that the natural mode of snowpack on the slope, which was determined in Chapter 3, holds true even for snow depths greater than those in our vibration experiment, the horizontal seismic coefficient a' of an arbitrary snow depth H' and frequency f is expressed by Equation (3).

$$a' = a\{(SR_{10Hz} - 1)IR + 1\} \quad (3)$$

Where, SR<sub>10Hz</sub> is the DMF (Table 1, Fig. 5) at the height H' of the snowpack on the slope at 10 Hz, and IR is the rate of increase in DMF for a given frequency when the DMF for 10 Hz is regarded as 100% (Table 3).

In addition to proposing Equation (2), Matsuzawa et al. (2007) [10] proposed the stability index of the snowpack

on the slope as SI<sub>E</sub>’. In their proposal, the seismic motion and the tensile force (tensile fracture strength) that is generated by the bonding between snow particles and that acts in the slab on the slope were considered. The index was proposed based on the safety factor assessment of embankment slopes at the time of earthquake. In this study, we propose Equation (4), in which the above-mentioned horizontal seismic coefficient a' is considered in the equation of Matsuzawa et al.

$$SI_E' = \frac{CL + \sigma_n L(\cos \psi - a' \sin \psi) \tan \varphi + \Sigma_t D}{\sigma_n L(\sin \psi + a' \cos \psi)} \quad (4)$$

Where, L is the length of snowpack on the slope (m), Σ<sub>t</sub> is the tensile fracture strength of the snow (N/m²), and D is the thickness of the snow layer above the weak layer (m).

## 4.3 Stability of the snowpack on the slope at avalanche occurrence

To investigate the occurrence conditions of earthquake-induced avalanches, Matsushita et al. (2014) examined earthquake-induced avalanches for which snow cover observation data, including avalanche type and snowpack data were available [8]. In the current study, we use four cases of surface avalanches from the report of Matsushita et al. whose peak ground acceleration (gal) and frequency (Hz) were measured by the Strong Motion Seismograph Networks (K-net, KiK-net) of Japan's National Institute for Earth Science and Disaster Prevention, and we perform a trial calculation for the stability of snowpack on the slope. (See Introduction and Table 4.) The snow layers above the sliding surface in the earthquake-

<figure><img src="image_216.png" /><figcaption>Figure 7. Forces acting on the snowpack on the slope.</figcaption></figure>

(a) Load and strength of the snowpack on the slope (b) Horizontal seismic coefficient of seismic motion.

<table><caption>Table 4. Conditions of the earthquake-induced avalanches.</caption><thead><tr><th rowspan="2">Cases</th><th rowspan="2">Date</th><th colspan="3">Earthquakes</th><th rowspan="2">Distance between the Strong Motion Observatory and the avalanche (km)</th><th rowspan="2">Avalanche type</th><th rowspan="2">Slope Angle</th><th rowspan="2">Length</th><th colspan="2">Bed surface</th><th colspan="4">Snowpack above the bed surface</th><th colspan="3">Snowpack below the bed surface</th><th rowspan="2">Snow-profile observation station</th></tr><tr><th>Magnitude (M)</th><th>Horizontal seismic coefficient a ≈ 1</th><th>Frequency f (Hz)</th><th>Grain type</th><th>Cohesive force</th><th>Load</th><th>backbone</th><th>Depth</th><th>Density</th><th>Main snow type</th><th>Depth</th><th>Density</th><th>Main snow type</th></tr></thead><tbody><tr><td>2011 northern Nagano Prefecture #1</td><td>2011/3/12</td><td>6.7</td><td>0.72</td><td>11.2</td><td>15</td><td>Surface Slab</td><td>35</td><td>10<sup>9/2</sup></td><td>Melt Forms</td><td>1580</td><td>1178</td><td>0.61</td><td>0.75</td><td>197<sup>9/2</sup></td><td>Melt Forms</td><td>1.62</td><td>197<sup>9/3</sup></td><td>Wet melt forms</td><td>Mochi, Tokamachi, Nügata Prefecture</td></tr><tr><td>2011 northern Nagano Prefecture #2</td><td></td><td></td><td></td><td></td><td>6</td><td>Surface Slab</td><td>25</td><td>10<sup>9/2</sup></td><td>Melt Forms</td><td>780</td><td>1813</td><td>0.56</td><td>0.62</td><td>330<sup>9/2</sup></td><td>Melt Forms</td><td>2.04</td><td>330<sup>9/3</sup></td><td>Wet melt forms</td><td>Nonaka, Tokamachi, Nügata Prefecture</td></tr><tr><td>2013 northern Tochigi Prefecture #1</td><td>2013/2/25</td><td>6.3</td><td>1.25</td><td>4.3</td><td>4</td><td>Surface Slab</td><td>42</td><td>9.5</td><td>Faceted Crystals</td><td>1510</td><td>744</td><td>0.41</td><td>0.55</td><td>185</td><td>Rounded Grains</td><td>0.99</td><td>316</td><td>Dry melt forms</td><td>Okukimi, Nikko, Tochigi Prefecture</td></tr><tr><td>2014 northern Tochigi Prefecture #2</td><td></td><td></td><td></td><td></td><td></td><td>Surface Slab</td><td>42</td><td>9.5</td><td>Faceted Crystals</td><td>2890</td><td>1626</td><td>0.74</td><td>1.00</td><td>224</td><td>Rounded Grains</td><td>0.66</td><td>224</td><td>Dry melt forms</td><td>Okukimi, Nikko, Tochigi Prefecture</td></tr></tbody></table>

Note 1: The horizontal seismic coefficient a is the maximum value, and the frequency f is that for the max.

Note 2: The value is from the figure in the reference [1].

Note 3: The density of the snow in the upper section of the snowpack was used as a substitute, because we did not have the data for this [1].
induced avalanche in northern Nagano Prefecture were melt forms [1], and those in northern Tochigi Prefecture were mainly rounded grains [2]. Therefore, we use a relational equation of tensile fracture strength (Eqs. (5a), (5b)) for each snow type developed by Watanabe (1977) [25].

$$ \Sigma_t = 7.78 \times 10^{-3} \rho^{2.60} \quad (\text{Melt Forms}) \quad (5a) $$

$$ \Sigma_t = 3.40 \times 10^{-4} \rho^{3.24} \quad (\text{Rounded Grains}) \quad (5b) $$

The stability of snowpack on the slope for each of the earthquake-induced avalanches was calculated based on the above. The left graph in Fig. 8 compares the stability $SI_E'$ in which the horizontal seismic coefficient $a'$ is used to consider the DMF (Equation (4)) and the stability $SI$ in which the $a'$ is not considered (Equation (1)). There were cases with $SI$ (shown as the horizontal axis) of 1.5 or lower (i.e., the mid to high degree of avalanche occurrence risk). It was thought that, in these cases, the snowpack on the slope was in an unstable condition before the earthquake [1]. For cases with $SI$ of 1.5 or higher, the $SI_E'$ is 1.5 or lower. It is thought that, in these cases, the snowpack had been relatively stable before the earthquake, and when the external force of the earthquake motion was exerted on the snowpack, it became unstable and avalanches occurred. The right graph in Fig. 8 compares the $SI_E'$ determined by using Equation (4) and the $SI_E'$ determined by substituting the horizontal seismic coefficient $a$ of the seismic motion into Equation (4) (i.e., the equation of Matsuzawa et al. [10]). The difference between the two is shown as the difference between the calculated $SI_E'$ values. In the earthquake in northern Nagano Prefecture (frequency: 11.2 Hz), the decrease in #1 was 0.95 and the decrease in #2 was 0.99. In the earthquake in northern Tochigi Prefecture (frequency: 4.3 Hz), the decrease in #1 was 0.15 and the decrease in #2 was 0.11. The above demonstrates that, for seismic motion with short periods and for snowpack whose snow depth $H'$ below the weak layer is great (Table 4), the decrease in stability is great.

## 4.4 Trial calculation of the risk of avalanche occurrence based on snow depth and horizontal seismic coefficient

When expressing avalanche occurrence likelihood, the use of the stability of snowpack on the slope as an index is not practical, because it is necessary to estimate the stability of snowpack on the slope based on the snowpack observation data each time. We had assumed that the risk of avalanche occurrence increases in stages when the snowpack stability of the slope $SI_E'$ is 2.0 or lower, and we conducted a case study to show the avalanche occurrence risk at the time of earthquake. In our study, which is based on Equation (4) and considers snow type, the depth of snowpack on the slope $H_s$ (m) (Fig. 7 (a)) is used as an index. Here, we assume a dry surface slab avalanche. The model slope shown in Fig. 7 (a) is used, and the conditions are as follows: A layer of faceted crystals that can become a weak layer exists in the rounded grains layer, and the rounded grains layer above the weak layer starts to flow when seismic motion (horizontal motion) with a frequency of 10 Hz causes a shear fracture in the weak layer. The density of 325 kg/m³ is set for rounded grains (Sec. 3.3), and the density of 200 kg/m³ is set for faceted crystals by selecting the smallest value for the density range of 200 - 400kg/m³ [15]. For the slope gradient $\psi$, 40° is selected, because many dry surface slab avalanches occur on slopes with gradients of 30° - 45°, and the frequency peak of the occurrence distribution is near 40° [26]. In addition to the above settings, the variables for Equation (4) are set as follows.

The cohesive force of snow particles in the weak layer $C$ (N/m²) is determined by using the relational equation for shear strength $\Sigma_s$ which was measured by Watanabe [25] in the condition without snow load above the weak layer (Equation (6)). In this case, it is possible to regard $\Sigma_s = C$, because Equation (1) gives 0 for the snow load $\sigma_n$.

$$ C = 3.56 \times 10^{-5} \rho^{3.36} \quad (\text{Faceted Crystals}) \quad (6) $$

$C = 1918$ N/m² is obtained by substituting a density of

<figure><img src="image_217.png" /><figcaption>**Figure 8.** The calculation results for the stability of the snowpack on the slope.</figcaption></figure>

Left: Comparison between the stability index $SI_E'$ determined by using horizontal seismic coefficient $a'$, for which DMF was considered (Eq. (4)), and the stability index $SI$ determined by not considering the $a'$ (Eq. (1)).

Right: Comparison between the $SI_E'$ obtained by using Eq. (4) into which $a'$ was substituted and the $SI_E'$ obtained by using Eq. (4) into which $a$ was substituted.
200 kg/m³ for faceted crystals into Equation (6).

The internal friction coefficient tanφ is considered in the cases for precipitation particles, decomposing and fragmented precipitation particles, and faceted crystals, but is negligible for the other types of snow [17], [27]. In a previous study, the internal friction coefficient for precipitation particles tanφ = 0.46 - 2.0 [10], [28], that for decomposing and fragmented precipitation particles tanφ = 0.23 [11], and that for faceted crystals tanφ = 0.21 [29]. In our current study, to assume faceted crystals as the weak layer, 0.21 was set as the internal friction coefficient.

For the thickness of the snow layer above the weak layer D (m), we used 0.60 m, which was roughly the average value of the occurrence frequency distribution of slab thicknesses (i.e., the thickness above the weak layer) obtained from 200 dry snow slab avalanches examined by McClung and Schaerer [26].

For the snow load above the weak layer per unit volume σn (N/m²), σn = 1911 N/m² was obtained by multiplying 325 kg/m³, which is the density of the rounded grains layer, by the thickness of the snow layer above the weak layer D = 0.6 m and the gravitational acceleration g. The tensile fracture strength Σt (N/m²) for the snow above the weak layer was obtained as Σt = 46771 N/m², based on Equation (5b), which is an equation for the tensile fracture strength of rounded grains.

For the length of the snowpack on the slope L (m), 10 m and 20 m were used. The arbitrary snow depth H' and the horizontal seismic coefficient a' at frequency f were obtained by using Equation (3), Table 1 and Fig. 5. The arbitrary snow depth H' can also be determined by using Equation (7).

$$H'=H_s-\frac{D}{\cos \psi} \quad (7)$$

Fig. 9 and Fig.10 show the relationships between the snow depth Hs and the horizontal seismic coefficient of seismic motion a when the stability of snowpack on the slope SI_E' is 1.0, 1.5 and 2.0. The length of the snowpack on the slope in Fig. 9 is 10 m, and that in Fig. 10 is 20 m.

<figure><img src="image_218.png" /><figcaption>**Figure 9.** Relationship between snow depth Hs and horizontal seismic coefficient a of seismic motion. (Calculation conditions: frequency: 10 Hz; slope gradient: 40°; The snowpack: rounded grains; weak layer: faceted crystals; internal friction coefficient of the weak layer: tanφ = 0.21; length of the slope: 10 m.)</figcaption></figure>

The regions above the blue curves in the figures indicate the regions where the stability of snowpack on the slope SI_E' is 2.0 or lower. The curves in the graphs show that the snow depth and the likelihood of avalanche occurrence are proportional, and that when SI_E' becomes 1.0 or lower (red curves and the regions above them), the avalanche likelihood becomes very high. With increases in the horizontal seismic coefficient and the length of the snowpack on the slope, the snow depth that has to be considered in the risk assessment of avalanche occurrence became smaller.

From the above, it is thought that the risk of earthquake-induced avalanche occurrence can be indicated by using snow depth Hs and horizontal seismic coefficient a as indexes and by using Equation (4), for which variables are set based on the concept of stability of snowpack on the slope. When the weak layer is assumed to be near the bottom of the snowpack, it is also thought that risk assessment of the occurrence of full-depth avalanches is possible by giving an appropriate friction coefficient to the equation.

## 5. Summary and future issues

In this study, a vibration experiment was conducted that focused on the seismic response of the snowpack on the slope. The experiment found that the seismic response to short-period vibration increases proportionally to the increase in the snowpack height, that the DMF of melt forms is greater than that of rounded grains, and that the DMF of wet snow is greater than that of dry snow. We propose a method for estimating the risk of avalanche occurrence in which the snow depth at the time of earthquake Hs and the horizontal seismic coefficient a are used as indexes.

The conditions of the snowpack on the slope at the time of earthquake in winter vary depending by location and period in winter. When the snowpack on the slope is dry new snow due to snowstorm, it tends to be unstable and it may slide even when the snow depth and the horizontal seismic coefficient are smaller than the estimated index values. In the future, we will further investigate the risk of avalanche occurrences for various

<figure><img src="image_219.png" /><figcaption>**Figure 10.** Relationship between snow depth Hs and horizontal seismic coefficient a of seismic motion. (Calculation conditions: frequency: 10 Hz; slope gradient: 40°; The snowpack: rounded grains; weak layer: faceted crystals; internal friction coefficient of the weak layer: tanφ = 0.21; length of the slope: 20 m.)</figcaption></figure>
snow types and vibration frequencies, and further our
studies on the risk of earthquake-induced avalanches in
winter.

Acknowledgments

In writing this paper, we used data provided by the
Strong Motion Seismograph Networks (K-net, KiK-net)
of the National Institute for Earth Science and Disaster
Prevention. In conducting the vibration experiment, we
were assisted by many people. We would like to express
our gratitude to all who assisted us in this study.

References

1. Kamiishi, I., et al., *J. Jpn. Soc. Snow Ice* (Seppyo), **74**, 159-169(2012). [In Japanese with English abstract]

2. Matsushita, H., et al., *Proceedings of International Snow Science workshop Grenoble*, 1122-1129 (2013)

3. Podolskiy, E. A. et al., *Journal of Glaciology*, **56**, 431-446 (2010).

4. Kazakov, N. A., et al., *Mater. Glyatsiol. Issled*. (Data of Glaciological Studies), **88**, 102–106 (2000). [In Russian with English abstract]

5. Higashiura, M. et al., *Rep. Nat. Res. Cent. Disaster Prev.*, **21**, 103-112 (1979). [In Japanese with English abstract]

6. Ogura, Y., et al., *Ann. Rep. Saigai-ken, Niigata Univ.*, **23**, 9-15 (2001). [In Japanese with English abstract]

7. Pérez-Guillén, C., et al., *Cold Regions Science and Technology*, **108**, 149-162 (2014).

8. Matsushita, H., et al., *Monthly Report of Civil Engineering Research Institute for Cold Region*, **733**, 39-44, (2014). [In Japanese]

9. Chernouss, P., E., et al., *Proceedings of the International Snow Science Workshop Penticton, BC, Canada*, 25-30, (2002).

10. Matsuzawa, M., et al., *Proceedings of Snow Ice Hokkaido* (Hokkaido no Seppyo), **26**, 95-98 (2007).
[In Japanese]

11. Chiba, T., et al., *Proceedings of Snow Ice Hokkaido* (Hokkaido no Seppyo), **34**, 91-94 (2015). [In Japanese]

12. Podolskiy, E. A. et al., *Journal of Glaciology*, **56**, 447-458 (2010).

13. Kamiishi, I., et al., *Proceedings of 2012 cold region technology conference*, **28**, 49-52 (2012). [In Japanese]

14. Matsushita, H., et al., *Proc. JSSI & JSSE Joint Conference 2013 in Kitami*, 222 (2013). [In Japanese]

15. Japanese Society of Snow and Ice (Eds.), *Japanese Dictionary of Snow and Ice, 2nd edition*, 246 (2014).
[In Japanese]

16. Maeno, N. and Kuroda, T., *Structure and physical properties of snow and ice*, 162-164 (1986). [In Japanese]

17. Yamanoi, K. and Endo, Y., *J. Jpn. Soc. Snow Ice* (Seppyo), **64**, 443-451 (2002). [In Japanese with English abstract]

18. Kinoshita S., Low Temperature Science, Ser., **A19**, 135-146 (1960). [In Japanese with English abstract]

19. Japanese Society of Snow and Ice (Eds.), Guidebook of snow observations, 31-53, 79-96 (2010). [In Japanese]

20. Nishimura K., et al., *Proceedings of 2005 cold region technology conference*, **21**, 244-248 (2005).

21. Nakamura, T., et al., *J. Glaciol.*, **56**(196), 333-338 (2010).

22. Perla, R., *Canadian Geotechnical Journal*, **14-2**, 206-213 (1977).

23. Japan Road Association: *Verification of Embankment Stability against the Action of Seismic Motion, The Guideline for Road Earthwork: Embankment works*, 2010 Edition, 119 – 127 (2010). [In Japanese]

24. Japan Road Association: *Using Static Verification to Verify Aseismic Capacity, Specifications for Highway Bridges and in "Chapter V. Aseismic Design" of the Associated Instruction Manual*, 57 – 108 (2010). [In Japanese]

25. Watanabe, Z., *Sci. Rep. Fukushima Univ.*, **27**, 27-35 (1977).

26. McClung, D.M. and P. Schaerer, *The avalanche handbook, 3rd Edition* (2006).

27. Jamieson, J. B., and Johnston C. D., *Ann. Glaciol.*, **26**, 296-302 (1998).

28. Seki, S., et al., Preprints of the 1998 conference Japanese Society of Snow and Ice, 183 (1998).

29. Zeidler, A. and B.Jamieson, *Cold Regions Science and Technology*, **44**, 194-205 (2006)
# Friction along a slider on snow

Werner Nachbauer<fnref n="1" />,ª, Sebastian Rohm<fnref n="1" />, Christoph Knoflach<fnref n="1" />, Joost van Putten<fnref n="2" />, Michael Hasler<fnref n="2" />

**Abstract.** The purpose of this study was to analyse friction along a slider on snow. The experimental approach consisted of friction measurements of four identical sliders of different lengths of 10, 25, 75 and 100 cm at 5 different velocities. The measurement took place on a linear tribometer located inside a cold laboratory. The normal force was adapted to assure the same contact pressure for the 4 sliders. The local friction force along the theoretical 100 cm slider model was calculated by subtracting the friction force of the next shorter slider from the slider and assigning this difference to the slider segment by which they differed. The local coefficient of friction along the slider revealed high friction in the front part (0-10 cm), considerably lower one in the second segment (10-20 cm), constant or slightly increased one in the middle part (25-75 cm), and towards the end increased friction (75-100 cm). This may be attributed to dry friction in the front part, subsequently smaller friction due to lubrication by meltwater and finally increased capillary drag towards the end.

## 1 Introduction

Based on [1] it is known that a liquid water film between ski and snow accounts for the low friction of skis on snow. [2] proposed different friction mechanisms that are caused by different thicknesses of the melt water layer. Without a water film, or even with a very thin one, dry friction dominates. With increasing water film thickness, mixed friction (dry/wet) comes up and finally with even thicker water films, wet friction dominates. [3] described friction along a ski, where the water film builds up as the ski slides and hence different friction mechanisms dominate along the ski. Later findings about a quasi-liquid layer on ice crystals [4] suggest that additional mechanisms may be involved. No experimental results are available to support models of different friction regimes on a sliding ski. The purpose of this study was to develop a method to measure the local coefficients of friction (CoF) of different sections along a ski.

## 2 Method

The study was carried out on the linear tribometer of the Department of Sport Science of the University of Innsbruck (Fig. 1). It is a friction test device which allows the measurement of the friction force of whole skis [5]. The tribometer is located in a cooling chamber. Air and snow temperature were set to -8°C. Fresh artificial snow was filled into the tribometer trough, compressed and evened. The snow had an average grain size of 0.2 mm and a density of 450 kgm-3.

To determine the CoF along a ski, wooden sliders of 10, 25, 75 and 100 cm running lengths with a width of 4.9 cm were produced. Tips and tails were 7 cm long with a radius of 20 cm. The sliders were flat without camber. Ultra-high molecular weight polyethylene ski bases were attached to the sliders by double-sided tape. All bases were grinded and waxed with liquid paraffin in the same way. The normal force was adapted for the different sliders to assure a homogeneous contact

pressure of 5.7 kPa. Each slider was tested at speeds of 0.5, 1, 2, 5 and 10 ms⁻<fnref n="1" />. For each test, series of 15 consecutive runs were measured, with each series starting on a fresh snow track. The friction force was averaged over the third through the fifteenth run. The local friction force for the 100 cm slider model was calculated by subtracting the friction force of the next shorter slider from the slider and assigning this difference to the slider segment by which they differed. The local CoF values along the slider model were calculated as quotient of friction force and normal force on the segments.

<figure><img src="image_220.png" /><figcaption>**Figure 1.** The linear tribometer of the Department of Sport Science of the University of Innsbruck</figcaption></figure>

## 3 Results

Figure 2 shows the friction forces of the 4 sliders of different lengths at the velocities 0.5, 1, 2, 5 and 10 ms⁻<fnref n="1" />. Friction stayed nearly constant at velocities from 0.5 – 2 ms⁻<fnref n="1" /> and increased at higher velocities for all slider lengths.

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

<footnote n="1">¹Department of Sports Science, University of Innsbruck, Austria</footnote>

<footnote n="2">²Centre of Technology of Ski and Alpine Sports, University of Innsbruck, 6020 Innsbruck, Austria</footnote>

<footnote>ª Corresponding author: werner.nachbauer@uibk.ac.at</footnote>
<figure><img src="image_222.png" /><figcaption>**Figure 2.** Friction force $F_f$ vs. velocity for the 4 sliders of different lengths</figcaption></figure>

In Figure 3 the CoF for the 4 different slider lengths are presented. A considerable decrease of the CoF was measured by increasing the slider length from 10 to 25 cm. For the 75 cm slider no further reduction of the CoF and for the 100 cm slider a slight increase was detected. These observations were more pronounced for the higher velocities of 5 and $10 \text{ ms}^{-1}$.

<figure><img src="image_223.png" /><figcaption>**Figure 3.** Friction coefficient $\mu$ vs. slider length for different velocities</figcaption></figure>

In Figure 4 the local CoF along the theoretical slider model are shown for the five tested velocities. The CoF was high in the first segment (0-10 cm), decreased strongly in the second segment (10-20 cm) and then stayed about constant or increased slightly in the third segment (25-75 cm), and increased in the end segment (75-100 cm).

The overall CoF of the 100 cm slider calculated by integrating the local CoF over the slider length increased with increasing speed: 0.041 (0.5 $\text{ms}^{-1}$), 0.042 (1 $\text{ms}^{-1}$), 0.043 (2 $\text{ms}^{-1}$), 0.050 (5 $\text{ms}^{-1}$), and 0.058 (10 $\text{ms}^{-1}$).

## 4 Discussion

To our knowledge the only study on the influence of the runner length on friction is from the year 1955 [6]. In accordance with our results a clear reduction of the CoF with increasing length from about 10 to 40 cm was found. The author suggested that it is likely that a smaller runner has more contact per unit than a larger one. Another explanation may lie in the melting of snow by frictional heating. Since the dry friction portion in front of the slider does not depend on the slider length the impact of

dry friction on the total friction increases with decreasing slider length.

<figure><img src="image_224.png" /><figcaption>**Figure 4.** Local friction coefficients $\mu$ along the theoretical slider model at different velocities. Slider tip is at 0 cm and endpoint at 100 cm.</figcaption></figure>

Interestingly, the local CoF along the slider revealed a similar pattern as suggested by the ideas and theoretical considerations described in [3]. High friction in a small portion of the front of the slider due to dry friction, subsequently smaller friction due to lubrication by meltwater and higher friction towards the end due to the increase of capillary drag may cause the measured friction distribution. It has to be noted that contrary to a real ski the measured slider had a uniform longitudinal load distribution.

## References

1. F.P. Bowden, T.P. Hughes, The Mechanism of Sliding on Ice and Snow. Proceedings of the Royal Society of London Seriea A, Mathematical and Physical Sciences 172:280-298 (1939).

2. S.C. Colbeck, The kinetic friction of snow. Journal of Glaciology 34:78-86 (1988).

3. S.C. Colbeck, A review of the processes that control snow friction. US Army Cold Regions Research and Engineering Laboratory Monograph 92-2 (1992).

4. J.G. Dash, H.Y. Fu, J.S. Wettlaufer, The premelting of ice and its environmental consequences. Reports on Progress in Physics 58:115-167 (1995).

5. W. Nachbauer, P. Kaps, M. Hasler, M. Mössner, Friction Between Ski and Snow. In: Braghin F, Cheli F, Maldifassi S, Melzi S, Sabbioni E (eds.) The Engineering Approach to Winter Sports, pp. 17-32. Springer New York, Place Springer New York (2016).

6. R. Ericksson, Friction of Runners on Snow and Ice. US Army Cold Regions Research and Engineering Laboratory Report (1995).
# Gliding friction of back country climbing skins

Michael Hasler<fnref n="1" />,ª, Sebastian Rohm<fnref n="1" />, Christoph Knoflach<fnref n="1" />, Joost van Putten<fnref n="2" /> and Werner Nachbauer<fnref n="1" />,<fnref n="2" />

**Abstract.** The purpose of this study was to analyse the gliding friction behaviour of climbing skins for back country skiing. On a linear tribometer positioned inside a low temperature laboratory, 2 different skins were tested over a wide range of normal forces at two velocities. The friction test showed that the type of skin has a major influence on the gliding performance. The coefficient of friction decreased with increasing normal load and was not affected by the velocity.

## 1 Introduction

Climbing skins are called skins because the first back country skiers in Norway used marine seals' skin and fur attached to their ski bases to prevent sliding backwards while walking uphill. Nowadays they consist of a textile construction of synthetic and/or mohair material. The outer surface of the climbing skin has hairs or scales that grab the snow, preventing backward movement of the skis. When the skis are moved forward, these surfaces flatten out to allow some glide. There is a huge variety of skins on the market. Beside geometric cut and mechanical fixation mechanism to the ski, they vary in material (nylon, mohair and combinations), fibre length, mechanical binding of the fibres to the supporting textile (v- or w-shape) and additional (hydrophobic) surface treatment. Little scientific evidence exists on the quantitative differences of friction under realistic conditions. [1] measured the friction of small samples under laboratory conditions and performed gliding tests in the field. The coefficients of friction (CoF) were between 0.1 and 0.2 and the correlation of the skin rankings between field and lab tests was good, but not for all specimens. The performance of the different materials and construction techniques is commonly only assessed by experts of the field.

During the gliding phase of a step, normal force on and velocity of the ski vary strongly. According to Coulombs law the friction force is assumed to be linear proportional to the normal force and thus the proportional factor $\mu$ is constant. However, from prior studies it is known that CoF can depend on velocity (e.g. polyethylene [2]) and normal force (e.g. textiles [3]). Hence, in the present study, the gliding friction of skins was measured at different normal forces and different velocities.

## 2 Method

The study was carried out on the linear tribometer of the Department of Sport Science of the University of Innsbruck (Fig. 1). It is a friction test device which allows the measurement of the friction force of whole skis over a

length of up to 24 m [2]. The tribometer is located in a cooling chamber. Air and snow temperature was set to -5°C. Fresh artificial snow was filled into the tribometer trough, compressed and evened. The snow had an average grain size of 0.2 mm and a density of 470 kgm⁻³.

<figure><img src="image_225.png" /><figcaption>**Figure 1.** The linear tribometer of the Department of Sport Science of the University of Innsbruck.</figcaption></figure>

Two different skins were tested. One was designed for regular ski tourers with a 70% mohair and 30% nylon structure (Climb, Pomoca Sa, Suisse) and the other for competitors and experienced ski tourers (Race, Pomoca Sa, Suisse). The skins were attached to a wooden slider with a running length of 100 cm and a width of 5 cm. Tips and tails were 7 cm long with a radius of 20 cm. The sliders were flat without camber. The skins were tested in gliding direction at two speeds (2.5 ms⁻<fnref n="1" /> and 5 ms⁻<fnref n="1" />) at normal forces of 8, 17, 58, 145, and 266 N.

For each test, series of 5 consecutive runs were measured, with each series starting on a fresh snow track. The friction force was averaged over the second through the fifth run. For further analysis, the measurements from 2 to 5 were used to calculate a mean value. The CoF were calculated as quotient of friction force $F_f$ and normal force $F_n$.

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

<footnote n="1">¹Centre of Technology of Ski and Alpine Sports, University of Innsbruck, Austria</footnote>

<footnote n="2">²Department of Sports Science, University of Innsbruck, Austria</footnote>

<footnote>ª Corresponding author: michael.hasler@uibk.ac.at</footnote>
## 3 Results

Figure 1 shows the friction force $F_f$ versus the normal force $F_n$. The progressions showed a clearly non-linear behaviour for both skins. $F_f$ of the Race-skin was always lower than $F_f$ of the Climb-skin, whereat the difference increased with increasing normal load.

<figure><img src="image_227.png" /><figcaption>**Figure 2.** Friction force $F_f$ vs. normal load $F_n$ of two skins (Climb-skin: circle; Race-skin: triangle) at the velocities of 2.5 ms⁻¹ (full line) and 5 ms⁻¹ (dashed line).</figcaption></figure>

Figure 2 shows the effect of the normal load on CoF for the Climb- and Race-skin for the two measured velocities. The gliding friction of the Race-skin was clearly superior to the Climb-skin at all velocities and normal forces. Except for the lowest loads, CoF decreased with increasing normal load at both velocities. Regarding the different velocities, no differences in CoF were observed.

<figure><img src="image_228.png" /><figcaption>**Figure 3.** Coefficient of friction $\mu$ vs. the normal load $F_n$ (Climb-skin: Circle; Race-skin: Triangle) for the velocities of 2.5 ms⁻¹ (full line) and 5 ms⁻¹ (dashed line).</figcaption></figure>

## 4 Discussion

There is a clear difference of the CoF between the two tested skins. The climb skin is designed for regular ski tourers with a mixture of comparatively long mohair and nylon fibres, whereas the race skin uses exclusively short mohair fibres. The length of the fibre is empirically known to affect the friction of skins on snow. With

decreasing fibre length, the friction is thought to be reduced [1] which is in accordance with our measurements. However, the underlying mechanism needs to be investigated.

Except for very low loads, the CoF for both skins decreases with increasing normal load. The initial lower values of the CoF at a normal load of 8 N may be attributed to an incomplete contact of the fibres with the snow. With further increasing load, the fibres might align better to the motion direction due to the pressure and the motion relative to the snow and hence CoF decreases.

The velocity showed no effect on the CoF of the tested climbing skins. Velocity effects were studied for other materials like polyethylene and PTFE [4,5] or steel [6] and revealed different results at similar snow temperatures as we studied. Usually an increase of CoF with velocity is explained by an increase of a melt water film. The velocity independence and the high CoF may be an indication that in the friction of back country skins on snow dry friction is dominating.

## Acknowledgement

This study was carried out within the “Competence Centre Sports Textiles” funded by the Standortagentur Tirol (Innsbruck, Austria)

## References

1. M. Auer, *Die Alpen*, 12, 27 (2006)

2. W. Nachbauer, P. Kaps, M. Hasler, M. Mössner, *The Engineering Approach to Winter Sports*, (Springer New York) 17 (2016)

3. M. Hasler, K. Schindelwig, W. Nachbauer, *Book of Abstracts of the 21st International Congress on Ski Trauma and Skiing Safety 8.-13.3.2015 68* (2015)

4. M. Shimbo, Friction on snow of ski soles, unwaxed and waxed. Scientific study of skiing in Japan, 99 (1971)

5. E. Spring, Valtion teknillinen tutkimuskeskus. Metallilaboratorio (1988)

6. S. Rohm, M. Hasler, C. Knoflach, J. van Putten, S. Unterberger, K. Schindelwig, R. Lackner, W. Nachbauer, Tribology Letters 59:1 (2015).
Session 4

# Snow physics / Snow drift II
# Characterizing the snowpack stratigraphy and its mechanical stability with hardness profiles measured by the Avatech SP1

Pascal Hagenmuller¹,ª, Thibault Pilloix¹

¹Météo-France/CNRS, CNRM-GAME/CEN, 1441 rue de la piscine, 38400 St Martin d'Hères, France

**Abstract.** Hardness has long been recognized as a good predictor of snow mechanical properties and thus as an indicator of snowpack stability. However, the usual measurement of hardness with a ramsonde lacks in resolution and misses snowpack features which are critical to assess its stability. The new snow penetrometer Avatech SP1 appears as a promising alternative by combining handy manipulation and millimetric vertical resolution. In this paper, we evaluate the SP1 accuracy in terms of penetration resistance and depth by comparing the SP1 profile with the profiles measured by the ramsonde and a reference penetrometer, the SMP. The capacity of the SP1 to detect unstable snowpack configurations is further assessed by modelling the snow pack stability directly from the measured hardness profile, which is then compared to stability tests. A method to aggregate numerous hardness profiles of similar snowpack configurations is also presented and provides a practical tool to extract synthetic information for the avalanche forecaster out of the large potential amount of data provided by the SP1.

## 1 Introduction

<figure><img src="image_229.png" /><figcaption>**Figure 1:** Different snow penetrometers: (a) ramsonde, (b) Avatech SP1 and (c) Snow MicroPenetrometer (SMP).</figcaption></figure>

The ramsonde (Fig. 1a) is used since the 1930's by avalanche forecasting services worldwide to provide quantitative information about snow stratigraphy. It measures hardness, defined as the resistance to penetration of an object into snow. However, due to the limited resolution of the instrument, a ramsonde profile does not include information about very thin snowpack features such as weak layers, which are critical to assess the snowpack stability. A ramsonde profile has thus to be complemented by time-consuming and possibly subjective stability tests and manual stratigraphy.

Recently, new portable high-resolution digital penetrometers such as the Snow MicroPenetrometer (SMP)[1] and the Avatech SP1 (www.avatech.com) were developed (Figs. 1b, 1c). They benefit from very sensitive force sensors and are expected to capture thin snowpack features invisible for the ramsonde (Fig. 2). On the contrary to the SMP whose price and fragility limit its usage to research, the SP1 is particularly easy to use and to transport in the field. It has the potential to become a common tool among snow professionals such as ski patrollers, mountain guides and avalanche forecasters.

The goal of this study is to assess the pros and cons of the SP1 in the framework of avalanche risk forecasting and to present methods to exploit the data provided by this new instrument. The presented work is based on preliminary measurements conducted in the French Alps

during winter 2014-2015. It will be complemented by measurements with the latest version of the SP1 (SP2) and feedbacks of regional forecasting centres of Météo-France during winter 2015-2016.

## 2 Methods and results

### 2.1 Evaluation of SP1 in terms of penetration resistance and depth

<figure><img src="image_230.png" /><figcaption>**Figure 2:** Example of hardness vertical profiles measured by the SMP, the ramsonde and the SP1. The arrows indicate some depth lags.</figcaption></figure>

In order to evaluate the accuracy and resolution of the SP1 in terms of penetration resistance and depth, co-located SP1, SMP and ramsonde profiles were measured and compared. The layers detected by the SP1 appear to be translated with a variable depth lag in comparison to the layers detected by the SMP (Fig. 2). The SMP depth accuracy is better than 1 cm. These depth lags are thus due to errors of the SP1 depth sensor. The alternating soft and hard layers are correctly reported by the SP1. The higher vertical resolution of the SP1 enables to capture detailed features that are not visible in the ramsonde profile.

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

<footnote>ª Corresponding author: pascal.hagenmuller@meteo.fr</footnote>
The difference between two hardness profiles can be decomposed into depth lags and residual hardness differences. To allow quantitative comparison between hardness profiles, we used an algorithm that matches two profiles by adjusting the layer thicknesses that maximize the hardness correlation [2]. Applied to 14 co-located SMP/SP1 profiles, the algorithm estimates the SP1 depth error in the range [-20, 20] cm with a standard deviation of 7.8 cm and a mean relative hardness error of 30%. Further tests will be conducted to evaluate the impact of the snow stratigraphy on the SP1 accuracy.

## 2.2 Evaluation of SP1 in terms of modelled failure initiation and propagation indices

<figure><img src="image_232.png" /><figcaption>**Figure 3:** Modelled distribution of the ratio between shear stress due to a skier point load (arrow) and shear strength, in a heterogeneous 38°-inclined slab. This ratio represents the likelihood of failure initiation.</figcaption></figure>

In the framework of avalanche forecasting, the SP1 capability of detecting unstable snowpack configurations cannot be directly derived from the evaluation of the SP1 in terms of penetration resistance and depth. It must be specifically evaluated. Moreover, it is necessary to develop an automatic method to assess the stability of a given hardness profile in order to point out the most critical situations among numerous profiles to the avalanche forecaster [3].

To this end, we adapted the method developed by [4] to SP1 data. The main idea of this method is to derive the mechanical properties of each layer from the hardness and to feed a finite element slab model with these layer properties (e.g. Fig. 3). This model then computes the likelihood of failure initiation and propagation.

On going combined stability tests and hardness profiles measured with the SP1 and SMP will help to

evaluate the capability of the SP1 to detect avalanche-prone stratigraphy and the benefit of the very high resolution of the SMP. Indeed, force fluctuations measured with the SMP can be linked to microstructural parameters and are not accessible with the SP1.

## 2.3 Spatial variability

For the purpose of real time analysis of the snowpack stability, the important amount of data collected by digital penetrometers, such as the SP1, must be quickly synthesized into a characterization representative of the test site of given altitude and aspect ranges.

At decameter scale, most snow layers are continuously present in the snowpack and the observed variability mainly originates from the varying thickness of each layer [5]. Assuming that depth deviations explain most of the observed variability and do not affect intensive layer properties, the algorithm used for the SP1/SMP comparison (see 2.1) is adapted to match multiple SP1 profiles by minimizing their intra-set hardness variability [2].

We successfully applied this algorithm on SP1 profiles and tracked centimetric layers over 20 m with little hardness variability (Fig. 4). The apparent spatial variability due either to real spatial variability or depth measurement errors must thus not hinder the development of high vertical resolution point snow measurements, which can provide information relevant at larger scale. Moreover, the method provides a practical way to synthetize numerous hardness profiles into one representative profile and estimates of the spatial variability of the snowpack in terms of layer depth and hardness.

## References

1. M. Schneebeli. Cold Reg. Sci. Technol. 30, 101–114 (1999).

2. P. Hagenmuller, T. Pilloix. A new method to compare snow hardness profiles. Submitted.

3. G. Giraud. In ISSW 97–104 (1992).

4. B. Reuter, J. Schweizer, A. van Herwijnen. The Cryosphere 9, 837–847 (2015).

5. Sturm, M., Benson, C. Ann. Glaciol. 38, 253–260 (2004).

<figure><img src="image_233.png" /><figcaption>**Figure 4:** The intra-set variability of the initial profiles (a) is minimized by adjusting the layer thicknesses on each (b). The transformations to go from the transformed profiles (b) to the initial profiles (a) are shown on subfigure (c).</figcaption></figure>
# Improvement of requirements for modeling snowdrifts in wind tunnels based on the measurements at Harbin

Qingwen Zhang<fnref n="1" />, <fnref n="2" />, Guolong Zhang<fnref n="1" />, <fnref n="2" /> and Feng Fan<fnref n="1" />, <fnref n="2" />

**Abstract.** Snowdrifts around cuboid models were investigated in wind tunnel and outdoors, considering the serious problems caused by snowdrifts around buildings when formed on undesirable places. Outdoor measurements were carried out to investigate the characteristics of drifting snow and to obtain data for comparison with wind tunnel results. The outdoor measurements were conducted at Harbin China, with model dimensions 750: 600: 750mm. Wind direction and speed were measured every 1 minute. Snowdrift heights around models were surveyed in approximately 400 points and contour maps produced. Then, a simulation of wind induced snowdrifts around reduced scale models, by using sodium bicarbonate in a direct wind tunnel, was described. The measurement at Harbin was selected as the prototype. Experimental procedure was proposed based on Kind's similarity parameters. With the results comparison between models and prototypes, the required length of upstream fetch was discussed.

## 1 Introduction

As an appropriate profiler to simulate snowdrift problems around buildings, wind-tunnel experimental methods have been explored for many years. However, due to the limitation of the experimental instrumentation and the complexity of the two-phase flow, similarity criteria cannot all be satisfied simultaneously. Compromises must be made based on the objective of the experiment. Related research has been undertaken in several studies (e.g., Kind; Anno; Iversen; Delpech).

Kind [1, 2] identified the important similarity criteria by analyzing the saltation process and discussed the compromises in modelling procedures. Naaim et al. [3] reviewed the conditions for a correct simulation of snowdrifts in wind tunnels, or outdoor, by comparing requirements proposed by four scientists (i.e., Kind, 1975, 1985; Iversen, 1979, 1980, 1981, 1982, 1984; Anno, 1984; Tabler, 1980). Delpech [4] overviewed the snowdrift mechanisms, discussed the experimental procedures based on similarity parameters proposed in the scientific literature. The simulation of wind induced snow drifting around buildings by using artificial snow, in the Jules Verne climatic wind tunnel, was carried out later, based on the proposed criteria. Okaze et al. [5] carried out a series of wind tunnel measurements over a loose snow surface to investigate the characteristics of drifting snow phenomena in a non-equilibrium boundary layer and to obtain data for examining the numerical model of snow transport coupled with computational fluid dynamics (CFD). Furthermore, as the previous wind tunnel experiments on snowdrifts usually used artificial snow as a substitute of natural snow, which cannot fully reflect the motion of natural snow particles. LÜ XiaoHui [6] conducted a series of experiments in wind tunnel to investigate the motion of natural snow (fresh snow and old snow). The results indicated that the threshold velocity of fresh snow is less than that of old snow, and that the mass flux rates of the two kinds of snow

exponentially decrease with height whereas the snow transport rate increases exponentially with wind velocity.

Although various studies dealing with snowdrift simulation have been carried out in past years, more details need determination, e.g. the right length of correctly modelled terrain, complete with particles, upstream of the region of interest. It is set to be roughly one or more meters by experience.

## 2 Outdoor experiments of snowdrifts around cuboid models

### 2.1 Study site and setup

The outdoor experiments were undertaken on an open flat base within the campus of Harbin Institute of Technology (HIT), China (Figure 1). It was part of a large snow engineering program. The wind pattern in this area has been previously found to be constant during the winter (Figure 2). The uniform wind direction makes the area ideal for snowdrift measurements.

<figure><img src="image_234.png" /><figcaption>**Figure 1.** Locations where outdoor measurements were conducted</figcaption></figure>

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

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

<footnote n="2">²Key Lab of Structures Dynamic Behavior and Control of China Ministry of Education, Harbin Institute of Technology, Harbin 150090, China</footnote>
<figure><img src="image_236.png" /><figcaption>Figure 2. Wind frequency diagram in the winter of 2012-2013</figcaption></figure>

The models were of rectangular shape with dimensions 750×600×750mm, which were set in the north direction.

The models were installed on December 15 on flat bare ground (Figure 3). In the following 3 months, Snow density and height were measured. Wind speed and direction were measured at 1.8m on a PC-4 automatic weather station, and the sampling interval was 1 minute.

<figure><img src="image_237.png" /><figcaption>Figure 3. Model setup</figcaption></figure>

## 2.2 Results

The results on 2014/12/29 have been selected for this investigation. Figure 4 shows the snowdrift geometry measured on 29th December, 2014. The incidence angle of the wind was 22.5° on the windward wall. The mean wind speed in the dominant direction was approximately 0.7m/s (at 1.8m) and the maximum speed was 1.6m/s, as indicated in Figure 5.

<figure><img src="image_238.png" /><figcaption>Figure 4. Photo of snowdrift geometry on 2014/12/29</figcaption></figure>

<figure><img src="image_239.png" /><figcaption>Figure 5. Wind speed and rose map for the whole period on 2014/12/29 at Harbin</figcaption></figure>

The mean snow thickness on the flat ground was 30mm. Approximately 400 points of snow heights around models were used to interpolate a surface, as indicated in Figure 6.

<figure><img src="image_240.png" /><figcaption>Figure 6. Contour map on 2014/12/29 at Harbin</figcaption></figure>

# 3 Experimental reduced scale snowdrift simulations

Considering the uncertain details in wind tunnel experiments, a series of experiments were conducted for comparison with the field measurements.

## 3.1 Kind's simulation requirements

Kind's approach [1] is an essential theory in wind tunnel experiments. It consists of topography modeling, flow modeling, saltation modeling, drifting pattern modeling and time scales.

### 3.1.1 Geometric similarity

The model and prototype length scale should be determined first, according to the cross area of wind tunnel. The traditional limit is 5%.

$$ \frac{L_M}{L_P} = \text{constant} $$

Here, $L_M$ is the reference length of the model and $L_P$ is the reference length of the prototype. Unlike model scale, the length scale ratio of the model and prototype
snow particle size can be ignored, according to Kind's analysis [1].

### 3.1.2 Similarity of surface particle motion

If it is assumed that drag is the only aerodynamic force acting on a free-flying particle. The requirements for dynamic similarity of two-phase flow involving free-flying particles are thus shown to be:

$$ 
\begin{aligned}
(u_{\text{th}}^*/V)_M &= (u_{\text{th}}^*/V)_P \\
(V_{\text{TER}}/V)_M &= (V_{\text{TER}}/V)_P \\
(V^2/Lg)_M &= (V^2/Lg)_P
\end{aligned}
 $$

Here, $u_{\text{th}}^*$ is the threshold shear velocity, $V$ is the reference velocity of the flow, $V_{\text{TER}}$ is the terminal velocity of the particles, $L$ is the reference length, $g$ is the gravitational acceleration, subscript M denotes models and subscript P denotes prototypes.

### 3.1.3 Reynolds number limits

$u_{\text{th}}^*/2gV$ should be greater than 30 for both model and prototype, or the Reynolds number should be matched in the model and prototype cases.

### 3.1.4 Time scale

Time scaling of the snowdrift is an important similarity requirement in relating experimental results to prototype conditions. For mass accumulations the ratio of model to prototype time scales is given by the usual relationship, namely

$$ \left( \frac{t_M^*}{t_P^*} \right)_{\text{mass accumulation}} = \frac{L_M / V_M}{L_P / V_P} $$

<table><caption>**Table 1.** Experimental parameters based on the prototypes at Harbin</caption><thead><tr><th>NO.</th><th>Prototype</th><th>Particle</th><th>Scale</th><th>Length of upstream (m)</th><th>Thickness (mm)</th><th>Speed (m/s)</th><th>Duration</th></tr></thead><tbody><tr><td>1-1</td><td>Harbin</td><td>NaHCO<sub>3</sub></td><td>1/8</td><td>1</td><td>5</td><td>5.5</td><td>15min</td></tr><tr><td>1-1.5</td><td>Harbin</td><td>NaHCO<sub>3</sub></td><td>1/8</td><td>1.5</td><td>5</td><td>5.5</td><td>15min</td></tr><tr><td>1-2</td><td>Harbin</td><td>NaHCO<sub>3</sub></td><td>1/8</td><td>2</td><td>5</td><td>5.5</td><td>15min</td></tr></tbody></table>

But for volume accumulation the ratio is given by the modified relationship

$$ \left( \frac{t_M^*}{t_P^*} \right)_{\text{volume accumulation}} = \frac{\rho_P / \sigma_P}{\rho_M / \sigma_M} \cdot \frac{L_M / V_M}{L_P / V_P} $$

Here, $\rho$ is the density of fluid and $\sigma$ is the density of particles.

Also, the model and prototype flow fields must be kinematically similar. The particles used for simulation of snow drifting should lie within a fairly narrow size range. If possible, the angle of repose of the model and prototype materials should be equal.

## 3.2 Experimental parameters

A direct type wind tunnel, with a test section that was 400cm long, 50cm wide, and 50cm high, was used in this experiment. Figure 7 shows a schematic view of the wind tunnel and model location.

<figure><img src="image_241.png" /><figcaption>**Figure 7.** Schematic view of the wind tunnel measurements</figcaption></figure>

Since little guidance can at present be offered as to the required length of the upstream fetch, the length was set to be more than one meter as Kind suggested. The outdoor measurements, conducted by authors at Harbin China, was selected as the prototype data.

The experimental parameters were determined based on the prototype parameters discussed above, according to Kind's similarity criteria. Table 1 shows the experimental modeling parameters based on the prototypes at Harbin.
## 3.3 Comparison of model and prototype snowdrifts

Figure. 8 compares the snow depth distributions obtained from the experimental model and the Harbin prototype. The darkest part in the figure indicates deep snow coverage. Snow depth is normalized by the mean snow thickness on flat ground. The coordinates of distribution on the X and Y axes are also normalized by the size of the model for the same axis.

Although the observed snow accumulation shows some differences to the predicted particle build-up, it does seem that the model satisfactorily predicts the accumulation location and scale. The first zone is the windward deposition before the cuboid, the second zone is the deposition in the leeward recirculation zone behind the cuboid and the third zone is erosion windward at a distance from the cuboid.

Compared with the upstream fetch lengths of 1.0 and 2.0m, the results for 1.5m correspond in greater detail with the prototype, i.e. the deposition in zone 1 and the deposition pattern in zone 2. Besides, the partial distribution is predicted at corners lateral to the cuboid. Therefore, this suggests that a reasonable length of upstream fetch is close to 1.5m.

<figure><img src="image_242.png" /><figcaption>a) Contour map of Harbin prototype</figcaption></figure>

<figure><img src="image_243.png" /><figcaption>b) Models with 1.0m lengths of upstream fetch</figcaption></figure>

<figure><img src="image_244.png" /><figcaption>c) Models with 1.5m lengths of upstream fetch</figcaption></figure>

<figure><img src="image_245.png" /><figcaption>d) Models with 2.0m lengths of upstream fetch</figcaption></figure>

**Figure 8.** Contour map of Harbin prototype and models with 1.0, 1.5 and 2.0m lengths of upstream fetch

## 4 Conclusions

This study has investigated the characteristics of drifting snow in outdoor measurements and wind tunnel tests. The following results were obtained. If upstream fetch is used to supply particles, it is suggested that the required length of the upstream fetch is 1.5m and a slight change is needed according to the specific snowdrift pattern. If the main pattern is deposition, this length needs to be increased. If the main type is erosion, reasonably short upstream fetches will probably suffice in all cases.

## References

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

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

3. N. bouvet. J. Surv. Geophys. **5-6**, 711 (1995)

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

5. T. Okaze, A. Mochida, Y. Tominaga, M. Nemoto, T. Sato, Y. Sasaki, K. Ichinohe. J. Wind. Eng. Ind. Aerodyn. **104-106**, 532 (2012)

6. X. LÜ, H. Ning, T. Ding. J. Technol. Sc. **55(4)**, 927 (2012)
# A new method for predicting snowdrift on flat roofs

Luyang Kang¹, Xuanyi Zhou¹,ª and Ming Gu¹

¹State Key Lab of Disaster Reduction in Civil Engineering, Tongji University, Shanghai, China

**Abstract.** Currently Eulerian method is employed by most scholars for the research on snow drifting around a building or on building roofs through CFD techniques. However, there is no scholar taking account of angle of repose of snow on the numerical simulation of snow drifting on roofs. In this paper, a new method for numerical simulation of snow drifting on flat roofs is employed, in which repose angle of snow is taken into account. The one-way coupling assumption is made for the interaction of wind and snow. The erosion or accumulation on snow surface is obtained by the calculation of snow flux on snow surface. The results of snow transports rates on flat roofs obtained by numerical simulation are analyzed in detail. Fitting equations are given to reflect the influence of roof span and wind velocity on transport rates of flat roofs.

## 1 Introduction

In recent years, with the development of computer technology, more and more scholars begin to study snowdrift through numerical simulation method [1, 2]. Based on the difference for treating snow phase, numerical simulation of snowdrift can be divided into two categories: Lagrangian method and Eulerian method. The difference lies in that the former directly simulate the trajectory of snow particles, while the latter studies the movement of snow particles based on space. Through Lagrangian method, the movement mechanism of snow particles can be studied more intuitively from the micro perspective. However, due to the requirement to track a large number of snow particles, it is difficult to promote this method in engineering applications. Currently, Eulerian method is employed by most scholars for the research on snow drifting around a building or on building roofs through CFD techniques.

Different from snowdrift around buildings, when snow accumulates on roofs, it will hold a certain repose angle at the edge of roofs, which will affect the wind-snow motion on roofs. However, there is no scholar accounting for this factor on the numerical simulation of snow drifting on roofs. To solve this problem, a new method for numerical simulation of snow drifting on flat roofs is applied in this paper. During numerical simulation, repose angle of snow is taken into account in the initial pattern of snow cover and one-way coupling steady method is applied.

Due to the fact that there is few field data about the mechanism of snow drifting on building roofs, when determining drift snow loads on building roofs, empirical relations for snow transport rate on building roofs, which are referred to field measurement on the ground, are usually used [3]. This method cannot take account into the influence of complex flow field of wind around the roof on the snow transport. Whereas, numerical simulation could study snow transport on roofs in detail. Hence, the new numerical method is used herein to simulate snow transport on flat roofs. The results of snow transport rates obtained by numerical simulation are analyzed in detail. Fitting equations are given to reflect the influence of roof span and wind velocity on transport rates of flat roofs.

The simulated results are compared with field measurement data on the ground.

## 2 Method of numerical simulation

### 2.1 Numerical method

The turbulent air field is modelled by employing the Reynolds averaged Navier-Stokes (RANS) equations of motion, in which the realizable k-ε model is adopted. For simulation of snow phase, Eulerian method is applied. The one-way coupling assumption is made for the interaction of wind and snow, namely snow drifting occurs under the action of wind whereas the migration of snow does not affect the movement of air. The erosion or accumulation on snow surface is obtained by the calculation of snow flux on snow surface.

The equations that govern the snow phase and adopted in the paper were referenced from Tominaga [2],

$$ \frac{\partial(\rho_s f)}{\partial t} + \frac{\partial(\rho_s f u_j)}{\partial x_j} = \frac{\partial}{\partial x_j} \left( v_t \frac{\partial \rho_s f}{\partial x_j} \right) + \frac{\partial}{\partial x_j} \left( -w_f \rho_s f \right) \quad (1) $$

where $u_j$ is the wind velocity which is calculated through the equations governing the air phase, $v_t$ is the turbulent kinematic viscosity, $\rho_s$ is the density of snow, $f$ is the snow volume fraction, $\rho_s f$ is the drifting snow concentration, and $w_f$ is the settling velocity of snow particles.

The erosive or accumulative snow fluxes on the snow surface of boundary condition are calculated using the formulas [1, 4] below,

$$ q_{\text{dep}} = \phi w_f \frac{u_{s_t}^2 - u_s^2}{u_{s_t}^2}, \quad u_s < u_{s_t} \quad (2) $$

$$ q_{\text{ero}} = -A_{\text{ero}}(u_s^2 - u_{s_t}^2), \quad u_s > u_{s_t} $$

where $\Phi$ is the snow concentration, which can be obtained from the numerical results of snow field. $A_{\text{ero}}$ is a coefficient that represents the snow pack bonding

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

<footnote>ª Corresponding author:zhouxytj@tongji.edu.cn</footnote>
strength. A value of $A_{ero} = 7 \times 10^{-4}$ was employed for the present analysis [1, 4].

The total snow mass flux $q_{total}$ on the snow surface is expressed as,

$$q_{total} = q_{ero} + q_{dep} \quad (3)$$

## 2.2 Numerical model

For in most current load codes and standards snow load on a flat roof is one of important content about calculating snow loads, this paper choose flat roof as the research object. The span of the flat roof studied is 12 m - 108 m (L) and the height is 6 m (H). Fig.1 demonstrates the dimensions of the flat roof studied.

<figure><img src="image_247.png" /><figcaption>Figure 1. Dimensions of the flat roof</figcaption></figure>

The two-dimensional snow transport on the prototype flat roof is modelled in the study. The model for numerical simulation is full scale and the calculation domain is 16L×30H. Structured progressive mesh was employed as grid scheme, growth factor of grid was 1.1 and the total number of quadrilateral grid cell was 60-455 thousands. The snow was supposed to be paved on the roof initially, which means only snow drifting without snowfall has been simulated.

Physical properties of snow particles used in numerical simulation are shown in Table 1. This study does not measure the physical properties of natural snow particles, and their values are determined with related literature [2, 5].

<table><caption>Table 1. Physical properties of snow particles</caption><thead><tr><td>Property</td><td>Value</td></tr></thead><tbody><tr><td>Diameter d</td><td>0.15 mm</td></tr><tr><td>Density ρ</td><td>150 kg m<sup>-3</sup></td></tr><tr><td>Threshold friction velocity u<sub>*t</sub></td><td>0.15 m s<sup>-1</sup></td></tr><tr><td>Angle of repose θ</td><td>50°</td></tr><tr><td>Setting velocity w<sub>f</sub></td><td>0.20 m s<sup>-1</sup></td></tr></tbody></table>

Table 2 shows the some key simulation parameters for the calculation conditions. Five different wind velocities are studied, and the total simulated conditions are 25. The initial snow depth is supposed to be 100 mm. The settings of boundary conditions and wind profiles refer to Zhou [6]. Thus, they are not given here in detail.

<table><caption>Table 2. Calculation conditions</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Roof height H</td><td>6 m</td></tr><tr><td>Roof span L</td><td>12 m, 24 m, 48 m, 72 m, 108 m</td></tr><tr><td>Wind velocity at the roof height U<sub>H</sub></td><td>4 m/s, 6 m/s, 8 m/s, 10 m/s, 12 m/s</td></tr><tr><td>Snow depth</td><td>100 mm</td></tr></tbody></table>

## 2.3 Validation of the Numerical Method

<figure><img src="image_248.png" /><figcaption>Figure 2. Snow distribution on roof surface [6]</figcaption></figure>

Scaled model test was conducted by Hu [7] at TJ-1 wind tunnel of Tongji University for the snow drifting on flat roofs, in which high-density silica sand was used and the geometric scale ratio of the test was 1:25. Zhou [6] made a full comparison of the numerical results using this method with those from wind tunnel test and field measurement. For in this paper, only initial snow transport rate on flat roofs are studied, the comparison of redistribution of snow in phase I are demonstrated in Fig. 2. In the figure, dimensionless snow depth is the ratio of the end-time snow depth of each phase on roof surface to initial snow depth $S_0$ (prototype 50 cm; model 20 mm). The snow distribution on the roof surface obtained from the numerical simulation was similar to that of the wind tunnel test in Phase I. Snow on the entire roof surface was eroded, and the erosion in the middle region was relatively uniform. The erosion amount obtained from the numerical simulation was slightly larger than that from the wind tunnel test in the middle roof. Whereas, in the front and rear regions the erosion amount of numerical simulation was slightly smaller than that from the wind tunnel test.

# 3 Results of numerical simulation

Field measurement of snow transport rate is usually conducted on the ground. There is few data about the snow transport rate on roofs. Due to the complex flow field around the roof, snow drifting on a roof could be significantly different from that on the ground. Whereas numerical method is a good way to obtain these prototype
data. In this section, snow transport rates from numerical simulation are discussed in detail.

## 3.1 Flow field of wind velocity

Fig. 3 shows the streamline of wind velocity around a flat roof, which is representative for the streamline of wind is almost the same for different roof spans. Before the wind begins blowing, the snow deposited on the roof is supposed to be paved uniformly on the roof with a repose angle. Due to this repose angle of snow at the roof edge, significant flow separation cannot be observed around the roof. The separation strength has been weakened significantly. The direction of wind velocity near snow surface is approximately parallel to the snow surface. From Fig. 3, we could deduce that snow particles are blown off the roof along the snow surface under the action of wind. Test results of Hu [7] proved this phenomenon.

<figure><img src="image_249.png" /><figcaption>Figure 3. Flow field around the flat roof (L = 48 m)</figcaption></figure>

## 3.2 Snow transport rate on flat roofs

In current load standards and codes, exposure coefficient, which is the ratio of roof snow loads to ground snow loads, is commonly used to reflect the effects of wind or other environment factors. Combined with meteorological data, transport rate could be used to determine snow loads on flat roofs. Hence, transport rate is analyzed here in detail.

<figure><img src="image_250.png" /><figcaption>Figure 4. Transport rate of snow on a flat roof</figcaption></figure>

Snow transport rate for a flat roof is defined here as the total snow mass blown off the roof per unit time and unit width, as shown in Fig. 4. Specifically, snow transport rate on roof surface is the product of bulk particle density and the integration of snow depth change

per unit width per unit time along the entire roof span. The calculation formula is as follows,

$$Q = \rho_b \int_L \Delta h(x) dx \quad (4)$$

Fig. 5 shows the influence of roof span on the snow transport rate for flat roofs. As one might expect, snow transport rate Q increases with the increasing roof span. The increasing amplitude is large for small span roof, and the amplitude decreases as roof span becomes larger. Logarithmic law is the most suitable for the relation between snow transport rate Q and roof span L. Fitting equations using logarithmic relation at different wind velocities are also given in Fig. 5.

<figure><img src="image_251.png" /><figcaption>Figure 5. Q as a function of L</figcaption></figure>

We could find from Fig. 5 that without snowfall snow transport on flat roofs are limited to the upwind fetch distance. Different from the relation given in Fig. 5, the reduction factor for snow transport rate on the ground given by Takeuchi [8] appears to be proportional to the square root of the fetch distance $\sqrt{F/210}$. Whereas based on field measurement, Tabler [9] gave this reduction factor as $(1-0.14^{FT})$, where T is the minimum distance for fully developed snow transport.

<figure><img src="image_252.png" /><figcaption>Figure 6. Q as a function of $V_H$</figcaption></figure>

Fig. 6 shows the influence of wind velocity at the roof height $V_H$ on snow transport rate Q. Similar to the empirical formulae obtained from field measurement on the ground (Table 3), the relation between Q and $V_H$
meets the power law. In addition, transport rates are increasing functions of wind velocity at the roof height.

Fitting functions of power law are given in Fig. 6. For the roof profile affects the wind flow around the roof, the growth rate of $Q$ for flat roofs is lower than that on the ground. The power of the fitting functions is between 2.4 and 2.7 for flat roofs. Whereas, the power of the empirical formulae obtained by ground field measurement, which are shown in Table. 3, is between 3 and 4. The order of $Q$ from numerical simulation is the same as the field data, which could be seen in Fig. 7.

<table><caption>**Table 3.** Empirical formulae for snow transport rate $Q$ (gm⁻¹s⁻¹), [8, 9]</caption><thead><tr><td>Investigator</td><td>Formula</td></tr></thead><tbody><tr><td>Khrgian</td><td>$Q = -5.8 + 0.267U_1 + 0.123U_1^2$</td></tr><tr><td>Inanov</td><td>$Q = 0.0295U_1^3$</td></tr><tr><td>Mel'nik</td><td>$Q = 0.092U_1^3$</td></tr><tr><td>Dyunin</td><td>$Q = 0.0334(1 - (4/U_1))U_1^3$</td></tr><tr><td>Komarov</td><td>$Q = 0.011U_1^{3.5} - 0.67$</td></tr><tr><td>Budd et al.</td><td>$\log Q = 1.22 + 0.0859U_1$</td></tr><tr><td>Kobayashi et al.</td><td>$Q = 0.03U_1^3$</td></tr><tr><td>Takeuchi</td><td>$Q = 0.0029U_1^{4.16}$ (dry snow)</td></tr><tr><td>Tabler</td><td>$Q = 0.00427U_{10}^{3.8}$</td></tr></tbody></table>

According to the square-root relationship used by O'Rourke [3], when roof span is larger than 210 m, the reduction factor would be larger than 1.0. Whereas, Tabler [9] used a better model to reflect the increase of transport rate with fetch distance, where the upper limit for the reduction factor is 1.0. Thus, the following equation used by Tabler is applied to compare the simulated results.

$$Q(U_{10}) = 4.27 \times 10^{-6} U_{10}^{3.8} (1 - 0.14^{L/210}) \quad (5)$$

Here, the unit of $Q$ is kgm⁻¹s⁻¹. For fresh snow, the minimum fetch distance for fully developed snow transport is set as 210 m, which was also used by O'Rourke [3]. Wind velocity at 10 m height is supposed to be the same as that at the roof height $V_H$.

The comparison of calculated results is shown in Fig. 7. To make it clear, only results of three different spans are shown. Comparison results for the other two spans have the same trend. At the lower wind velocity, the snow transport rate $Q$ on the ground is lower than that for roof. At a certain wind velocity, snow transport rate achieves the same value between that on the roof and that on the ground. This wind velocity is approximately 8 m/s for $L = 108$ m, 9 m/s for $L = 48$ m, and 10.5 m/s for $L = 12$ m. Here, we call this velocity as equilibrium velocity. When $V_H$ is larger than the equilibrium velocity, due to a larger power of the empirical formula, snow transport

rate on the ground is larger than that on the roof. At 12 m/s, compared with $Q$ on the roof, snow transport rate on the ground is 18.7%, 52.3%, 53.8% respectively larger for the roof span 12 m, 48 m and 108m. This means that when the empirical formula on the ground is used to predict snow loads on flat roofs, a lower wind exposure coefficient would be obtained, which is unfavourable for roof structures.

<figure><img src="image_253.png" /><figcaption>**Figure 7.** Comparison of snow transport rate $Q$</figcaption></figure>

## 4 Conclusions

In this paper, a new simulation method of snow drifting on flat roofs, where repose angle of snow is taken into account, is applied to calculate snow transport rates on flat roofs. The simulated results are compared with field measurement data on the ground. The main conclusions are as follows:

(1) The initial pattern of snow cover on roof surface has great influence on the wind-snow motion on roof surface. The existence of angle of repose of snow will alter the position of flow separation at the roof edge and weaken the strength of flow separation thus the snow redistribution on the roof is greatly affected. The direction of wind velocity near snow surface is approximately parallel to the snow surface.

(2) Snow transport rate $Q$ increases with the increasing roof span. Logarithmic law is the most suitable for the relation between snow transport rate $Q$ on the roof and roof span $L$. Whereas the relation between $Q$ and $V_H$ meets the power law.

(3) Compared with empirical formulae obtained on the ground, the power of the fitting equations for $Q$ on the roof is lower. At lower wind velocity, the snow transport rate $Q$ on the ground is lower than that for roof. When $V_H$ is larger than the equilibrium velocity, due to a larger power of the empirical formula, snow transport rate on the ground is larger than that on the roof.
References

1. J.H.M. Beyers, P.A. Sundsbø, T.M. Harms, Numerical simulation of three-dimensional, transient snow drifting around a cube, J. Wind Eng. Ind. Aerodyn., **92**(9): 725-747, (2004).

2. Y. Tominaga, T. Okaze, A. Mochida, CFD modeling of snowdrift around a building: An overview of models and evaluation of a new approach, J. Bui. Env., **4**:899-910, (2011).

3. O'Rourke, M., DeGaetano, A., Tokarczyk, J.D., Analytical Simulation of Snow Drift Loading. J. Struct. Eng. **131**, 660–667, (2005).

4. Naaim, M., Naaim-bouvet, F., Martinez, H., Numerical simulation of drifting snow: erosion and deposition rn.odels. Ann. Glaciol. **26**, 191–196, (1998).

5. D. Kuroiwa, Y. Mizuno, M. Takeuchi, Micromeritical properties of snow, *Proc. of Intl. Con. on Low Temp. Sci.: Phys. of Snow and Ice*, **1**(2): 751-772, (1967).

6. Zhou, X., Kang, L., Gu, M., Qiu, L., Hu, J., Numerical simulation and wind tunnel test for redistribution of snow on a flat roof. J. Wind Eng. Ind. Aerodyn. **153**, 92–105, (2016).

7. J. Hu. *Experimental study on wind-induced snow motion on building roof surface*, Master Thesis, Tongji University, (2013).

8. Takeuchi, M., Vertical profile and horizontal increase of drift-snow transport. J. Glaciol. **26**, 481–492, (1980).

9. Tabler R D., *Controlling blowing and drifting snow with snow fences and road design*. National Cooperative Highway Research Program Project, NCHRP Project 20-7(147), (2003).
# Improved design relations for roof snow drifts

Michael O'Rourke<fnref n="1" /> and John Cocca<fnref n="2" />

**Abstract.** The paper will present new relations for roof snow drifts which include a measure of winter windiness. The relation was developed by determining the annual maximum drift loads at 42 locations across the US through physics based numerical simulation using 19 years of weather records for both snowfall and wind speed at each location. Annual maximum were then used to determine the 50 year MRI drift for various upwind fetch distances. The new building code suitable design drift relation involving a map-able wind parameter, 50 year ground snow load and upwind fetch was then developed using multiple regression analysis.

## 1. Introduction

Structural loads due to drifted snow continue to be one of the key considerations for roof design in North America. As an example, a FEMA building collapse investigation team surveyed damage to school roofs in the Greater Boston region resulting from record snow in February 2015. For all four schools visited, collapse was due to localized roof snow drifts.

In North American building codes, snow drift loads were originally taken to be a multiple of the ground snow load. Current relations for the height of the roof snow drift are based upon an analysis of insurance company case histories, in which the height is presented as a function of the ground snow load and the fetch (i.e. length) of the upwind snow source area. Absent from the current relations is any measure of the local wind speed during the winter.

## 2. Current U.S. roof step drift relations

Current U.S. building codes utilize a roof snow drift relation based upon a statistical analysis and multiple regression analysis of empirical case history information. The location, ground snow load and roof geometry were the typical available case history parameters. However wind speed and direction for the case history drifts was generally not available. Hence the current U.S. empirical relation bases drift size on ground snow load (50 year MRI value in ASCE 7) and upwind fetch.

Specifically for a roof step in Figure 1 with wind from left to right, the height of the nominally triangular roof snow drift above the balanced snow, h_d, is given by:

$$h_d = 0.43(L_u)^{\frac{1}{3}}(P_g + 10)^{\frac{1}{4}} - 1.5 \quad (1)$$

where h_d is in feet, the upwind fetch, L_u, is in feet and the 50 year ground snow load, P_g, is in lbs/ft<fnref n="2" />. As one might expect, the drift height is an increasing function of the upwind snow source region characterized by the fetch and ground snow load. However the increase is not linear. If the upwind fetch or the ground snow load were to double, the drift height would increase but not double.

<figure><img src="image_254.png" /><figcaption>Figure 1. Roof Snowdrift at Leeward Roof Step</figcaption></figure>

The horizontal extent of the roof step drift is taken as four times its height. The 1:4 rise to run is thought to be related to the angle of repose for drifted snow. Finally the snow density γ (lbs/ft³) needed to convert drift height into drift load is given

$$\gamma = 0.13p_g + 14 \quad (2)$$

The snow density increase with ground snow load is due to compaction (densification) near the bottom of the snow pack, due to the weight of the snow near the top.

The relations in equation (1) and (2) have been in use in the U.S. for the past 20 years or more. Although intuitively appealing, they do not include wind which theoretically should influence roof snow drift size.

It can be argued that improvements in drift load estimation procedures in building codes and load standards require an ability to numerically simulate roof snow drifts.

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

<footnote n="1">¹Rensselaer Polytechnic Institute, New York, USA</footnote>

<footnote n="2">²Wiss, Janney, Elstner Associates, Inc., Connecticut, USA</footnote>
### 3. Simulation of roof snow drifts

Fortunately Tabler (1994) provides guidance in his pioneering work related to snow fences at highway cuts. Tabler considered a two-step physics based process wherein the snow flux or transport rate of snow is first established and then a trapping efficiency is applied to estimate drift size. As one might expect, Tabler's transport rate for an unlimited fetch $T_r(V)$ (lbs/foot width perpendicular to the wind direction/hr) for ground snow is a function of wind speed

$$T_r(V) = (4.8 * 10^{-4})V^{3.8} \quad (3)$$

where V is the wind speed in miles/hr at a height of 32.8 ft (10 m). Based upon full scale measurements by Takeuchi (1980), an upwind fetch of roughly 750 ft (225 m) is needed for full (unlimited fetch) transport. For shorter fetch distances, a correction factor is utilized

$$T_r(V) * (L/750)^{0.5} \quad (4)$$

where L is the upwind fetch in feet.

The wind flow patterns and trapping efficiency for a snow fence upwind of a highway cut are much different than those for a roof step in Figure 1. Based upon water flume tests, O'Rourke et al. (2005) have shown that the trapping efficiencies at a leeward roof step ranged between 47 and 58%, and a 50% trapping efficiency is used herein. That is, for a leeward roof step the initial growth rate (prior to complete filling of the aerodynamic shade region) is half the transport rate in equation (2) for unlimited fetch distances or equation (3) for fetch distances less than 750 ft.

Cocca (2006) used these procedures to determine the maximum annual leeward roof step drift for the years 1977 through 1996 (19 years) at 46 locations across the U.S. The 19 samples of the maximum annual drift load were fit with a Gumbel Distribution to determine the 50-year MRI drift load again at 46 locations across the U.S. The drift load was characterized by a dimensionless quantity, the drift ratio, DR. DR was defined as the 50-year drift load (half the drift height times drift length times snow density) divided by the product of the 50-year ground snow load, $P_g$ and the roof upwind fetch, $L_u$. Since the trapping efficiency is taken as 50%, the maximum possible value for DR is 0.5. This would correspond to a winter with a 50-year ground snow load and enough wind to blow all the snow on the upper level roof towards the roof step.

As one might expect, the 50 year drift ratio is a decreasing function of the upwind fetch and the 50-year ground snow load. For given amounts of ground snow and winter wind, one obtains smaller drift ratios for increasing upwind fetch distances. Similarly, for given upwind fetch and winter winds, one obtains smaller drift

ratios for increasing ground snow loads. In both cases, if the snow source area (proportional to product of fetch and ground snow load) is increasing, a fixed amount of winter wind will be able to transport only a smaller percentage of the upwind snow source.

Somewhat surprisingly, the DR also varied significantly with location. Figure 2 shows the 50-year drift ratios for an upwind fetch of 250 feet. Notice that areas with particularly large drift ratios are the upper Midwest (the Dakotas and Nebraska) and leeward of Lake Erie (Erie PA and Buffalo NY). Areas with particularly small drift ratios are located in the Pacific Northwest (Helena MT, Landers WY, Spokane WA and Yakima WA). For the five locations in the U.S. with a 50 year ground snow load in the 30 to 25 psf range, the drift ratios, again for a 250 foot upwind fetch, varied from 0.014 to 0.205, more than an order of magnitude difference. Clearly there is significant variation of winter windiness across the U.S. Note that the current ASCE 7 (ASCE, 2010) drift relations in equations (1) and (2) are functions of the upwind fetch and ground snow load, but not local winter windiness. As such they underestimate the snow drift hazard for some locations (upper Midwest) while overestimating it for other locations (Pacific Northwest).

<figure><img src="image_256.png" /><figcaption>Figure 2. Simulated 50-year drift ratios for upwind fetch $L_u = 250$ ft.</figcaption></figure>

As shown above, the site location in general and the local winter windiness in particular, have a significant influence on the magnitude of roof snow drifts. In order to incorporate this influence into a relationship suitable for use in building codes and load standards, a number of parameters were tested to determine a mapable value to characterize “Winter Windiness”. Specifically three parameters were considered. The first, W1, was the mean wind speed over the winter months. The second was W2, the percent time when the wind speed was greater than the snow transport threshold of 10 MPH. The third, W3, was the mean value of the wind speed to the 3.8 power when the wind speed was greater than 10 MPH. It turned out that W2 provided the best fit to the physics based simulated 50-year drift ratio values. Figure 3 is a plot of the winter windiness parameter W2 for various regions in the U.S. For the 46 locations across the U.S., W2 ranged
from 0.167 to 0.725 with a mean of .484 and a standard deviation of 0.119. As one might expect, W2 was high for the Upper Midwest (.498 to .659) and lee of Lake Erie (0.626 to 0.659), while W2 was low for the Pacific Northwest (0.152 to 0.353).

<figure><img src="image_257.png" /><figcaption>Figure 3. U.S. map of Winter Windiness Parameter W2</figcaption></figure>

Multiple regression analysis resulted in the following building code suitable predictive formula for the 50-year drift ratio.

$$DR = \frac{4.89 W_g^{1.81}}{P_g^{2.54} L_u^{2.96}} \quad (5)$$

Note that the 50-year ground snowload, $P_g$ (psf) and the upwind fetch, $L_u$ (ft.) are both in the denominator, that is DR is a decreasing function of both $P_g$ and $L_u$. However since the drift size (lbs/ft) is the product of DR, $P_g$ and $L_u$, drift size is an increasing function of $L_u$ and $P_g$ which is generally consistent with Equation (1) for the drift height, $h_d$.

## 4. Comparison between physics based and code suitable relations

Figure 4 shows a comparison of DR from the physics based numerical simulation and the code suitable value from Equation (5). The origin constrained best fit straight line thru the data points has a slope of 1.05: 1.0 and a correction coefficient $R^2$ of 0.43.

<figure><img src="image_258.png" /><figcaption>Figure 4. DR Simulated vs. DR Predicted Using W2</figcaption></figure>

Tables 1 thru 3 present the minimum, maximum, mean and standard deviation for the 50 year drift ratio from the physics based numerical simulation $DR_{sim}$, drift ratio from the code suitable predictive equation $DR_{pred}$, and from the current ASCE 7 procedure $DR_{ASCE}$, for upwind fetch distances of 100, 250 and 500 feet. Also presented are the ratios of Pred/Sim and ASCE/Sim.

Note that Figure 4 and Tables 1 thru 3 suggest that “on average” the code suitable predictive equation (Equation 5) provides a reasonable match to the expected physics based 50 year drift ratio. In Figure 4, the best fit straight line has a slope of 1.05:1.0. This is consistent with the mean values in Tables 1 thru 3 for the ratio Pred/Sim for upwind fetch distances of 100, 250, and 500 feet (1.04, 1.04, and 1.06 respectively). Also the range mean and standard deviation for $DR_{pred}$ (minimum value of 0.013, maximum value of 0.259 mean of 0.118 and standard deviation of 0.051 for $L_u$ = 250 feet) are similar to those for $DR_{sim}$ (0.013, 0.294, 0.125, and 0.066 again for $L_u$ = 250 feet).

However although “on average” $DR_{pred}$ nominally match $DR_{sim}$, Figure 4 shows significant differences in $DR_{pred}$ and $DR_{sim}$ values. As shown in Tables 1 thru 3, the ratio $DR_{pred}/DR_{sim}$ can be as low as roughly 0.5 and as high as roughly 2.0. That is, W2 is not the perfect parameter to characterize winter windiness. First of all, it is not the typical wind speed, simply the percentage of time the speed is above the threshold. Similarly, it does not provide a direct measure of the wind speed during and immediately after significant snowfalls. Rather it includes periods well before and after significant snowfall events. These “well before” and “well after” periods would have minimal impact on drift formation. Finally, drift size would be longer if the winter wind was typically out of the same direction, as opposed to out of multiple directions over the course of a given winter. Note in this regard that consistency of direction is not measured by W2.

For all its drawbacks, W2 is an easy to evaluate parameter and correlates reasonably well with ASCE 7 maps for the 3-second gust speed concurrent with freezing rain, which one uses to determine atmospheric icing loads in ASCE 7.

<table><caption>Table 1. 50 year Drift Ratios for $L_u$ = 100 ft.</caption><thead><tr><th></th><th>DR<sub>Sim</sub></th><th>DR<sub>Pred</sub></th><th>DR<sub>ASCE</sub></th><th><u>Pred</u><br><u>Sim</u></th><th><u>ASCE</u><br><u>Sim</u></th></tr></thead><tbody><tr><th>Min.</th><td>0.016</td><td>0.017</td><td>0.133</td><td>0.516</td><td>0.406</td></tr><tr><th>Max.</th><td>0.344</td><td>0.34</td><td>0.256</td><td>1.78</td><td>9.01</td></tr><tr><th>Mean</th><td>0.158</td><td>0.155</td><td>0.164</td><td>1.04</td><td>1.46</td></tr><tr><th>St. Dev.</th><td>0.075</td><td>0.068</td><td>0.03</td><td>0.316</td><td>1.47</td></tr></tbody></table>
<table><caption>**Table 2.** 50 year Drift Ratios for $L_u = 250$ ft.</caption><thead><tr><th></th><th>DR<sub>sim</sub></th><th>DR<sub>Pred</sub></th><th>DR<sub>ASCE</sub></th><th><u>Pred</u><br><u>Sim</u></th><th><u>ASCE</u><br><u>Sim</u></th></tr></thead><tbody><tr><th>Min.</th><td>0.013</td><td>0.013</td><td>0.114</td><td>0.460</td><td>0.427</td></tr><tr><th>Max.</th><td>0.294</td><td>0.259</td><td>0.250</td><td>2.06</td><td>10.576</td></tr><tr><th>Mean</th><td>0.125</td><td>0.118</td><td>0.151</td><td>1.04</td><td>1.82</td></tr><tr><th>St. Dev.</th><td>0.066</td><td>0.051</td><td>0.033</td><td>0.341</td><td>1.895</td></tr></tbody></table>

<table><caption>**Table 3.** 50 year Drift Ratios for $L_u = 500$ ft.</caption><thead><tr><th></th><th>DR<sub>sim</sub></th><th>DR<sub>Pred</sub></th><th>DR<sub>ASCE</sub></th><th><u>Pred</u><br><u>Sim</u></th><th><u>ASCE</u><br><u>Sim</u></th></tr></thead><tbody><tr><th>Min.</th><td>0.012</td><td>0.011</td><td>0.099</td><td>0.409</td><td>0.442</td></tr><tr><th>Max.</th><td>0.270</td><td>0.211</td><td>0.230</td><td>2.12</td><td>10.37</td></tr><tr><th>Mean</th><td>0.103</td><td>0.096</td><td>0.136</td><td>1.06</td><td>2.05</td></tr><tr><th>St. Dev.</th><td>0.06</td><td>0.042</td><td>0.031</td><td>0.388</td><td>2.089</td></tr></tbody></table>

## 5. Comparison between simulated and ASCE

Figure 5 shows a comparison of DR from the physics based numerical simulations and the corresponding drift ratios from the current ASCE 7 relations in Equation (1) and (2). For each of the 46 locations, drift ratios are determined for the three values of $L_u$; 100, 250 and 500 feet.

The most striking characteristic of Figure 5 is the relatively wide range of physics based $DR_{sim}$ values in comparison to the relatively narrow range for $DR_{ASCE}$. The $DR_{sim}$ values range from roughly 0.015 to roughly 0.30 while the $DR_{ASCE}$ values only range from roughly 0.10 to roughly 0.25. Hence for many sites the expected physics based DR is much lower than the DR from the current ASCE 7 provisions. Specifically for three fetch distances considered, $DR_{ASCE}$ was larger than the expected value (i.e. $DR_{sim}$), for about 75% of the sites (68% for $L_u = 100$ ft, 78% for $L_u = 250$ ft, and 83% for $L_u = 500$ ft).

The apparent overestimation associated with current ASCE drift relations is also evident by the ratios in Tables 1 thru 3. On average for the three fetch distances considered, the mean of the $DR_{ASCE}/DR_{sim}$ ratio was about 1.78 (1.46 for $L_u = 100$ ft, 1.82 for $L_u = 250$ ft, and 2.05 for $L_u = 500$ ft), with the overestimation increasing with increasing fetch distance.

Hence although W2 is an imperfect wind parameter, the drift size from the code suitable predictive equation utilizing W2 provides a much better estimation of the expected drift size from the physics based numerical simulation than the current ASCE relations which include no winter wind parameter.

<figure><img src="image_259.png" /><figcaption>**Figure 5:** DR ASCE vs. DR Simulated</figcaption></figure>

## 6. Comparison between code suitable and ASCE

Based upon the proceeding discussion, one could argue that drift relations based upon the code suitable predictive equation utilizing the winter wind parameter W2 (equation 5) are an improvement over the current ASCE 7 drift relations. The effects of this change would vary from site to site.

For example W2 for Lexington KY happens to be the mean value of 0.484, while Boston MA has the highest W2 (0.725) and Yakima WA has the lowest (0.167). For an upwind fetch of 250 feet, the 50 year drift ratios predicted by Equation (5) are 0.133, 0.217, and 0.013 for Lexington, Boston, and Yakima respectively. The 50 year drift ratios predicted by the ASCE procedure in Equation (1) and (2) for the same three cities are 0.190, 0.130 and 0.135 respectively. Hence if the current drift relations (equations 1 and 2) were replaced by the predictive relation (equation 5) the mean 50-year drift (lbs/ft) would decrease by roughly 30% (.133/.190 = .7). For particularly windy places such as Boston the 50-year drift (lbs/ft) would increase by about 70% (.217/.130 = 1.67), while for particularly calm locations such as Yakima, the 50 year drift ratio decreases by as much as 90% (.013/.135 = .096). That is, if the predictive relation in Equation 5 were used, drift loads for some locations would increase but for most locations they would decrease. Note that these values are based upon site specific values for both $P_g$ and W2. Due to the “rounding up” that typically occurs when parameters are mapped, the net effect for typical locations would be a decrease of something less than 30%.

## 7. Conclusions

Herein current roof step drifts from the ASCE 7 provisions are compared with physics based simulated and code suitable predicted roof step drifts. The ASCE 7 drift provisions were based upon an analysis of drift load case histories from insurance company files in which wind speed information was not available. The simulated drifts were based upon the physics of drift formation
(snow transport from the snow source area, and trapping
at the aerodynamic shade region), utilizing 19 years of
snowfall and wind speed information at 46 sites across
the US. Finally various easily mapped winter wind
parameters were investigated in order to determine a
predictive relation suitable for use in building codes and
load standards.

The code suitable predictive relationship was found to
provide a reasonable estimate to the expected drift size
based upon the physics of drift formation (i.e., simulated
drift). Comparison of the ASCE drifts to the simulated
drift suggest that the inclusion of a winter windiness
parameter is quite important. That is, expected drift from
numerical simulation showed much more variability
(both larger and smaller drift loads) that those from the
current ASCE 7 provisions.

Finally if the building code suitable predictive relation
were to replace the current ASCE 7 relations, one would
expect that on average design drift sizes (lbs/ft) would be
reduced by roughly 30%. For some locations the
reduction would be about 90%. At location which
experience much stronger than average winter wind the
design drift would increase by about 70%.

References

1. ASCE (American Society of Civil Engineers), 2010, Minimum design loads for building and other structures ASCE/SEI 7-10, ASCE, Reston VA, USA.

2. Cocca, J. (2006) Mathematical Simulation of 50-Year Snow Drift Ratios. MS Thesis, Rensselaer Polytechnic Institute.

3. M. O'Rourke, A. DeGaetano and J. D. Tokarczyk, 2005: Analytical Simulation of Snow Drift Loading, *Journal of Structural Engineering*. ASCE, Vol. 121. No. 4, April.

4. R. D. Tabler, 1994: Design Guidelines for the Control of Blowing and Drifting Snow, *Strategic Highway Research Program SHRP-H-381, National Research Council*, pp. 364.

5. M. Takeuchi 1980: *Vertical Profile and Horizontal Increase of Snow – Drift Transport*, J. Glaciology Vol. 26 pp. 481 – 492.
Session 5

# Structural loading II
# Specification of the design value of the ground snow load considering measurements of the snow height – part 2: regional approach

Michael Kasperski, Benjamin Czwikla<sup>a</sup>

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

**Abstract.** A regional model for the snow density is required to exploit the large number of stations which only monitor the snow height. Starting point are simultaneous observations of the water equivalents and the snow heights. Since the ensembles are confined, the station-wise input data have an inherent uncertainty. Furthermore, the observation period may vary. The shorter the observation period, the larger will be the uncertainties. The paper explains how to study the uncertainties based on simulations and presents a new strategy for considering the uncertainties in the identification of the model parameters. As example of application, the correlation of snow loads is studied for 14 stations in Mecklenburg-Vorpommern, Germany.

## 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 $h_s$ to snow loads s:

$$s = \rho_s \cdot g \cdot h_s \qquad (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 10 larger than the number of stations measuring the water equivalent (fig.1). Thus, interpolation between stations becomes easier, especially if additional parameters are to be considered as e.g. the altitude. Exploiting this huge amount of data requires a regional model of the snow density. This paper discusses the strategies to develop a regional model. Starting point is the approach for single stations as described in [1].

A naïve approach simply uses the identified parameters of the single stations as input for the geostatistics. Then, only data sets of the same length and quality should be used, thus avoiding biasing effects from different statistical stability of the input data. However, for the snow data in a specific region, the length of the observation periods might vary considerably. Then, it is required to consider the different level of uncertainties in the input data.

This paper develops a basic strategy for a refined approach. Chapter 2 explains how uncertainties can be studied based on simulations. Basic tool is the NORTA-approach which allows generating correlated data for any combination of probability distributions. In chapter 3, the refined approach is illustrated based on a simple example of a linear law. Chapter 4 presents results for a spatial correlation model of the snow loads.

<figure><img src="image_260.png" /><figcaption>**Figure 1:** Distributions of stations measuring the water equivalent and stations measuring only the snow height. (Germany, 2006)</figcaption></figure>

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

<footnote><sup>a</sup> Corresponding author: benjamin.czwikla@rub.de</footnote>
## 2 Stastical stability of the identified parameters of the snow density model

The developed snow density model [1] assumes that both the water equivalents and the snow heights follow an extreme value distribution type I. The snow density is described by a two parameter model as follows [1]:

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

The limiting density $\rho_{\infty}$ is obtained as the ratio of the standard deviations of the snow load and the snow heights:

$$ \rho_{\infty} = \frac{\sigma_s}{\sigma_h \cdot g} \qquad (3) $$

with g – gravitational acceleration ≈ 9.81 m/s²

The term $\Delta\rho$ is obtained as the weighted difference of the two mean values:

$$ \Delta\rho = \frac{\rho_{\infty} \cdot m_h - m_s / g}{h_{\text{ref}}} \qquad (4) $$

The reference height, which is required for consistency in the dimension, is an arbitrary value and can be assumed as $h_{\text{ref}} = 1$ m.

For confined ensembles, the statistical uncertainties in the estimated mean value and standard deviation are discussed in [1]. As main conclusion it is stated that even for observation periods beyond 100 years, the estimated parameters might deviate from the true parameters. Since the parameters of the snow density model are obtained as combinations of uncertain input values, the uncertainties in the final results will be larger than in the input values.

A consistent analysis of the uncertainties can be based on simulations. Since the snow loads and the snow heights are correlated, the simulations have to be based on the NORTA-strategy [2], which

- first generates correlated variables $x_n$ and $y_n$ following the standardized normal distribution,

- then transfers these values to correlated uniformly distributed variables u and v with $u = \Phi(x_n)$ and $v = \Phi(y_n)$, where $\Phi$ is the cumulative probability distribution for the standardized normal distribution, using e.g. the error function

- and finally transfers u and v into two correlated variables $t_1$ and $t_2$ which may follow any target probability distributions (NORTA = Normal To Anything)

The generation of correlated variables $x_n$ and $y_n$ is based on two uncorrelated standardized normally distributed variables $x_n$ and $z_n$ using the following transformation:

$$ y_n = r \cdot x_n + \sqrt{1-r^2} \cdot z_n \qquad (5) $$

with r – target correlation between $x_n$ and $y_n$

The transformation into two correlated uniformly distributed variables u and v leads to a distortion in the correlation, i.e. the obtained correlation between u and v differs from the input correlation of $x_n$ and $y_n$ [3]:

$$ r_{uv} = \frac{6}{\pi} \arcsin \left( \frac{r_{xy}}{2} \right) \qquad (6) $$

While in the first transformation step from normally distributed variables to uniformly distributed variables the whole solution space is kept with values for the correlation $r_{uv}$ from -1 to + 1, the further transformation to the final target distributions may lead to a limited solution space.

The final joint probability density of $t_1$ and $t_2$ is obtained from the joint normal probability density $f_{xy}$ as follows:

$$ f_{t_1 t_2}(t_1, t_2) = \frac{f_{xy}(x_n(u), y_n(v), r_{in})}{f_x(x_n(u)) \cdot f_y(y_n(v))} \cdot \frac{\partial u}{\partial t_1} \cdot \frac{\partial v}{\partial t_2} \qquad (7) $$

where the joint normal probability density is given as:

$$ f_{xy}(x_n, y_n) = \frac{\exp\left(-\frac{x_n^2 + 2 \cdot x_n \cdot y_n + y_n^2}{2 \cdot (1 - \rho_{in}^2)}\right)}{2 \cdot \pi \cdot \sqrt{1 - \rho_{in}^2}} \qquad (8) $$

The probability densities $f_x$ and $f_y$ correspond to the normal probability density:

$$ f_x(x_n) = \frac{1}{\sqrt{2 \cdot \pi}} \exp\left(-\frac{x_n^2}{2}\right) \qquad (9) $$

The terms $\partial u / \partial t_1$ and the $\partial v / \partial t_2$ are the partial derivatives of the respective target cumulative probability distributions. Finally, the correlation $r_{out}$ results from a double integration over $t_1$ and $t_2$ supposing that the two integration variables are used in normalized form, i.e. have a zero mean and a standard deviation of 1:

$$ r_{\text{out}} = \iint_{t_1, t_2} f_{t_1 t_2}(t_1, t_2) dt_1 dt_2 \qquad (10) $$

For $t_1$ and $t_2$ following the extreme value distribution type I, the range of possible correlations is limited on the negative side with $r_{min} = -0.8858$. The required input correlation $r_{in}$ for the generation of correlated normally distributed variables in equation (5) can be approximated by a polynomial of 3rd order as follows:

$$ r_{in} = \frac{r_{out}^3}{209.6} - \frac{r_{out}^2}{14.55} + \frac{r_{out}}{0.9397} \qquad (11) $$

with $r_{out}$ – target correlation between $t_1$ and $t_2$

Figure 2 presents examples of joint-probability densities for varying correlations of two variables following the type I extreme value distribution. The contour lines are oriented at the major diagonal and clearly show the skewness to the right tails.
<figure><img src="image_262.png" /><figcaption>**Figure 2:** Examples of the joint-probability density of two correlated random variables for different correlations</figcaption></figure>

For confined ensembles, the observed correlation may deviate from the true correlation. This is shown in figure 3 for two correlated variables following the extreme value distribution type I varying the ensemble size K and the true correlation r.

For the 225 analysed stations, the range of observed correlations is shown in figure 4. The majority of stations show large positive correlations; the average observed correlation is 0.82. It is reasonable to assume that the larger the correlation the higher will be the statistical stability of the density model.

The analysis of the statistical stability of the snow density model is more or less based on the results for Aachen. The correlation is set to 0.95, the threshold value for the snow loads is 0.05 kN/m², the scale is 0.15 kN/m². The corresponding values for the snow heights are 0.01 m for the threshold value and 0.1 m for the scale. The average number of events is set to 1.5, the observation period is 56 years. The snow density model then becomes:

$$ \rho_s = 152.91 - \frac{3.568 \cdot h_{\text{ref}}}{h_s} \quad [\text{kg/m}^3] \qquad (12) $$

with $h_{\text{ref}} = 1$ m

Figure 5 shows the scatter in the identified parameters of the snow density model. The two parameters show a strong negative correlation, i.e. for larger values of $\rho_\infty$ the parameter $\Delta\rho$ decreases. Figure 6 presents the cumulative probability distributions of the estimated snow density parameters. For both parameters, the randomly estimated values follow approximately a normal distribution.

The parameter $\rho_\infty$ has a variation coefficient of about 7%, the 90% confidence interval leads to values from 137.4 kg/m³ to 164 kg/m³, which correspond to a relative range from 0.9 to 1.07 around the true value of 1. The variation coefficient of the randomly estimated parameter $\Delta\rho$ is about 25%. The 90% confidence interval ranges from 2.42 kg/m³ to 4.67 kg/m³. Normalisation with the true value of 3.57 kg/m³ leads to a relative range of 0.68 to 1.31. The statistical stability of $\Delta\rho$ is considerably smaller than the statistical stability of $\rho_\infty$.
<figure><img src="image_263.png" /></figure>

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

<figure><img src="image_265.png" /><figcaption>**Figure 3:** Range of randomly estimated correlations for two variables following the extreme value distribution type I, r – true correlation value, K – number of observations</figcaption></figure>

<figure><img src="image_266.png" /><figcaption>**Figure 4:** Range of observed correlations for snow load and snow height for 225 German stations</figcaption></figure>

<figure><img src="image_267.png" /><figcaption>**Figure 5:** Random scatter in the identified parameters of the snow density model assuming that the annual non-exceedance probabilities of the snow height and the snow load follow the extreme value distribution type I input values:</figcaption></figure>

snow loads: $x_{th}$ = 0.05 kN/m², $sc$ = 0.15 kN/m²;
snow heights: $x_{th}$ = 0.01 m, $sc$ = 0.1 m; correlation 0.95,
average number of events 1.5 per year, 56 years
observation period

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

<figure><img src="image_269.png" /><figcaption>**Figure 6:** Cumulative probability distributions of the randomly estimated parameters of the snow density model</figcaption></figure>

Although in the above example the correlation between the snow height and the snow load has been assumed to be very large, the statistical stability of the estimated parameters of the snow density model is not convincing. Therefore, the station-wise identified parameters do not form a sound basis for the development of a regional model if the uncertainties are not considered.
### 3. Consideration of statistical uncertainties in the input data

A more reasonable approach uses for each individual station the probability density of the randomness in the uncertain parameters as a weighting factor. This strategy is illustrated in figure 7 for a sample of four points which are assumed to follow a linear law. The typical least square fit minimizes the sum of the squared differences between the theoretical model $y = a \cdot x + b$ and the observed values. Within the new strategy, for each support point a different statistical uncertainty is assumed, here represented by normal distributions with different standard deviations. The optimum fit is obtained as the straight line leading to the maximum sum of the probability densities for the four entries:

$$ \sum_{i} f_{y}(y_{i}) = \max \qquad (13) $$

In case of a linear approach and normal distributions for the uncertainties, the final equation becomes:

$$ \sum_{i} \frac{1}{\sigma_{i} \cdot \sqrt{2\pi}} \exp \left( -0.5 \cdot \left( \frac{ax_{i}+b-m_{i}}{\sigma_{i}} \right)^{2} \right) = \max \qquad (14) $$

with $m_i$ and $\sigma_i$ - mean value and standard deviation for point i

Table 1 summarizes the chosen input values for the example. These input values and the classical least square fit solution are shown in figure 7. The probabilistic optimum considering the demand from equation (14) leads to $y = 0.8951 \cdot x + 0.2246$.

<table><caption>**Table 1:** Uncertainties of the example input values</caption><tr><td>x</td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr><td>ymean</td><td>1.1</td><td>2.1</td><td>2.9</td><td>3.8</td></tr><tr><td>ysdev</td><td>0.1</td><td>0.15</td><td>0.12</td><td>0.09</td></tr></table>

<figure><img src="image_270.png" /><figcaption>**Figure 7:** Illustrating example of considering uncertainties for the identification of the theoretical model</figcaption></figure>

It is important to note, that the distribution of the uncertainties may differ from the normal distribution, e.g. the uncertainties may be skewed or limited. If no theoretical density is available, the uncertainties may be modelled based on simulations.

### 4. Example of application

For the development of a spatial model of the snow density it is vital to analyse the spatial correlation of the basic observation data for the snow heights and the snow loads. As example of application, the correlation structure of the snow loads in Mecklenburg-Vorpommern is studied. There are 14 stations with simultaneous observations of the water equivalent and the snow height, however with considerable different observation periods as shown in figure 8. The longest period covers 60 years, the shortest only 17 years. Therefore, uncertainties will vary considerably.

Basically, 14 stations lead to $(14^2-14)/2 = 91$ correlation coefficients. A minimum of 10 simultaneous observations is introduced for estimating the correlation for any pair of stations. Only observations are considered which occur at both stations in a time window of three days. These restrictions reduce the number of observed correlations to 64. As a first approach, the estimated correlations are plotted versus the distance between stations in figure 9.

<figure><img src="image_271.png" /><figcaption>**Figure 8:** Variation of the observation period in years for 14 stations in Mecklenburg-Vorpommern</figcaption></figure>

<figure><img src="image_272.png" /><figcaption>**Figure 9:** Observed correlations for 14 stations in Mecklenburg-Vorpommern for the water equivalent and the corresponding snow height</figcaption></figure>
Since the correlation is independent from the mean value and the standard deviation of the two variables, the simulation only needs to vary N and r<sub>true</sub>.

Basically, the uncertainties become smaller for increasing ensemble size. For an observed correlation of zero, the uncertain correlation is symmetric. If the absolute value of the observed correlation increases, the probability density of the uncertain input correlation becomes skewed and the scatter reduces.

The observed correlations in figure 9 suggest as basic simple law an exponential decay depending only on the distance between two stations:

$$
r(\Delta x) = \exp\left(-\frac{\Delta x}{L_{\text{ref}}}\right) \qquad (15)
$$

with Δx – distance between two station
L<sub>ref</sub> – reference value

There are two options for the classical approach of a least-square fit. The first step is to take the logarithm of both sides leading to:

$$
\ln(r(\Delta x)) = A \cdot \Delta x \qquad (20)
$$

The straight-forward approach uses the following sum:

$$
A = \frac{\sum_{i} \Delta x_{i}^{2}}{\sum_{i} \ln(r(\Delta x_{i})) \cdot \Delta x_{i}} \quad (21)
$$

This approach gives larger weight to small
observation values. Equal weight for all observations is
obtained if the summands in equation (21) are weighted
with the observation values:

$$
A = \frac{\sum_{i} r(\Delta x_i) \cdot \Delta x_i^2}{\sum_{i} r(\Delta x_i) \cdot \ln(r(\Delta x_i)) \cdot \Delta x_i} \quad (22)
$$

This approach, however, only works if all observed
correlations are larger or equal zero. Some of the
observed correlations in figure 9 do not meet this
demand. If the number of negative correlation entries
remains small, a less strict approach simply ignores these
values. A consistent approach is obtained if the sum of
the squared differences is minimized directly by varying
the reference value L<sub>ref</sub>:

$$
\sum_i \left( r(\Delta x_i) - \exp\left(-\frac{\Delta x_i}{L_{\text{ref}}}\right) \right)^2 \stackrel{!}{=} \min \qquad (23)
$$

Negative correlations randomly occur even for clearly
positive r<sub>true</sub> values if the ensemble is small. These
observations do not violate the assumed model which
only leads to positive values.

Figure 11 presents the final results in terms of the identified correlation model for the snow loads in Mecklenburg-Vorpommern. The best estimate of the reference length is 213.6 km. This result is obtained iteratively by solving equation (13) for 44 input values of Δx and 44 corresponding probability distributions of random correlations.

<figure><img src="image_273.png" /><figcaption>**Figure 11:** Identified model for the spatial correlation of the snow loads in Mecklenburg-Vorpommern considering uncertainties in the observed correlations</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. About 5500 stations measure only the snow
height. For codification of the design snow load, usually
maps are developed, i.e. a spatial model is required which
interpolates the obtained station-wise values.
Geostatistical tools allow interpolations considering
further influence parameters as e.g. the altitude. The
interpolation in a specific region obviously can be
improved if a snow density model is available. Then, also
the stations measuring only the snow height can be used
for the interpolation as well, i.e. the density of the input
grid is increased considerably.

A naïve approach simply uses the obtained station-
wise data. This strategy does not allow considering
different uncertainties in the input values from station to
station. A new strategy is developed which considers the
uncertainties by using the station-wise probability
densities as weighting coefficient for the fitting of the
theoretical model. The uncertainties can be analysed
based on simulations.

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. B. Czwikla, M. Kasperski: Specification of the design value of the ground snow load considering measurements of the snow height – part 1: single stations, ICSE 2016, Nantes

2. M.C. Cario, B.L. Nelson: Modelling and Generating Random Vectors with Arbitrary Marginal Distributions and Correlation Matrix, 1997, http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.48.281

3. S.T. Lee, J.L. Hammond: Generation of pseudo-random numbers with specified univariate distributions and correlation coefficients. IEEE Transactions on Systems, Man and Cybernetics, 557-561, 1975