# Snow engineering: recent advances

Philippe Delpech, Thomas K Thiis

► To cite this version:

Philippe Delpech, Thomas K Thiis. Snow engineering: recent advances. Delpech P.; Thiis T-K. ICSE 2016, 8th International Conference on Snow Engineering, Jun 2016, Nantes, France. CSTB - Centre scientifique et technique du bâtiment, 2016. hal-01426073v2

HAL Id: hal-01426073

https://cstb.hal.science/hal-01426073v2

Submitted on 18 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Recent advances

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

Philippe Delpech, Thomas K. Thiis

Editors
<figure><img src="image_4.png" /></figure>
## Organizers

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

Centre scientifique et technique du bâtiment, France

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

Norwegian University
of Life Sciences

Norwegian University of Life Sciences, Norway

## Conference chair

Philippe Delpech, Thomas K. Thiis

## Scientific committee

Perry Bartlet (Switzerland)

Andrzej Flaga (Poland)

Paolo Formichi (Italy)

Teruyuki Fukuhara (Japan)

Pascal Hagenmuller (France)

Peter Irwin (Canada)

Christian Jaedicke (Norway)

Alan Jones (Canada)

Michael Kasperski (Germany)

Gregorz Kimbar (Poland)

Alex Klein-Paste (Norway)

Holger Koss (Denmark)

Vivian Meylosund (Norway)

Florence Naiim (France)

Michael O'Rourke (U.S.A)

Sergey Pichugin (Russia)

Zoltan Sadowsky (Slovakia)

Toru Takahashi (Japan)

Xuanyi Zhou (China)

Jerzy Zuranski (Poland)
Preface

Since 1988 seven International Conferences on snow Engineering (ICSE) were held in Santa Barbara, USA (1988, 1992), Sendai, Japan (1996), Trondheim, Norway (2000), Davos, Switzerland (2004), Whistler, Canada (2008), and Fukui, Japan (2012). The ICSE, also known as “SnowEng”, is now an established forum for academics, researchers and professional snow engineers, who present, share and interpret the results of their work. Unlike other conferences, ICSE is dedicated to the actual applications of snow sciences to industrial areas, urban planning, road management, building, avalanche engineering and civil engineering.

The 8th edition of the ICSE, was held in Nantes 14-17 June 2016. In response to the call-for-paper we received about 60 abstracts coming from the entire world and illustrating the maintained interest of snow engineering issues in the international community.

The scientific committee selected 43 abstracts for possible oral presentations at the seminar and 10 for posters. In addition, 4 keynote lectures were invited to present overviews on selected topics.

These led to 51 papers published in these proceedings. They were grouped into the following topics:

- Snow drift,

- Structural loading,

- Transportation,

- Buildings / Simulation,

- Avalanche / Snow physics.

These proceedings include papers from 139 authors from 18 countries: Australia, Austria, Canada, China, Denmark, France, Germany, Italy, Japan, Netherland, Norway, Poland, Romania, Russia, Sweden, Switzerland, Ukraine and United States of America.

The papers provide an international state-of-the-art knowledge in the field of snow engineering. They show that the need of expertise in this domain is still high and we feel confident that the ICSE will provide a leading forum for information exchange between researchers and engineers, regulators and other stakeholders in the future.

We would like to thank the members of the scientific committee who did a wonderful work in reviewing the submitted abstracts, the members of the organizing committee, the authors and the sponsoring organizations for their contribution to the success of the conference.

Philippe Delpech, Thomas K. Thiis

Co-chairman
Table of contents

<table><tr><td>Preface</td><td>1</td></tr></table>

Session 1: Snow drift I

<table><tr><td>Keynote Lecture</td><td></td></tr><tr><td>Numerical simulation of snowdrift around buildings: past achievements and future perspectives</td><td>10</td></tr></table>

<table><tr><td>Yoshihide Tominaga</td><td>Niigata Institute of Technology, Department of Architecture and Building Engineering, Japan</td></tr></table>

<table><tr><td>Improvements of the viscous treatment of the snow phase in two-way coupled Eulerian-Eulerian simulations of drifting snow</td><td>20</td></tr></table>

<table><tr><td>Ziad Boutanios<sup>1,2</sup> Hrvoje Jasak<sup>2</sup></td><td></td></tr><tr><td><sup>1</sup> Binkz Incorporated, Canada</td><td></td></tr><tr><td><sup>2</sup> CFD Lab, FSB, University of Zagreb, Croatia</td><td></td></tr></table>

<table><tr><td>CFD prediction of snowdrift in a building array</td><td>26</td></tr></table>

<table><tr><td>Tsubasa Okaze<sup>1</sup>, Saeka Kato<sup>2</sup>, Yoshihide Tominaga<sup>3</sup>, Akashi Mochida<sup>4</sup></td><td></td></tr></table>

<table><tr><td><sup>1</sup> Tokyo Institute of Technology, Japan</td><td></td></tr></table>

<table><tr><td><sup>2</sup> Takenaka Corporation, Japan</td><td></td></tr></table>

<table><tr><td><sup>3</sup> Niigata Institute of Technology, Japan</td><td></td></tr></table>

<table><tr><td><sup>4</sup> Tohoku University, Japan</td><td></td></tr></table>

<table><tr><td>Outdoors experiments of snowdrift on typical cubes based on axial flow fan matrix in Harbin</td><td>30</td></tr></table>

<table><tr><td>Mengmeng Liu<sup>1,2</sup>, Qingwen Zhang<sup>1,2</sup>, Feng Fan<sup>1,2</sup></td><td></td></tr><tr><td><sup>1</sup> School of Civil Engineering, Harbin Institute of Technology, China</td><td></td></tr></table>

<table><tr><td><sup>2</sup> Key Lab of Structures Dynamic Behavior and Control of China Ministry of Education, Harbin Institute of Technology, China</td><td></td></tr></table>

<table><tr><td>Numerical simulation of snowdrift on a membrane roof and wind-induced response analysis under coupled wind and snow loads</td><td>34</td></tr></table>

<table><tr><td>Sun Xiaoying, He Rijin, Wu Yue</td><td>Key Lab of Structure Dynamic Behaviour and Control of the Ministry of Education, Harbin Institute of Technology, China</td></tr></table>

<table><tr><td>Developing experimental method for investigating snow deposition around buildings using snow substitutes</td><td>40</td></tr></table>

<table><tr><td>Jennifer Fiebig, Hans Holger Hundborg Koss</td><td>Technical University of Denmark (DTU), Dep. of Civil Engineering, Denmark</td></tr></table>

<table><tr><td>Use of numerical simulations of snow drift in planning of infrastructure – A case study from Northern Norway</td><td>47</td></tr></table>

<table><tr><td>Thomas K. Thiis<sup>1,2</sup></td><td></td></tr><tr><td><sup>1</sup> Norwegian University of Life Sciences, Norway</td><td></td></tr><tr><td><sup>2</sup> Multiconsult ASA, Norway</td><td></td></tr></table>
**Seismic responses of snowpack on a slope in a vibration experiment** 121

Yusuke Harada ¹, Wataru Takahashi ¹, Satoshi Omiya ¹, Hiroki Matsushita ², Takahiro Chiba ³, Masaru Matsuzawa ¹

¹ PWRI, Civil Engineering Research Institute for Cold Region, Japan

² PWRI, Snow Avalanche and Land Slide Research Center, Japan

³ Hokkaido University of Science, Japan

**Friction along a slider on snow** 130

Werner Nachbauer ¹, Sebastian Rohm ¹, Christoph Knoflach ¹, Joost van Putten ², Michael Hasler ²

¹ Department of Sports Science, University of Innsbruck, Austria

² Centre of Technology of Ski and Alpine Sports, University of Innsbruck, Austria

**Gliding friction of back country climbing skins** 132

Michael Hasler ¹, Sebastian Rohm ¹, Christoph Knoflach ¹, Joost van Putten ², Werner Nachbauer ¹, ²

¹ Centre of Technology of Ski and Alpine Sports, University of Innsbruck, Austria

² Department of Sports Science, University of Innsbruck, Austria

## Session 4: Snow Physics / Snow drift II

**Characterizing the snowpack stratigraphy and its mechanical stability with hardness profiles measured by the Avatech SP1** 135

Pascal Hagenmuller, Thibault Pilloix

Météo-France/CNRS, CNRM-GAME/CEN, France

**Improvement of requirements for modeling snowdrifts in wind tunnels based on the measurements at Harbin** 137

Qingwen Zhang ¹, ², Guolong Zhang ¹, ², Feng Fan ¹, ²

¹ School of Civil Engineering, Harbin Institute of Technology, China

² Key Lab of Structures Dynamic Behavior and Control of China Ministry of Education, Harbin Institute of Technology, China

**A new method for predicting snowdrift on flat roofs** 141

Luyang Kang, Xuanyi Zhou, Ming Gu

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

**Improved design relations for roof snow drifts** 146

Michael O'Rourke ¹, John Cocca ²

¹ Rensselaer Polytechnic Institute, USA

² Wiss, Janney, Elstner Associates, Inc., USA

## 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** 152

Michael Kasperski, Benjamin Czwikla

Ruhr-Universität Bochum, Research Team EKIB, Germany
**Collapse process of pipe-framed greenhouses under snow loading**
Kazuya Takahashi, Yasushi Uematsu
Department of Architecture and Building Science, Tohoku University, Japan
160

**Study on evaluation of roof snow load considering rain-on-snow surcharge: statistical evaluation of snow cover and precipitation in winter in Japan**
Masaya Otsuki¹, Toru Takahashi², Yoshihiko Saito¹, Takuya Tsutsumi³, Kikitsu Hitomitsu⁴
¹ Yukiken Snow Eaters Co., Ltd., Japan
² Department of Architecture, Chiba University, Japan
³ Northern Regional Building Research Institute, Japan
⁴ National Institute for Land and Infrastructure Management, Japan
166

**Structural damage caused by rain-on-snow load in Japan**
Toru Takahashi¹, Takahiro Chiba², Kazuki Nakamura³
¹ Department of Architecture, Chiba University, Japan
² Department of Architecture, Hokkaido University of Science, Japan
³ National Research Institute for Earth Science and Disaster Prevention, Japan
173

**A new method of predicting slide snow load for sloped roofs**
Xuanyi Zhou, Jialiang Li, Peng Huang, Ming Gu, Lulu Sun
State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, China
179

**Experimental study of the distribution of snow deposits on the surface of structures with complex three-dimensional shape of the roof**
Olga Poddaeva, Pavel Churin
Moscow State University Of Civil Engineering (Mgsu), Russian Federation
184

**Theoretical and experimental study of ice accretion due to freezing rain on an inclined cylinder**
Krzysztof Szilder
Aerospace, National Research Council, Canada
188

**Parametric approach for assessing risks due to falling ice and snow**
Jan Dale, Scott Gamble, Albert Brooks, Jill Bond
Rowan Williams Davies & Irwin Inc., Canada
198

## Session 6: Building / Simulation

**Keynote Lecture**
**Prediction of snow loads: past, present and future**
Peter A. Irwin
Rowan Williams Davies and Irwin Inc., Canada
206

**CFD simulation of drift snow loads for an isolated gable-roof building**
Yoshihide Tominaga¹, Tsubasa Okaze², Akashi Mochida³
¹ Niigata Institute of Technology, Japan
² Tokyo Institute of Technology, Japan
³ Tohoku University, Japan
214
**Analysis of snow drifts on arch roofs** 221

Michael O'Rourke ¹, Jan Potac ², Thomas K. Thiis ³

¹ Rensselaer Polytechnic Institute, USA

² Multiconsult ASA, Norway

³ Norwegian University of Life Sciences, Norway

**Falling snow and ice from buildings and structures: risk assessment and mitigation – two case studies** 228

Stefan Margreth

WSL Institute for Snow and Avalanche Research SLF, 7260 Davos Dorf, Switzerland

**Capture of windward drift snow** 234

Jan Potac ¹, Michael O'Rourke ², Thomas K. Thiis ³

¹ Multiconsult AS, Norway

² Rensselaer Polytechnic Institute, USA

³ Norwegian University of Life Sciences, Norway

**Wind tunnel tests and analysis of snow load distribution on three different large size stadium roofs** 238

Andrzej Flaga ¹, Łukasz Flaga ²

¹ Prof.D.Sc.Eng. Andrzej Flaga, Wind Engineering Laboratory, Faculty of Civil Engineering, Cracow University of Technology, Poland

² Ph.D.Eng. Arch. Łukasz Flaga, Faculty of Civil Engineering, Department of Technology of Building and Materials Processes, Częstochowa University of Technology, Poland

**Comparison of physical snow accumulation simulation techniques** 246

Albert Brooks, Scott Gamble, Jan Dale, Jill Bond

Rowan Williams Davies and Irwin (RWDI), Guelph, Ontario, CA

**Snowdrifts on two-level building roofs and modeling of snow density at Harbin** 256

Guolong Zhang ¹, ², Yu Zhang ¹, ², Feng Fan ¹, ²

¹ School of Civil Engineering, Harbin Institute of Technology, China

² Key Lab of Structures Dynamic Behavior and Control of China Ministry of Education, Harbin Institute of Technology, China

**Advantages and features of four different snow utilizing facilities** 260

Seiji Kamimura ¹, Yoshiomi Ito ², Junki Zen ³

¹ Nagaoka University of Technology, Department of Mechanical Engineering, Japan

² Yuki-daruma (snowman) Foundation, Japan

³ Nagaoka University of Technology, Graduate School of Engineering, Japan

## Session 7: Transport

**Modelling the thermal conductivity of melting snow layers on heated pavements** 269

Anne Nuijten ¹, Knut Vilhelm Høyland ¹,², Cor Kasbergen ³, Tom Scarpas ³

¹ NTNU, Department of Civil and Transport Engineering, Norway

² Sustainable Arctic Marine and Coastal Technology (SAMCoT), Centre for Research-based Innovation (CRI), Norwegian University of Science and Technology, Norway

³ TU Delft, Department of Structural Engineering, the Netherlands

**Snow engineering questions related to road and rail vehicles** 276

Jean-Paul Bouchet, Sylvain Aguinaga, Pierre Palier, Philippe Delpech

Centre Scientifique et Technique du Bâtiment (CSTB), CAPE Department, France
**Anti- and de-icing of walking and cycle paths – Field trials of new follow-up techniques for quantifying salt amount and resulting ice quality** 282

Göran Blomqvist, Bengt Lindström, Ida Järlskog, Emelie Karlsson, Anna Niska

Swedish National Road and Transport Research Institute (VTI), Sweden

**Shallow geothermal switch point heating system** 284

Lars Staudacher¹, Damian Schink², Dr. Roman Zorn³, Dr. Hagen Steger⁴

¹Bavarian Center for Applied Energy Research, Germany

²Pintsch Aben geotherm GmbH, Germany

³European Institute f. Energy Research (EIFER), Germany

⁴Karlsruher Institut für Technologie (KIT) Institut für Angewandte Geowissenschaften, Germany

**A method for estimating road friction coefficients with ice film subjected to melting by de-icing agents** 289

Akihiro Fujimoto¹, Shunsuke Tanaka¹, Kenji Sato¹, Roberto Tokunaga¹, Naoto Takahashi¹, Tateki Ishida¹,

Kiyoshi Takeichi²

¹PWRI, CERI, Japan

²Hokkai-Gakuen University, Japan

**Airplane braking friction on dry snow, wet snow or slush contaminated runways** 297

Alex Klein-Paste

NTNU, dept. of Civil and Transport Engineering, Winter Maintenance Research Group, Norway

**Engineered pavements of snow and ice** 302

Adrian McCallum, Greg White

University of the Sunshine Coast, Queensland, Australia

**Performance of remote road surface sensor on different pavement types** 312

Naoto Takahashi, Kenji Sato, Roberto Tokunaga

Civil Engineering Research Institute for Cold Region, Traffic Engineering Research Team, Japan

## Poster session

**The protection of roads from blizzards** 320

Tatiana Samodurovaa, Olga Gladysheva, Jurij Baklanov, Konstantin Panferov

Voronezh State University of Architecture and Civil Engineering, Russia

**Influence of sunshine hours in fine weather on the rate of wintry accidents** 325

Akira Saida, Masayuki Hirasawa, Naoto Takahashi, Tateki Ishida

Civil Engineering Research Institute for Cold Region, Public Works Research Institute, National Research and Development Agency, Japan

**Roof snow slide-off experiments using membrane deformation** 330

Hiroaki Terasaki, Teruyuki Fukuhara

University of Fukui, faculty of engineering, Japan

**Accuracy of snow depth measurements on roods measured with photogrammetry** 334

Takahiro Chiba¹, Thomas K. Thiis²

¹Department of Architecture, Faculty of Engineering, Hokkaido University of Science, Japan

²Department of Mathematical Science and Technology, Norwegian University of Life Science, Ås, Norway
A new ring-shaped wind tunnel facility to study wind-packing of snow 339

Christian G. Sommer ¹,², Michael Lehning ¹,², Charles Fierz ¹

¹ WSL Institute for Snow and Avalanche Research SLF, Switzerland

² CRYOS, School of Architecture, Civil and Environmental Engineering, EPFL, Switzerland

The use of sheet piles as measures against rapid mass flows 342

Árni Jónsson ¹, Guðmundur Heiðreksson ², Torfi B. Jóhannsson ³, Magnús Steinarsson ³

¹ Norwegian Geotechnical Institute (NGI), Norway

² Icelandic Road and Coastal Administration (IRCA) Iceland

³ MogT Engineering, Iceland

In flight wet snow particles characterization 344

Philippe Delpech ¹, Guy Febvre ², Christophe Gourbeyre ², Dominique Lenoir ¹, Fabrice De Oliveira ¹

¹ Centre Scientifique et Technique du Bâtiment, France

² Laboratoire de Météorologie physique, UMR 6016, France

CR 1-1-3/2012 - the snow loads code in Romania 351

Alexandru Aldea, Sorin Demetriu, Dan Lungu, Cristian Neagu, Radu Vacareanu, Cristian Arion

Technical University of Civil Engineering Bucharest, Romania

Authors index 359
Session 1

# Snow drift I
# Numerical simulation of snowdrift around buildings: past achievements and future perspectives

Yoshihide Tominaga¹,ᵃ

¹Niigata Institute of Technology, Department of Architecture and Building Engineering, 1719 Fujihashi, Kashiwazaki, Japan

**Abstract.** This paper reviews the current status of computational fluid dynamics (CFD) modeling of snowdrift around buildings and discusses the review findings to give insight into future applications. First, an overview of past achievements in CFD simulations of snowdrift around buildings is provided by reviewing previous research to clarify their applicability and limitations. Next, recent application examples are presented with respect to several topics. Finally, the remaining development issues and new challenges are presented from an engineering viewpoint. Clearly, CFD has considerable potential, as established in this review. However, to use it appropriately, sufficient attention should be paid to modeling theory and implementation. Because snowdrift is a highly complex phenomenon, we should not excessively rely on numerical simulation. It is important to use data from actual phenomena for all aspects of study.

## 1 Introduction

Under specific conditions of snowfall and wind, a large amount of snowdrift forms on and around buildings. This often leads to the collapse or cracking of buildings because of unbalanced snow loads and overhanging snow on roofs. Moreover, snowdrift on and around buildings not only is very troublesome to remove but also causes difficulties for inhabitants, pedestrians, and vehicular traffic. Therefore, to develop suitable measures for handling difficulties and hazards due to snow in built-up environments, it is necessary to predict the spatial distributions of snow coverage around buildings with high accuracy.

For predicting snowdrift around buildings, considerable experimental research has been carried out over the years using wind or water tunnels. The advantage of such methods is that conditions such as the wind speed, wind direction, and characteristics of snow particles (often modelled by artificial particles) are very controllable. Therefore, these methods have been recognized as useful tools for the practical design for snowdrift around buildings (e.g., [1-5]). However, wind/water tunnels are not always available and are usually expensive and time consuming to use. Furthermore, they have serious limitations with regard to the similarity law. On the other hand, in recent decades, dramatic increases in computational power together with advances in computer software have allowed engineers/researchers to more accurately simulate many types of problems across various engineering realms. Especially, the computational fluid dynamics (CFD) technique has become a strong tool for solving fluid-related engineering problems. CFD simulation has some particular advantages over experimental (full-scale and reduced-scale) testing; it can provide detailed information on the relevant flow variables in the entire calculation domain under well-controlled conditions and without similarity constraints [6]. As a result, the CFD-based numerical simulation technique is expected to be a powerful design tool for addressing the problem of snowdrift around buildings in snowy regions. Hereafter in this paper, “numerical simulation” is used to mean “CFD simulation” only.

This paper reviews the current status of CFD modeling of snowdrift around buildings and discusses the review findings to give insight into future applications. First, in section 2, an overview of past achievements in CFD simulations of snowdrift around buildings is provided by reviewing previous research to clarify their applicability and limitations. Next, in section 3, recent application examples are presented with respect to several topics. Finally, in Section 4, the remaining development issues and new challenges are presented from an engineering viewpoint.

Readers are forewarned that the papers were chosen for this review based on relevance and familiarity to the author, and it is likely that additional and perhaps better examples are available. The author apologizes to all researchers whose valuable contributions are not included in this paper.

## 2 Overview of past achievements

Figure 1 shows the year-wise chart of the number of papers returned by Google Scholar for searches using the keywords “Numerical simulation” and “Snowdrift”. It is clearly seen that the number of relevant papers that include these words has increased markedly in recent years.

<figure><img src="image_7.png" /><figcaption>**Figure 1.** Search results for “Numerical simulation” and “Snowdrift” on Google Scholar in March 2016.</figcaption></figure>

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

<footnote>ª Corresponding author: tominaga@abe.niit.ac.jp</footnote>
<table><caption>**Table 1.** Overview of CFD studies on snowdrift published previously in archival journals.</caption><thead><tr><th>Authors (year)</th><th>Configuration</th><th>Turbulence modeling</th><th>Suspension</th><th>Salvation</th><th>Validation</th></tr></thead><tbody><tr><td>Uematsu et al. (1991) [7]</td><td>3D</td><td>RANS (0-eq)</td><td>EA<sup>a</sup></td><td>MT [21]</td><td>Y (FO)</td></tr><tr><td>Liston et al. (1993) [8]</td><td>2D/fence</td><td>RANS (SKE)</td><td>NA</td><td>MT [22]</td><td>N</td></tr><tr><td>Bang et al. (1994) [9]</td><td>3D/buildings</td><td>RANS (SKE)</td><td>EA<sup>b</sup></td><td>NA</td><td>Y (FO)</td></tr><tr><td>Sundsbø (1998) [10]</td><td>2D/building step</td><td>RANS (1-eq)</td><td>EA<sup>c</sup></td><td>TE</td><td>N</td></tr><tr><td>Naaïm et al. (1998) [11]</td><td>2D/fence</td><td>RANS (MKE)</td><td>EA<sup>d</sup></td><td>EF [23]</td><td>Y (WT)</td></tr><tr><td>Sekine et al. (1999) [12]</td><td>2D/building</td><td>NA</td><td>LA</td><td>MT [21]</td><td>Y (WT)</td></tr><tr><td>Tominaga and Mochida (1999) [13]</td><td>3D/actual building</td><td>RANS (MKE)</td><td>EA<sup>a</sup></td><td>NA</td><td>N</td></tr><tr><td>Thiis (2000) [14]</td><td>3D/model buildings</td><td>RANS (SKE)</td><td>EA<sup>b</sup></td><td>NM</td><td>Y (FO [24])</td></tr><tr><td>Beyers et al. (2004) [15]</td><td>3D/isolated cube</td><td>RANS (SKE)</td><td>EA</td><td>TE</td><td>Y (FO [25])</td></tr><tr><td>Beyers et al. (2008) [16]</td><td>3D/buildings</td><td>RANS (SKE)</td><td>EA</td><td>MT [22]</td><td>Y (FO)</td></tr><tr><td>Tominaga et al. (2011) [17]</td><td>3D/isolated model building</td><td>RANS (MKE)</td><td>EA<sup>a</sup></td><td>EF [23]</td><td>Y (FO [26])</td></tr><tr><td>Tominaga et al. (2011) [18]</td><td>3D/isolated model building</td><td>RANS (MKE)</td><td>EA<sup>a</sup></td><td>EF [23]</td><td>Y (FO)</td></tr><tr><td>Thiis and Ferreira (2015) [19]</td><td>3D/array of pillars</td><td>RANS (MKE)</td><td>EA<sup>b</sup></td><td>NM</td><td>N</td></tr><tr><td>Okaze et al. (2015) [20]</td><td>2D/drifting snow boundary layer</td><td>RANS (MKE)</td><td>EA<sup>a</sup></td><td>EF [23]</td><td>Y (WT [27])</td></tr></tbody></table>

RANS: Reynold-averaged Navier-Stokes; 0-eq = zero-equation (constant eddy viscosity) model; 1-eq = one-equation (k) model; SKE = Standard k-ε model; MKE = Modified k-ε models; EA = Eulerian approach; LA = Lagrangian approach; MT = Mass transport rate model; EF = Erosion-flux model; TE = Transport equation for saltation region; Y = yes; N = no; WT = Wind tunnel experiment; FO = Field observation; NA = Not applicable; NM = Not mentioned

a Transport equation of drifting snow density

b Drift-flux model

c Transport equation of snow volume fraction

d Transport equation of particle concentration

## 2.1 Papers published in scientific journals

Table 1 lists previous CFD studies on snowdrift around obstacles published in archival journals. Several meteorological studies in regional scales that do not solve the flowfield directly were excluded based on an engineering point of view (e.g., [28-30]). Papers that do not provide enough information about the snowdrift calculation are also excluded from the list.

In the 1990s, applications of CFD to wind engineering problems increased (e.g., [6, 31]), and this technique was being applied to various environmental problems including snowdrift. A pioneering study on numerical simulations of snowdrift around structures by CFD was conducted by Uematsu et al. [7]. In this paper, they proposed an approach in which the velocity field was predicted by a turbulence model, the suspension was modeled by the transport equation of drift density, and the saltation was modelled as a function of friction velocity. This approach has also been employed as a basic model in many other studies. An example of the result obtained using this approach is shown in Figure 2. Shortly thereafter, Liston et al. [8] analyzed the snowdrift around a two-dimensional fence using the standard k-ε model. Here, only the saltation was analyzed and the suspension was not. Bang et al. [9] proposed a snowdrift model of two-phase flow called “a generalized drift-flux model” and applied it to three-dimensional snow deposition around a group of houses. They did not consider the influence of saltation. Sundsbø [10] employed the volume-of-fluid method to create a transient development model of snowdrift and conducted a numerical study on two-dimensional snowdrift around a porous fence in flat terrain

using a one-way coupled model, where the air flow controlled the snowdrift within an Eulerian frame of air

<figure><img src="image_9.png" /><figcaption>**Figure 2.** Computed wind scoop around a hut [7].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 3.** Snow accumulation in a building step obtained by numerical simulation (a) without and (b) with a wind deflection fin [10].</figcaption></figure>
<table><caption>**Table 2.** Papers pertaining to CFD studies on snowdrift presented previously at the ICSE.</caption><thead><tr><th>Number</th><th>Year</th><th>Location</th><th>Contributions (Target: Authors)</th></tr></thead><tbody><tr><td>3</td><td>1996</td><td>Sendai, Japan</td><td>Buildings: Kawakami et al. [32], Waechter et al. [33]<br>Step: Sato et al. [34]<br>Fence: Sundsbø and Hansen [35]</td></tr><tr><td>4</td><td>2000</td><td>Trondheim, Norway</td><td>Building group: Thiis and Jaedicke [36]<br>Two cubical obstacles: Thiis and Jaedicke [37]<br>Complex terrain: Jaedicke et al. [38], Sundsbø and Bang [39]<br>Slope and road: Thordarson and Norem [40, 41]</td></tr><tr><td>5</td><td>2004</td><td>Davos, Switzerland</td><td>Mountain hut: Schaelin et al. [42]<br>Elevated structures: Beyers et al. [43]</td></tr><tr><td>6</td><td>2008</td><td>Whistler, Canada</td><td>In-cloud icing on tall buildings: Beyers et al. [44]<br>Building with curved roof: Thiis and Ramberg [45]<br>Cubical obstacle: Tominaga et al. [46]<br>Boundary layer: Okaze et al. [47]</td></tr><tr><td>7</td><td>2012</td><td>Fukui, Japan</td><td>Fences: Potac [48]<br>Boundary layer: Zhou et al. [49], Okaze et al. [50]<br>Four buildings: Ito et al. [51]<br>Model building: Tominaga et al. [52]</td></tr></tbody></table>

velocities. Figure 3 shows the numerical prediction results for snow accumulation on a building step with and without a wind deflection fin. The effect of the deflection fin was clearly demonstrated. Naaim et al. [11] developed an erosion-deposition model in which mass exchange was defined in relation to the flow turbulence, threshold-friction velocity, and snow concentration. They introduced additional terms into the transport equations of $k$ and $\epsilon$ to take into account the effect of snow particles on turbulence. This model has often been used in subsequent studies. Other different approaches have also been proposed. For example, Sekine et al. [12] applied Lagrangian analysis of a snow particle to snowdrift around a two-dimensional building. Tominaga et al. [13] analyzed the snowdrift around an actual apartment building using the revised $k-\epsilon$ model. In this study, because snow blowing into an elevator hall of the building was the point of interest, the transport equation of drift density was solved only for snow particles suspended in the air. The computational results for the trajectory of snow particles are shown in Figure 4.

In the 2000s, more precise studies on the prediction of snowdrift around buildings were performed based on the accumulated knowledge of engineering CFD and snowdrift modeling. Thiis [14] applied the drift-flux model proposed by Bang et al. [9] to snowdrift around a simple building model and examined its performance by comparing the numerical results with full-scale measurements [24]. In this study, the snow depth was not computed directly but the computational results for the friction velocity of mixed fluid were extrapolated to the area where snowdrifts were expected. Beyers et al. [15] predicted transient snowdrift around a cubic structure using the standard $k-\epsilon$ model and compared the numerical results with experimental data obtained in Antarctica [25]. They concluded that the snowdrift predictions matched well with the measured values in terms of both the location and magnitude, but the simulated results for the accumulation near the cube walls were not generally in accordance with observations. Subsequently, Beyers and Waechter [16] obtained application results for snowdrift around more complex structures, as shown in Figure 5, and

<figure><img src="image_11.png" /><figcaption>**Figure 4.** Trajectories of snow particles near a building predicted by CFD [13].</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 5.** Snowdrift simulation around adjacent elevated buildings: (a) CFD simulation and (b) field observation [16].</figcaption></figure>
reported that the prediction results corresponded qualitatively with the actual snowdrift pattern. Tominaga et al. [17] presented the first review of CFD simulations applied to the prediction of snowdrift around buildings. They also presented the CFD prediction results for snowdrift around a cubic building model using a new approach involving a snowdrift model based on experimental and numerical studies. The numerical results were compared with data obtained from detailed field measurements [26] to confirm the accuracy of the new model. More recent studies [18-20] are mentioned in the later sections.

## 2.2 Papers presented at International Conference on Snow Engineering

Important papers on the numerical simulation of snowdrift around buildings have also been presented at previous editions of the International Conference on Snow Engineering (ICSE). The papers from the proceedings that the author felt were relevant are listed in Table 2. The details of the results cannot be explained because of space limitations. Please refer to the original references for further information.

## 3 Recent progress

In this section, noteworthy research topics from recent snowdrift simulation studies are chosen and briefly introduced.

### 3.1 Prediction of roof snow

In the prediction of roof snow loads on buildings, unbalanced snow accumulation due to wind flow is a difficult problem because it involves a complex interaction between snow particle motions and fluid flow affected by the building geometry. So far, only a few attempts have been made to apply CFD to roof snow issues. This et al. [45, 53] predicted the snow distribution on the curved roof of a sports hall located in Oslo and compared the results with measurement results. Figure 6a shows the simulated results for the overall snowdrift pattern on the roof. The location of the end effects can be seen close to the side edge of the roof. Generally, the simulated shape fits the measured shape (Figure 6b) well. This study mainly focused on the reproduction of snow distributions on an actual building under specific weather conditions. There are a few examples in literature wherein the applicability of CFD to snowdrift is examined for a generic configuration. The author recently studied the CFD prediction of snowdrift on a gable-roof building with different roof pitches [54]. In this study, the prediction accuracy of the velocity field around the buildings was validated by comparing the results with wind tunnel experimental data in detail [55]. It should be noted that the correct prediction of the flowfield around roofs with complicated shapes is necessary for the accurate prediction of snowdrift formed on the roof.

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Snow deposition on a curved roof: (a) simulated snowdrift pattern and (b) measured snow depth contours [45, 53].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 7.** Various factors influencing snowdrift formation [28].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 8.** Snow deposition pattern predicted by CFD and an FAE computer snow model [35].</figcaption></figure>
## 3.2. Consideration of meteorological influences

As reviewed in the foregoing sections, most previous studies were basically concerned with only the snowdrift formation due to wind during one drifting snow event. However, snow distributions in built-up environments are affected not only by wind but also by many other factors such as temperature, solar radiation, and artificial heat from buildings. These meteorological effects on snow cover have been investigated for predicting an avalanche (e.g., [56]), but very few studies have considered these effects in the civil and architectural fields. In addition, numerical simulations are expected to be used for predicting more multilateral snow accumulation phenomena, as illustrated in Figure 7, because different conditions can be freely considered. That is, the influence of temporal change in the weather conditions and the heat budget of the snow surface on snow accumulation can be studied.

Waechter et al. [35] predicted the snow distributions around a planned building in Antarctica by CFD and a method called the finite area element (FAE). The FAE is a calculation method in which computational domains are divided into finite elements and the temporal variation in the snow depth at each element is computed according to the weather condition. An example of the result is shown in Figure 8. Furthermore, the effects of the terrain, surroundings, and temporal variations in boundary conditions are very important in precisely predicting snow distribution. For example, the prediction of complicated snowdrift formations such as snow flumes as seasonal characteristics by numerical methods has been very difficult. This is because temporal variations in boundary conditions play an important role in such phenomena. That is, snowdrift simulations driven and updated by hourly meteorological data can be a very useful and important contribution to this field. Therefore, a simulation system is required that can predict and evaluate snow coverage distribution with high accuracy while taking into account the above effects. Later, Beyers et al. [44] proposed a methodology for the prediction and simulation of potential in-cloud icing on buildings. In their method, first, an advanced meteorological analysis of the North American Regional Reanalysis meteorological datasets is performed to predict the supercooled droplet availability, temperature, and prevailing wind direction and speed at multiple elevations above the sea level for a particular site. Then, the in-cloud supercooled droplet collection efficiency is determined through the analysis of the flowfield around a building obstacle by CFD. For this type of analyses, research and development results in the meteorological field are very useful. In addition, obtaining reliable experimental data is also important.

The above prediction methods were not examined by comparing their results with observational results. Tominaga et al. [18] developed a system combining a mesoscale meteorological model that predicted precipitation, including snowfall in an area, and a CFD model that predicted snow phenomena on a building scale. The performance of the system was examined by comparing its results with measured data. Figure 9a outlines the prediction system. It consists of a mesoscale

<figure><img src="image_16.png" /><figcaption>**Figure 9.** Outline of the system used for predicting snow distribution in built-up environments combining a mesoscale meteorological model and a CFD model: (a) schematic structure of the prediction system and (b) analysis procedure of building scale analysis [18].</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Distribution of snow depth: (a) observation and (b) prediction result [18].</figcaption></figure>

meteorological model, snowdrift model, and snowmelt model. Because the models operate independently, each of them can be selected arbitrarily and treated as a module. The analysis procedure of the building scale analysis is illustrated in Figure 9b. In the snowdrift part, the distribution of drifting snow around buildings is simulated by snowdrift analysis based on CFD under snowy weather conditions. In this part, only a transport equation for the drifting snow density is solved under fixed velocity distributions, which are obtained by scaling precomputed wind velocities using the ratio of a reference velocity to the weather conditions. Then, the change in the snow depth attributable to snowmelt is computed by snowmelt analysis based on heat balance analysis on the snow surface under weather conditions including thermal
components such as temperature. The final snow depth distributions can be obtained by summing the two snow depth changes obtained from the two analyses. The snowdrift patterns, i.e., erosion around the upwind corners and deposition in front of and behind a building, obtained by the model show good correspondence with those obtained from field observations, as shown in Figure 10. Later, Ito et al. [50] presented a more advanced system for predicting the spatial distribution of snow depth in urban areas by combining snowdrift and snowmelt models. Although the snowdrift model is the same as that in the previous system, the snowmelt model predicts the snowmelt rate by solving the heat balance equation for the snow surface. The heat balance calculation adopts a coupled approach of convection, conduction, and radiation, including multi-reflections of shortwave and longwave radiation, to accurately predict the heat balance on urban surfaces within the complex wind and thermal environments around buildings.

## 3.3 Application of CFD results to practical design

As described earlier, CFD can provide detailed information on relevant flow variables in the entire calculation domain under well-controlled conditions and without similarity constraints. Therefore, CFD is suitable for considering countermeasures and counterproposals to design problems. Okaze et al. [57] studied the contribution ratios of the snow particles falling from the sky and those from the ground surface to the total snowdrift. These ratios were calculated by solving two transport equations of the two types of snow drifting densities with different settling velocities. Figure 11 shows the normalized snow depth (Figure 11a) and the contribution ratio of snow from the ground surface to the snowdrift (Figure 11b) predicted by CFD around a cubic building. Although the values of the contribution ratio are large at locations laterally away from the building sides, the values are small near the building. This means that the snowdrift near the building is mainly caused by the snow falling from the sky, but that at the laterally distant regions are mainly caused by the snow blown from the ground surface. Such information is quite useful for providing an effective countermeasure to the snowdrift problem around buildings. As a different example, Thiis and Ferreira [19] studied the efficiency of windbreaks comprising three-dimensional arrays of pillars, which were used for snow deposition. The air flow and snow deposition around three different geometrical configurations of pillars were analyzed by CFD. As an example of the result, a predicted snow deposition pattern around windbreaks is shown in Figure 12. Based on the results, an arrangement can be proposed, which can be regarded as an alternative for landscape architects aiming to shelter a specific area. Based on the reliability of CFD simulations, it is important to utilize the extensive information provided by CFD in practical design.

<figure><img src="image_18.png" /><figcaption>**Figure 11.** Horizontal distribution of (a) normalized snow depth and (b) contribution ratio of snow from the ground surface to the total snowdrift predicted by CFD around a cubic building [57].</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 12.** Snow deposition pattern around windbreaks predicted by CFD [19].</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 13.** Schematic representation of the flow around an isolated low-rise building ([58]; modified from [59]).</figcaption></figure>
# 4 Future perspectives

## 4.1. Model development

The concepts and parameters used in snowdrift modeling have been mostly derived empirically from observations and experiments. Many of the observations and experiments were based on a simple drifting boundary layer developed on a flat plane. However, the flowfield around buildings includes accelerations and decelerations due to separation and recirculation flows, as shown in Figure 13. For example, an estimation model for the transport rate of drifting snow attributable to saltation $Q_{sal}$ proposed by Pomeroy and Gray [22] has often been used in many previous studies. The constant in the model was determined by observation in a drifting boundary layer. However, it should be noted that the proposed relation gives the saturated transport rate of drifting snow based on the assumption that drifting snow in the saltation layer has reached its equilibrium state, where neither wind velocity distribution nor $Q_{sal}$ changes in the streamwise direction. Okaze et al. [27, 47] investigated the change in the spatial distribution of $Q_{sal}$ for non-equilibrium flow in a developing saltation layer and its relationship with wind velocity by performing a wind tunnel experiment using actual snow particles [60]. This experiment was performed on a “loose” snow surface, which was covered with soft snow such as fresh snow (with almost no cohesion between particles); here, snow particles jumped out from the surface because of wind, and a saltation layer developed. Figure 14 provides an outline of the experiment. There were two experimental conditions under which the reference wind velocity $U_0$ was changed. The relationship between the fetch distance and $Q_{sal}$ is shown in Figure 15. The two solid lines parallel to the horizontal axis in the figure show the predicted value of $Q_{sal}$, which was obtained by substituting the presumed $U^*$ (based on the downstream velocity profile where the boundary layer was fully developed) into Pomeroy’s relation. The values obtained from Pomeroy’s relation are greatly overestimated in almost the entire region with $U_0 = 5.0$ m/s and in the windward region with $U_0 = 7.0$ m/s. Thus, it is problematic to apply Pomeroy’s relation (which is based on the equilibrium condition) to a flow around a building in which drifting snow is in non-equilibrium because of acceleration and deceleration. This suggests that the empirical equation of the snowdrift rate used in previous studies does not work appropriately for snowdrift accompanied by rapid acceleration and deceleration of flow such as the flow around buildings. It is necessary to further consider the applicability and validity of the modeling of sub-processes in a saltation layer at the developing stage (non-equilibrium state).

Furthermore, it should be mentioned that the CFD modeling of snowdrift applied previously is based on a Reynolds-averaged Navier-Stokes (RANS) equation model approach. This is expected to be useful for the practical application of CFD to snowdrift around a building. However, in actual situations, the deposition and erosion of snow are closely related to the real downburst and splash effects present in the surface boundary layer.

<figure><img src="image_21.png" /><figcaption>Figure 14. Outline of the experiment [27, 47].</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 15. Relationship between the fetch distance obtained from experimental results and $Q_{sal}$ [27, 47].</figcaption></figure>

That is, it is important to take the instantaneous turbulent flow structures that generate sweep and ejection events into account in snowdrift modeling. This may mean that RANS approaches are not suitable and large eddy simulation (LES) is required to extract non-isotropic three-dimensional velocity fluctuations. Further investigation should consider such effects on snowdrift modeling.

## 4.2. Validation database

The accuracy and reliability of CFD simulations are points of concern, and solution verification and validation studies are imperative. Verification is defined as the process of assessing numerical uncertainty and estimating the sign and magnitude of the numerical error itself and the uncertainty in that error estimate. On the other hand, validation is defined as the process of assessing modeling uncertainty by using benchmark experimental data and estimating the sign and magnitude of the modeling error itself. In the validation process, reliable measurement data, wherein the boundary conditions are clearly defined, are indispensable for examining prediction results. Unfortunately, only a few studies have provided such reliable measurement data [26, 61]. Field measurements and observations are important because they are based on actual phenomena. However, their results are affected by numerous parameters and the conditions are not controllable. Therefore, advanced experimental facilities, in which various influencing factors such as the wind velocity, temperature, and snow fall rate are fully controllable, have been developed (e.g., [60, 62, 63]). The extensive utilization of such facilities is expected to provide important information for snowdrift modeling. Furthermore, in the near future, it is expected that a validation database will be opened for researchers worldwide and will be available online. In addition,
numerical simulations have certain inherent limitations. The experimental results should be not only used as a validation database but also considered to be complementary to simulation results. Therefore, a "hybrid" approach, in which both methods are complementarily used, is very important to elucidate complex phenomena such as snowdrift.

### 4.3. Best practice guidelines

In CFD simulations, a large number of choices need to be made by the user. It is well known that these choices can have a very large impact on the results. It is of course apparent that the existing practical guidelines for wind engineering [64-67] and general CFD application [68-70] are important and useful for snowdrift modeling. The correct prediction of velocity fields is necessary for obtaining accurate snowdrift calculation results. However, because only few examples of snowdrift around buildings are available currently, the influence of parameters in a snowdrift model on the results is not clear. It is suggested that best practice guidelines for snowdrift calculations by CFD should be provided by accumulating comprehensive sensitive analysis results in the future.

## 5 Conclusions

This article reviewed previous CFD studies on snowdrift around buildings conducted with different research purposes, configurations, boundary conditions, and modeling approaches to provide a perspective on CFD modeling and give some insights into future directions for practical applications. This review provided an overview of CFD applications to snowdrift around buildings in recent years. Clearly, CFD has considerable potential, as established in this review. However, to use it appropriately, sufficient attention should be paid to modeling theory and implementation. Because snowdrift is a highly complex phenomenon, we should not excessively rely on numerical simulation. It is important to use data from actual phenomena for all aspects of study. Furthermore, it is also important to respect traditional approaches such as field observations and physical experiments and to seek a complimentary approach that combines the advantages of experimental and simulation techniques.

## Acknowledgements

The author is very grateful to Dr. Philippe Delpech and Prof. Dr. Ing. Thomas K. Thiis, who are the organizers of the Eighth International Conference on Snow Engineering to be held in Nantes, France, for their kind invitation to deliver a keynote lecture and for their efforts in aiming to make this conference a great success.

## References

1. J.D. Iversen, Comparison of wind-tunnel model and full-scale snow fence drifts, J. Wind Eng. Ind. Aerodyn. 8(3), 231-249 (1981)

2. F.D.M. Sant'Anna, D.A. Taylor, Snow drifts on flat roofs: wind tunnel tests and field measurements, J. Wind Eng. Ind. Aerodyn. 34(3), 223-250 (1990)

3. N. Isyumov, M. Mikitiuk, Wind tunnel model tests of snow drifting on a two level flat roof, J. Wind Eng. Ind. Aerodyn. 36(2), 893-904 (1990)

4. D.J. Smedley, K.C.S. Kwok, D.H. Kim, Snowdrifting simulation around Davis Station workshop, Antarctica, J. Wind Eng. Ind. Aerodyn. 50, 153-162 (1993)

5. M. O'Rourke, A. DeGaetano, J.D. Tokarczyk, Snow drifting transport rates from water flume simulation, J. Wind Eng. Ind. Aerodyn. 92(14-15), 1245-1264 (2004)

6. B. Blocken, 50 years of Computational Wind Engineering: Past, present and future, J. Wind Eng. Ind. Aerodyn. 129, 69-102 (2014)

7. T. Uematsu, T. Nakata, K. Takeuchi, Y. Arisawa, Y. Kaneda, Three-dimensional numerical simulation of snowdrift, Cold Reg. Sci. Technol. 20, 65-73 (1991)

8. G.E. Liston, R.L. Brown, J.D. Dent, A two-dimensional computational model of turbulent atmospheric surface flows with drifting snow, Ann. Glaciol. 18, 281-286 (1993)

9. B. Bang, A. Nielsen, P.A. Sundsbø, T. Wiik. Computer simulation of wind speed, wind pressure and snow accumulation around buildings (SNOW-SIM), Energ. Buildings 21, 235-243 (1994).

10. P.A. Sundsbø, Numerical simulations of wind deflection fins to control snow accumulation in building steps, J. Wind Eng. Ind. Aerodyn. 74-76, 543-552 (1998)

11. M. Naaim, F. Naaim-Bouvet, H. Martinez, Numerical simulation of drifting snow: erosion and deposition models, Ann. Glaciol. 26, 191-196 (1998)

12. A. Sekine, M. Shimura, A. Maruoka, H. Hirano, The numerical simulation of snowdrift around a building. Int. J. Comput. Fluid D. 12, 249-255 (1999).

13. Y. Tominaga, A. Mochida, CFD prediction of flowfield and snowdrift around a building complex in a snowy region, J. Wind Eng. Ind. Aerodyn. 81, 273-282 (1999)

14. K.T. Thiis, A comparison of numerical simulations and full-scale measurements of snowdrifts around buildings, Wind Struct. 3(2), 73-81 (2000)

15. 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, 725-747 (2004)

16. M. Beyers, B. Waechter, Modeling transient snowdrift development around complex three-dimensional structures, J. Wind Eng. Ind. Aerodyn. 96, 1603-1615 (2008)

17. Y. Tominaga, T. Okaze, A. Mochida, CFD modeling of snowdrift around a building: An overview of models and evaluation of a new approach, Build. Environ. 46(4), 899-910 (2011)
18. Y. Tominaga, A. Mochida, T. Okaze, T. Sato, M. Nemoto, H. Motoyoshi, S. Nakai, T. Tsutsumi, M. Otsuki, T. Uamatsu, H. Yoshino, Development of a system for predicting snow distribution in built-up environments: Combining a mesoscale meteorological model and a CFD model, J. Wind Eng. Ind. Aerodyn. 99, 460-468 (2011)

19. T. Thiis, A.D. Ferreira, Sheltering effect and snow deposition in arrays of vertical pillars, Environ. Fluid Mech. 15, 27-39 (2015)

20. T. Okaze, Y. Takano, A. Mochida, Y. Tominaga, Development of a new k-ε model to reproduce the aerodynamic effects of snow particles on a flow field, J. Wind Eng. Ind. Aerodyn. 144, 118-124 (2015)

21. J.D. Iversen, R. Greeley, B.R. White, J.B. Pollack, Eolian erosion of the Martian surface, Part 1; Erosion rate similitude. Icarus 26(3), 321-331 (1980)

22. J.W. Pomeroy, D.M. Gray, Saltation of snow, Water Resour. Res. 26(7), 1583-1594 (1990)

23. R.S. Anderson, P.K. Haff, Simulation of Eolian Saltation, Science 241, 820-823 (1988)

24. T.K. Thiis, Y. Gjessing, Large-scale measurements of snowdrifts around flat-roofed and single-pitch-roofed buildings, Cold Reg. Sci. Technol. 30(1-3), 175-181 (1999)

25. J.H.M. Beyers, T.M. Harms, Outdoors modelling of snowdrift at SANAE IV Research Station, Antarctica, J. Wind Eng. Ind. Aerodyn. 91(4), 551-569 (2003)

26. S. Oikawa, T. Tomabechi, Daily observation of snowdrifts around a model cube, *Snow Engineering*, *Hjorth-Hansen, Holand, Loset & Norem (eds)*, 369-375 (2000)

27. T. Okaze, A. Mochida, Y. Tominaga, M. Nemoto, T. Sato, Y. Sasaki, K. Ichinohe, Wind tunnel investigation of drifting snow development in a boundary layer, J. Wind Eng. Ind. Aerodyn. 104-106, 532-539 (2012)

28. G.E. Liston, M. Sturm, A snow-transport model for complex terrain, J. Glaciol. 44(148), 498-516 (1998)

29. E.M. Greene, G.E. Liston, R.A. Pielke Sr., Simulation of above treeline snowdrift formation using a numerical snow-transport model, Cold Reg. Sci. Technol. 30(1-3), 135-144 (1999)

30. G.E. Liston, R.B. Haehnel, M. Sturm, C.A. Hiemstra, S. Berezovskaya, R.D. Tabler, Simulating complex snow distributions in windy environments using SnowTran-3D, J. Glaciol. 53(181), 241-256 (2007)

31. S. Murakami, Computational wind engineering, J. Wind Eng. Ind. Aerodyn. 36(1), 517-538 (1990)

32. S. Kawakami, T. Uematsu, T. Kobayashi, Y. Kaneda, Numerical study of a snow wind scoop, *Snow Engineering: Recent Advances, Izumi, Nakamura & Sack (eds)*, 55-62 (1997)

33. P.-A. Sundsbo, E.W.M. Hansen, Modelling and numerical simulation of snow drift around snow fences, *Snow Engineering: Recent Advances, Izumi, Nakamura & Sack (eds)*, 353-359 (1997)

34. T. Sato, T. Uematsu, Y. Kaneda, Application of random walk model to blowing snow, *Snow Engineering: Recent Advances, Izumi, Nakamura & Sack (eds)*, 133-138 (1997)

35. B.F. Waechter, R.J. Cinclair, G.D. Schuyler, C.J. Williams, Snowdrift control design: Application of CFD simulation techniques, *Snow Engineering: Recent Advances, Izumi, Nakamura & Sack (eds)*, 511-516 (1997)

36. C. Jaedicke, T.K. Thiis, B. Bang, The snowdrift pattern around a small hill in the High Arctic, *Snow Engineering: Recent Advances, Izumi, Nakamura & Sack (eds)*, 75-80 (2000)

37. T.K. Thiis, C. Jaedicke, Changes in the snow drift pattern caused by building extension - Investigations through scale modelling and numerical simulations, *Snow Engineering, Hjorth-Hansen, Holand, Loset & Norem (eds)*, 363-368 (2000)

38. T.K., Thiis, C. Jaedicke, The snowdrift pattern around two cubical obstacles with varying distance - Measurements and numerical simulations, *Snow Engineering, Hjorth-Hansen, Holand, Loset & Norem (eds)*, 369-375 (2000)

39. P.A. Sundsbo, B.Bang, Snow drift control in residential areas - Field measurements and numerical simulations, *Snow Engineering, Hjorth-Hansen, Holand, Loset & Norem (eds)*, 377-382. (2000)

40. S. Thordarson, H. Norem, Simulation of two-dimensional wind flow and snow drifting application for roads: Part1, *Snow Engineering, Hjorth-Hansen, Holand, Loset & Norem (eds)*, 437-443 (2000)

41. S. Thordarson, H. Norem, Simulation of two-dimensional wind flow and snow drifting application for roads: Part2, *Snow Engineering, Hjorth-Hansen, Holand, Loset & Norem (eds)*, 445-452 (2000)

42. A. Schaelin, L. Ilg, M. Benesch, Snow deposition around mountain hut - design optimization by CFD and scaled water channel model and realization of solutions, *Snow Engineering V, Bartelt, Adams, Christen, Sack & Sato (eds)*, 147-155 (2004)

43. J.H.M. Beyers, T.M. Harms, P.A. Sundsbø, Numerical simulation of snow drift around an elevated structure, *Snow Engineering V, Bartelt, Adams, Christen, Sack & Sato (eds)*, 185-191 (2004)

44. J.H.M. Beyers, X. Qiu, R. Hakimi, Meteorological prediction and CFD simulation of in-cloud icing on tall buildings, *In: Proceedings of Snow Engineering VI, Whistler, Canada (2008)*

45. T.K. Thiis, J.F. Ramberg, Measurements and numerical simulations of development of snow drifts of curved roofs, *In: Proceedings of Snow Engineering VI, Whistler, Canada (2008)*

46. Y. Tominaga, T. Okaze, A. Mochida, T. Shida, H. Yoshino, CFD prediction of snowdrift around a cubic building model, *In: Proceedings of Snow Engineering VI, Whistler, Canada (2008)*

47. T. Okaze, A. Mochida, Y. Tominaga, M. Nemoto, Y. Ito, T. Shida, T. Sato, H. Yoshino, Modeling of drifting snow development in a boundary layer and its effects on wind field, *In: Proceedings of Snow Engineering VI, Whistler, Canada (2008)*

48. J. Potac, A 2D numerical modeling of snowdrift development, *In: Proceedings of the 7th International Conference on Snow Engineering (Snow Engineering VII)*, Fukui, Japan (2012)
49. X. Zhou, X. Li, M. Gu, Y. Zhang, P. Huang, A new two-equation model for numerical simulation on snow drifting, In: Proceedings of the 7th International Conference on Snow Engineering (Snow Engineering VII), Fukui, Japan (2012)

50. Y. Ito, A. Mochida, T. Okaze, Y. Sasaki, Y. Tominaga, Development of a prediction method for snow-depth distribution in urban areas based on the coupling of snowdrift and snowmelt models, In: Proceedings of the 7th International Conference on Snow Engineering (Snow Engineering VII), Fukui, Japan (2012)

51. Y. Tominaga, T. Okaze, A. Mochida, M. Nemoto, T. Sato, T. Tsutsumi, Y. Sasaki, Seasonal change of snow distribution around a building: comparison between numerical simulation and filed measurements, In: Proceedings of the 7th International Conference on Snow Engineering (Snow Engineering VII), Fukui, Japan (2012)

52. T. Okaze, A. Mochida, Y. Tominaga, Y. Takano, Numerical modeling of snow particles on turbulent flowfield based on canopy model, In: Proceedings of the 7th International Conference on Snow Engineering (Snow Engineering VII), Fukui, Japan (2012)

53. T. K. Thiis, J. Potac, J.F. Ramberg, 3D numerical simulations and full scale measurements of snow depositions on a curved roof, In 5th European & African Conference on Wind Engineering, Florence (2009)

54. Y. Tominaga, T. Okaze, A. Mochida, CFD simulation of drift snow loads for an isolated gable-roof building, 8th International Conference on Snow Engineering, Nantes, France (2016) (to be submitted)

55. Y. Tominaga, S. Akabayashi, T. Kitahara, Y. Arinami, Air flow around isolated gable-roof buildings with different roof pitches: Wind tunnel experiments and CFD simulations, Build. Environ. 84, 204-213 (2015)

56. P.B. Bartelt, M. Lehning, A physical SNOWPACK model for avalanche warning services. Part I: numerical model, Cold Reg. Sci. Technol. 35, 123-145 (2002)

57. T. Okaze, A. Mochida, Y. Tominaga, Y. Ito, H. Yoshino, CFD prediction of snowdrift around a cube using two transport equations for drifting snow density, In 5th International Symposium on Computational Wind Engineering (CWE2010), Chapel Hill, North Carolina, USA (2010)

58. B. Blocken, T. Stathopoulos, J. Carmeliet, J.L.M. Hensen, Application of computational fluid dynamics in building performance simulation for the outdoor environment: an overview, J. of Build. Perform. Simu. 4(2), 157-184 (2011)

59. R.P. Hosker, Flow and diffusion near obstacles. In:D. Randerson, ed. Atmospheric Science and Power Production. (Office of Scientific and Technical Information, United States Department of Energy 1984)

60. T. Sato, K. Kosugi, A. Sato, Saltation-layer structure of drifting snow observed in wind tunnel, Ann. Glaciol. 32, 203-208 (2001)

61. T. Tutsumi, T. Chiba, T. Tomabechi, Snowdrifts on and around buildings based on field measurement, In: Proceedings of the 7th International Conference on Snow Engineering (Snow Engineering VII), Fukui, Japan (2012)

62. P. Delpech, P. Palier, J. Gandemer, Snowdrifting simulation around Antarctic buildings, J. Wind Eng. Ind. Aerodyn. 74-76, 567-576 (1998)

63. P. Delpech, T. Thiis, Applications of “snowwind” engineering – climatic wind tunnel methods, Technical Transactions iss. 12. Civil Engineering iss. 2-B, 381-403 (2015)

64. J. Franke, C. Hirsch, A.G. Jensen, H.W. Krüs, M. Schatzmann, P.S. Westbury, S.D. Miles, J.A. Wisse, N.G. Wright, Recommendations on the use of CFD in wind engineering, In: van Beeck, J.P.A.J. (Ed.), Proceedings of the International Conference on Urban Wind Engineering and Building Aerodynamics. COST ActionC14, Impact of Wind and Storm on City Life Built Environment. Von Karman Institute, Sint-Genesius-Rode, Belgium (2004)

65. J. Franke, A. Hellsten, H. Schlünzen, B. Carissimo, (Eds.), Best practice guideline for the CFD simulation of flows in the urban environment (COST Office, Brussels, 3-00-018312-4, 2007)

66. Y. Tominaga, A. Mochida, R. Yoshie, H. Kataoka, T. Nozu, M. Yoshikawa, T. Shirasawa, AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings, J. Wind Eng. Ind. Aerodyn. 96(10-11), 1749-1761 (2008)

67. B. Blocken, Computational Fluid Dynamics for urban physics: Importance, scales, possibilities, limitations and ten tips and tricks towards accurate and reliable simulations, Build. Environ. 91, 219-245 (2015)

68. M. Casey, T. Wintergerste, ERCOFTAC Special Interest Group on Quality and Trust in Industrial CFD, Best Practice Guidelines (ERCOFTAC, Brussels, 2000)

69. W.L. Oberkampf, T.G. Trucano, C. Hirsch, Verification, validation, and predictive capability in computational engineering and physics. Appl. Mech. Rev. 57(5), 345-384 (2004)

70. C. Hirsch, The Development of a Framework for CFD Validation and Best Practice: The QNET-CFD Knowledge Base, Chinese J. Aeronaut. 19(2), 105-113 (2006)
# Improvements of the viscous treatment of the snow phase in two-way coupled Eulerian-Eulerian simulations of drifting snow

Ziad Boutanios<fnref n="1" />,<fnref n="2" />,⁸ and Hrvoje Jasak<fnref n="2" />

**Abstract.** This paper presents an extension to a viscous model of the snow phase developed by the authors. The original model is developed from first principles, based on analogy between Lagrangian and Eulerian representations of the snow phase in the creeping layer. Here, the original model is extended to the snowbed, the saltation and the suspension layers, and both original and new viscous models are compared and discussed. The original formulation consists of a two-way coupled four-equation Navier-Stokes system with all drifting snow layers resolved, for a single snow particle size. Validations of the original viscous model against controlled wind tunnel drifting snow experiments showed good agreement in the snow fluxes in the saltation layer, provided the separately computed contributions of several snow particle size classes are combined in a two-parameter Gamma particle size distribution. The extended viscosity model is seen to greatly improve the predictions of snow fluxes in the suspension layer, providing a good match to experimental measurements. No improvements are seen in the airflow velocity and turbulent kinetic energy profiles since the momentum coupling term between phases is not affected by the improvements. Recommendations are made for future work.

## 1 Introduction

Several drifting snow models are available in the literature, most of which are based on Reynolds Averaged Navier-Stokes (RANS) formulations in the Eulerian-Eulerian and Eulerian-Lagrangian frames regarding the air and snow phases. Both approaches can yield reasonable results compared to experiments but the Eulerian-Eulerian approach requires a lower computational effort since Lagrangian particle tracking requires a great deal of particles to yield statistically meaningful results [1].

Presently, Eulerian-Eulerian modelling of drifting snow is based on two main approaches: the transport of snowdrift density approach and the Volume of Fluid (VOF) approach. The snow density transport approach solves a scalar transport equation of the snow density in the suspension layer where the motion of the snow phase does not affect the airflow (one-way coupling), in addition to the one-way coupled airflow continuity and momentum equations. Instead of solving the snow momentum equation, the snow particle fall velocity is either assumed constant [2, 3], or calculated using a buoyancy-based formula derived from the equations of motion of spherical particle in steady-state airflow [4]. Another variation on the snow density approach solves the mixture continuity and momentum equations in two-way coupled form [3]. Different treatments of turbulence are used such as a Prandtl mixing layer [2], a standard k-epsilon model with a particle effect correction [3], and a Launder-Kato k-epsilon model with particle effect correction as well [4]. All these methods represent the saltation layer as a boundary condition on the snowbed surface obtained from an equilibrium empirical formulation of the transport rate of drifting snow in the saltation layer [5]. This empirical formulation was shown to overpredict the transport rate of snow in non-

equilibrium situations [6]. An exception is the approach of Tominaga et al. [4] that resolves the saltation layer with the same equations as the suspension layer. Another notable exception is a two-way coupled formulation which fully resolves both saltation and suspension layers without equilibrium empirical boundary conditions [7]. However, the snow saltation layer transport equations are based on self-similarity equilibrium relationships between both phases.

The VOF approach is similar to the snow density approach in that it solves another one-way coupled system of equations consisting of the mixture continuity and momentum equations, in addition to a scalar transport equation of the volume fraction of one of the phases. The latter equation is equivalent to a continuity equation and the drift velocity between the phases is set to different constant values in the different snow transport layers [8] or calculated with a drift flux term [9]. The drift flux term [10] is strictly valid for negligible particle inertia which applies only to the smallest particles transported by turbulence into the suspension layer, and not the larger ones in the saltation layer. All the approaches previously described calculate the snow depth by explicitly accounting for the contributions of saltation and suspension fluxes based on experimental equilibrium relationships.

A two-way coupled approach is taken by Boutanios and Jasak [11] where separate Navier-Stokes systems are solved for each phase and coupling between the phases is accounted for through drag and lift. The Navier-Stokes system solved is:

$$ \frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot \alpha_k \rho_k \mathbf{u}_k = 0 \quad (1) $$

$$ \frac{\partial \alpha_k \rho_k \mathbf{u}_k}{\partial t} + \nabla \cdot \alpha_k \rho_k \mathbf{u}_k \mathbf{u}_k = -\alpha_k \nabla p + \nabla \cdot \alpha_k \boldsymbol{\tau}_k + \alpha_k \rho_k \mathbf{g} + \mathbf{M}_k \quad (2) $$

Here, k=1,2 is the index of either phase, $\alpha_k$ is the phase volume fraction, $\rho_k$ is the phase density, $\boldsymbol{\tau}_k$ is the phase stress tensor, $\mathbf{u}_k$ is the phase velocity vector, p is the

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

<footnote n="1">¹Binkz Incorporated, Laval, Quebec, Canada</footnote>

<footnote n="2">²CFD Lab, FSB, University of Zagreb, Zagreb, Croatia</footnote>

<footnote>a Corresponding author: ziad@binkz.ca</footnote>
static pressure field, and g is the gravitational acceleration vector. M_k is the momentum exchange term between the phases which consists of the drag and lift forces, found by scale analysis to be the dominant forces for drifting snow [7]. Turbulence is accounted for with the standard k-ε model, without particle effects. The snow phase shear stress term in the momentum equation above requires a snow viscosity parameter, so one was developed from first principles based on an analogy between a Lagrangian rolling particle on the snowbed and a corresponding Eulerian flowing fluid. This analogy is shown in Figure 1. The two-way coupled viscous approach showed good agreement with experimentally measured snow fluxes in a controlled wind tunnel experiment [12], provided the distribution of particle sizes is accounted for by statistically combining the results of several single diameter simulations. Boutanios and Jasak [11] do not need to specifically compute the change in snow depth since convection and diffusion fluxes are already part of the Navier-Stokes system of equations.

<figure><img src="image_24.png" /><figcaption>Figure 1. Lagrangian particles in Eulerian control volume.</figcaption></figure>

In this paper, the derivation of the viscosity model is extended to the snowbed, and the upper saltation layer and the suspension layer. Comparisons of the numerical results of the extended viscosity model to the numerical results the original viscosity model, and to the experimental measurements are provided and discussed.

## 2 Extended viscosity model

The extended viscosity model consists of Equations (3), (4) and (5), applying to the snowbed, the creeping/saltation layer, and the suspension layer respectively.

$$v_s = \begin{cases} \frac{\tau_t}{2\rho_s \varepsilon_s} & \text{in snowbed} & (3) \\ -\frac{F_{d,v}d_p}{6\alpha_a\varepsilon_s} + \frac{\tau_t}{2\rho_s\varepsilon_s} & \text{in creeping/saltation layer} & (4) \\ C_t^2 v_t & \text{in suspension layer} & (5) \end{cases}$$

Here, $v_s$ is the snow phase kinematic viscosity, $\rho_s$ is the snow density, $\tau_t$ is the drifting threshold shear stress, $\varepsilon_s$ is the snow rate of strain, $F_{d,v}$ is the drag force per unit

volume, $d_p$ is the snow particle diameter, $\alpha_a$ is the air volume fraction, $v_t$ is the air turbulent viscosity, and $C_t$ is the turbulent viscosity multiplier as determined by Gosman et al. [13] for dilute particle-laden flows. Since the flow in the suspension layer is one-way coupled, we use here a constant $C_t=1$. The original viscosity model consists of Equation (4) which is strictly applicable in the creeping and saltation layers. However, tests have shown that it can be used in the snowbed as well by retaining only the threshold shear stress contribution to the viscosity as shown in Equation (3). In the suspension layer the flow is extremely dilute and one-way coupled, with transport of the snow provided by turbulent lift. There, the snow phase flow is essentially inviscid, but turbulent transport is accounted for in a Boussinesq fashion through $C_t$ the turbulent viscosity multiplier, as shown in Equation (5). Validation results are presented in the next section.

## 3 Validation against experimental data

### 3.1 Experimental setup

Validation of the extended viscosity model is done against the same experimental data used to validate the original model, as well as the original model itself. The experimental data was measured in a controlled wind tunnel drifting snow experiment carried out by the Tohoku university group [14] in Japan. In the Tohoku experiments, the profiles of mass flux of drifting snow, wind velocity, and turbulence statistics were measured under several different wind speed conditions, and at four different measurement stations in the wind tunnel working section shown in Figure 2. The measurement stations were located at X=3, 6, 9, and 11.5 m. Loose snow is deposited in a 2 cm deep gutter along the length of the tunnel floor which is preceded by a 1 m fetch of frozen non-drifting hard snow, in order to generate a non-equilibrium saltation layer. The measured airflow velocity and turbulence kinetic energy at station X=0m were used as boundary conditions for the numerical simulations. These profiles are shown in Figure 3 where $Z_t=0.2$ m and $U_r \approx 7.0$ m/s are the reference height and airflow velocity.

<figure><img src="image_25.png" /><figcaption>Figure 2. Side view of the mid-section of the wind tunnel working section (flow direction from left to right).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 3.** Numerical inlet boundary conditions, as measured in the wind tunnel at X=0 m.</figcaption></figure>

Samples of the snow particles used in the experiment are shown in Figure 4. These are quite irregular and bulky, rarely smaller than 0.50 mm in either length or width, which implies a particle distribution with an average diameter larger than 0.50 mm. It should be noted that particle sizes occurring in natural snow and mechanical breakage phenomena are best described with a two-parameter Gamma Probability Distribution Function (PDF), be it as aggregate on the ground [15] or drifting above it [16].

<figure><img src="image_27.png" /><figcaption>**Figure 4.** Sample of the experiment snow particles.</figcaption></figure>

## 3.2 Numerical setup

The 2D computational mesh used for the simulations is shown in Figure 5. The top part of the figure is a close-up of the mesh at the inlet of the computational domain of the wind tunnel at X=0m, with the loose snow layer in the gutter shown in white. The mesh is fully structured hexahedral with a transverse element size of 4 mm in the gutter and at the top of the tunnel. The longitudinal element size in the flow direction along the X-axis is about 6 cm. Tests were carried out with a mesh twice finer in both directions and the results varied by less than 15% so results obtained from the present mesh can be reasonably considered mesh-independent. Simulations were carried out at a threshold shear stress of 0.052 Kg/(m.s2) which is the minimum experimentally observed value for drifting.

<figure><img src="image_28.png" /><figcaption>**Figure 5.** Close-up of the numerical mesh near the tunnel inlet showing the snowbed in the gutter (above), and an overall view of the mesh (below).</figcaption></figure>

Seven transient simulations were carried out for constant particle diameters ranging from 0.1–1.3 mm. Boutanios and Jasak [11] show that a statistical combination of the individual solutions is required to account for the natural particle size segregation occurring normal to the snowbed in drifting snow [16]. Here, the individual simulations are accordingly combined into one average diameter solution, according to the particle class weights of a two-parameter Gamma distribution with an average diameter of 0.7 mm, shown in Figure 6 below.

<figure><img src="image_29.png" /><figcaption>**Figure 6.** Two-parameter Gamma distribution with average diameter of 0.7 mm.</figcaption></figure>

The average diameter solution was then compared to the experimentally measured profiles for both original and extended viscosity models. It should be pointed out that the experimental profiles were measured within a 30s window at all four measurement stations, under negligible bed erosion conditions. The same approach was adopted in extracting the numerical profiles from the simulation data for the comparisons to be valid.

## 3.2 Snow flux profiles

The snow flux profiles for the 0.7mm average diameter solutions of the original viscosity model at X=11.5m is shown in Figure 7 below.
<figure><img src="image_30.png" /><figcaption>**Figure 7.** Numerical snow flux profile vs. experimental measurements, for $d_{avg}=0.7$mm at $X=11.5$m; original viscosity model.</figcaption></figure>

The numerical results indicate that the original viscosity model does a good job in matching the experimental measurements in the saltation layer, which extends up to about $Z=4$cm. The simulation snow fluxes are also seen to be stable enough to provide a 30s measurement window during which a stable snowbed can be observed, thus satisfying the comparison requirements to the experimental measurements. The results also imply that the particle average size at $X=11.5$m in the experiments could be around 0.7mm. This makes sense given the snow particle sample shown in Figure 4. The agreement between experiment and simulations is not good in the suspension layer, above $Z=4$cm, since transport by suspension is due to turbulent lift and that effect is not included in the original viscosity model.

The same snow flux profiles for the 0.7mm average diameter solution for the extended viscosity model at $X=11.5$m is shown in Figure 8. The extended viscosity treatment into the suspension layer, based on the Gosman Ct model, is seen to produce substantial amounts of turbulent lift across both saltation and suspension layers. The 0.7mm average size solution is much closer to the experimental measurements and exhibits the stable conditions required to match the experimental measurement window.

<figure><img src="image_31.png" /><figcaption>**Figure 8.** Numerical snow flux profile vs. experimental measurements for $d_{avg}=0.7$mm at $X=11.5$m, extended viscosity model.</figcaption></figure>

## 3.2 Non-dimensional airflow velocity profiles

The non-dimensional airflow velocity profiles for the 0.7mm average diameter solutions of the original viscosity model at $X=11.5$m is shown in Figure 9 below.

<figure><img src="image_32.png" /><figcaption>**Figure 9.** Numerical non-dimensional airflow velocity profiles vs. experimental measurements, for $d_{avg}=0.7$mm at $X=11.5$m; original viscosity model.</figcaption></figure>

The simulation results exhibit the required stability for taking the experimental measurements and the shapes of the curves are very similar to the experimental curve shape. However, there is a large quantitative difference between simulations and experiment.

The non-dimensional airflow velocity profiles for the 0.7mm average diameter solutions of the extended viscosity model at $X=11.5$m are shown in Figure 10 below. The results of the extended viscosity model are very close to those of the original viscosity model, and no improvement is noticed as with the snow flux profiles. This discrepancy is due to a known deficiency of the $k-\epsilon$ turbulence model which is used here. The turbulent kinetic energy is overestimated by $k-\epsilon$, and this excess turbulent kinetic energy can only come from a velocity deficit based on the model's formulation. The excess
turbulent kinetic energy production can be assessed in Figure 11 and Figure 12 below.

<figure><img src="image_33.png" /><figcaption>**Figure 10.** Numerical non-dimensional airflow velocity profile vs. experimental measurements for davg=0.7mm at X=11.5m, extended viscosity model.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 11.** Numerical non-dimensional TKE profile vs. experimental measurements for davg=0.7mm at X=11.5m, original viscosity model.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 12.** Numerical non-dimensional TKE profile vs. experimental measurements for davg=0.7mm at X=11.5m, extended viscosity model.</figcaption></figure>

The black line with solid squares is the turbulent kinetic energy profile measured at X=11.5m for a layer of hard snow that does not drift. This was done in the experiment to assess the drifting particles effect on turbulent kinetic energy. The black line with x symbols is a measurement of the turbulent kinetic energy profile for drifting snow. Clearly there is production of turbulent kinetic energy by the drifting particles. Moreover, the experiment authors report that their ultra-sonic anemometer could not measure small-scale fluctuations at frequencies higher than 20Hz. It is then expected that the actual turbulent kinetic energy would be even higher than measured. However, the simulation results for both viscosity models are more than an order of magnitude larger than the experimental measurements so it is not likely that the small-scale component would make up the large differences seen in Figure 11 and Figure 12. The extended viscosity model has no effect on the airflow velocity and the turbulent kinetic energy profiles because the viscosity is part of the viscous stress tensor only in Equation (2) and does not affect the momentum coupling term which would impact the air momentum equation and the airflow velocity as a consequence.

## 4 Conclusions and future work

The extended viscosity model was found to provide a notable improvement in the suspension layer on the original one. This improvement provided a good match in the snow flux profiles across the saltation layer, and well into the suspension layer. The improvement is due to the modelling of the turbulent lift through the Gosman $C_t$ turbulent multiplier in a Boussinesq fashion.

There was no improvement in the airflow velocity profiles since the viscosity only affects the viscous stress tensor and not the momentum coupling term. The airflow velocity deficit is due to an excess production of the turbulent kinetic energy by the k-ε turbulence model, and there are different ways to remedy this deficiency in the k-ε formulation [17][18]. Other options include switching to the k-ω SST model which does a better job at predicting turbulent kinetic energy, or even the higher order DES and LES turbulence models. Work is already underway to this effect.

Another area of improvement is modifying the code to be able to handle several particle diameters simultaneously. This would eliminate the significant computational overhead required in computing several individual particle diameter simulations, and combining them in a post-processing exercise to get the particle size distribution effect. Handling several particle classes simultaneously can be done by solving a Probability Density Function transport equation, or solving the continuity and momentum balance equations of several diameter classes directly. The latter approach is being implemented and tested presently with three particle classes.
## Acknowledgments

The authors warmly thank Professors Akashi Mochida and Yoshihide Tominaga, as well as Dr Tsubasa Okaze for sharing their experimental results. In particular the patience and dedication of Dr Okaze to answering our numerous questions are gratefully acknowledged. Many thanks!

## References

1. B. Lee, J. Tu, C. Fletcher. On numerical modeling of particle-wall impaction in relation to erosion prediction: Eulerian versus Lagrangian method. *Wear* **252**.

2. T. Uematsu, T. Nakata, K. Takeuchi, Y. Arisawa, Y. Kaneda. Three-dimensional numerical simulation of snowdrift. *Cold Regions Science and Technology* **20**.

3. M. Naaim, F. Naaim-Bouvet, H. Martinez. Numerical simulation of drifting snow: erosion and deposition model. *Annals of Glaciology* **26**.

4. Y. Tominaga, T. Okaze, A. Mochida. CFD prediction of flowfield and snowdrift around a building complex in a snowy region. *Building and Environment* **46**.

5. J. Pomeroy, D. Gray. Saltation of snow. *Water Resources Research* **26**.

6. T. Okaze, A. Mochida, Y. Tominaga, M. Nemoto, Y. Ito, T. Shida. Modeling of drifting snow development in a boundary layer and its effect on windfield. *The Sixth Snow Engineering Conference, Whistler, B.C.*, Canada.

7. P. Gauer. Blowing and drifting snow in Alpine terrain: a physically-based numerical model and related field measurements. *PhD dissertation, ETH Zurich*.

8. P. Sundsbø. Numerical simulations of wind deflection fins to control snow accumulation in building steps. *Journal of Wind Engineering and Industrial Aerodynamics* **74-76**.

9. J.H.M. Beyers. Numerical modeling of the snowdrift characteristics surrounding the SANAE IV research station. *PhD dissertation, University of Stellenbosch*.

10. B. Bang, A. Nielsen, P. Sundsbo, T. Wiik. Computer simulation of wind speed, wind pressure and snow accumulation around buildings (SNOW-SIM). *Energy and Buildings* **21**.

11. Z. Boutanios, H. Jasak. Two-way coupled Eulerian-Eulerian simulations of drifting snow with viscous treatment of the snow phase. *Journal of Wind Engineering and Industrial Aerodynamics* (submitted for publication).

12. T. Okaze, A. Mochida, Y. Tominaga, M. Nemoto, T. Sato, Y. Sasaki, K. Ichinohe. Wind tunnel investigation of drifting snow development in a boundary layer. *Journal of Wind Engineering and Industrial Aerodynamics* **104-106**.

13. A.D. Gosman, R.I. Issa, C. Lekakou, M.K. Looney, S. Politis. Multidimensional modelling of turbulent two-phase flows in stirred vessels. *AIChE J., 38(12)*.

14. T. Okaze, A. Mochida, Y. Tominaga, M. Nemoto, T. Sato, Y. Sasaki, K. Ichinohe. Wind tunnel investigation of drifting snow development in a boundary layer. *J. Wind Eng. Ind. Aerodyn.*, **104-106**.

15. F. Budd. The drifting of nonuniform snow particles. *Stud. Antarctic Meteorol., Am. Geophys. Union, Antarctic Res. Ser.*, **9**.

16. R.A. Schmidt. Vertical profiles of wind speed, snow concentration, and humidity in blowing snow. *Boundary Layer Meteorology*, **23**.

17. P.A. Durbin. On the k-3 stagnation point anomaly. *Int. J. Heat and Fluid Flow*, **17**.

18. V. Guimet, D. Laurence. A linearised turbulent production in the k-epsilon model for engineering applications. *In Engineering Turbulence Modelling and Experiments*, **5**.
# CFD prediction of snowdrift in a building array

Tsubasa Okaze<fnref n="1" />,ᵃ, Saeka Kato<fnref n="2" />, Yoshihide Tominaga<fnref n="3" />, Akashi Mochida<fnref n="4" />

**Abstract.** This study aimed to investigate the performance of CFD prediction of snowdrift in and around a building array. The studied building array was located in Sapporo, Japan, and consisted of 12 one-story buildings built in three rows and four columns. The snowdrift model was based on two transport equations for two drifting snow densities corresponding to snow particles falling from the sky and those lifted up from the surface by strong wind. The predicted results were compared with field measurement results. The predicted snow depth distribution showed good agreement with measurement results, except for the region in front of windward buildings and the open space that allows flow in the prevailing wind direction.

## 1 Introduction

Approximately 60% of Japan consists of snowy regions; in such regions, snowdrift formation around buildings can considerably affect human activities. Computational fluid dynamics (CFD) is expected to be a useful tool for simulating snowdrift around buildings in such snowy regions.

Recently, Okaze et al. [1] proposed a new snowdrift model based on two transport equations of drifting snow densities for two different types of settling velocities corresponding to snow particles falling from the sky and those lifted up from the snow surface by strong winds. The proposed model was validated with detailed field measurements around a cubic building model in Sapporo, Japan by Oikawa et al. [2]. In almost all the previous studies, including those of the authors, the performance of the model was quantitatively validated only for snowdrift around an isolated building.

Snowdrift formation is more complicated for multiple buildings. Bang et al. [3] predicted snowdrift patterns (e.g., deposition/erosion) around three houses located in northern Norway by using a drifting flux model with two-phase flow equations. This [4] predicted snowdrift around two adjacent cubical buildings by using the model proposed by Bang et al. [3]. Beyers and Weacher [5] demonstrated the results of the application of snowdrift predictions for three elevated-houses. However, very few studies have dealt with the performance of the model for predicting snowdrift in and around multiple buildings because of the lack of detailed full-scale field measurements needed for validation. To take effective preventive measures against snowdrift in an actual urban environment, snowdrifts need to be predicted quantitatively around buildings under a complicated flow field. Recently, Tsutsumi et al. [6] conducted a detailed field measurement around twelve buildings by using a 3D laser scanner with a resolution of 0.5 m. This measurement data is useful for the validation of CFD models.

This study aims to investigate the performance of the snowdrift model by comparing CFD model results with

field measurements obtained by Tsutsumi et al. for snowdrift formation around a building array [6].

## 2 Outline of the snowdrift model

In this study, a CFD model based on two transport equations for two drifting snow densities [1] was applied to predict snowdrift formation around buildings. The drifting snow density was divided into two parts: the density of complex-shaped snow particles falling from the sky $\langle\Phi_{sky}\rangle$ and that of snow particles lifted up from the snow surface by strong winds $\langle\Phi_{surf}\rangle$, which are often broken apart by repeated collision with the snow surface. Thereafter, two transport equations for the two drifting snow densities, i.e., eqs.(1) and (2), were solved to predict the spatial distributions of both types of drifting snow densities with two different settling velocities.

$$ \frac{\partial \langle\Phi_{sky}\rangle}{\partial t} + \frac{\partial \langle\Phi_{sky}\rangle \langle u_i \rangle}{\partial x_i} + \frac{\partial \langle\Phi_{sky}\rangle \langle w_{surf} \rangle}{\partial x_3} = \frac{\partial}{\partial x_i} \left[ \frac{v_i}{\sigma_s} \left( \frac{\partial \langle\Phi_{sky}\rangle}{\partial x_i} \right) \right] \quad (1) $$

$$ \frac{\partial \langle\Phi_{surf}\rangle}{\partial t} + \frac{\partial \langle\Phi_{surf}\rangle \langle u_i \rangle}{\partial x_i} + \frac{\partial \langle\Phi_{surf}\rangle \langle w_{surf} \rangle}{\partial x_3} = \frac{\partial}{\partial x_i} \left[ \frac{v_i}{\sigma_s} \left( \frac{\partial \langle\Phi_{surf}\rangle}{\partial x_i} \right) \right] \quad (2) $$

To define the boundary conditions of the transport equations for $\langle\Phi_{sky}\rangle$ and $\langle\Phi_{surf}\rangle$ on the snow surface, the following assumptions were made: the snow particles falling from the sky break apart immediately upon collision with the snow surface, and no snow particles with their original complex shapes re-entrain from the surface. Based on these assumptions, the turbulent flux of $\langle\Phi_{sky}\rangle$ at the snow surface was set to zero.

$$ \frac{v_i}{\sigma_s} \left( \frac{\partial \langle\Phi_{sky}\rangle}{\partial x_3} \right) \Bigg|_{\text{surface}} = 0 \quad (3) $$

On the other hand, as mentioned earlier, all snow particles lifted up from the surface are assumed to immediately change to an approximately spherical shape owing to repeated collisions with the snow surface. The boundary condition for $\langle\Phi_{surf}\rangle$ is given by eq.(4) by using the erosion flux on surface $E_{surf}$ kg/(m<fnref n="2" />·s). Here, erosion occurs when the friction velocity on the surface $\langle u^* \rangle$ is larger than the threshold friction velocity ($\langle u^* \rangle > \langle u^*_{i} \rangle$).

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

<footnote n="1">¹ Tokyo Institute of Technology, 4259 Nagatsuta-cho, Midori-ku, Yokohama, Japan</footnote>

<footnote n="2">² Takenaka Corporation, 4-1-13 Hommachi, Chuo-ku, Osaka, Japan</footnote>

<footnote n="3">³ Niigata Institute of Technology, 1719 Fujihashi, Kashiwazaki, Japan</footnote>

<footnote n="4">⁴ Tohoku University, 6-6-11 Aoba, Aramaki Aza, Aoba-ku, Sendai, Japan</footnote>

<footnote>ª Corresponding author: okaze.t.aa@m.titech.ac.jp</footnote>
The relationship between friction velocity and erosion flux for blowing sand as proposed by Shao and Li [7] was used to model the erosion flux on the snow surface attributed to the shear stress exerted by wind.

$$E_{surf} = -\frac{\pi \zeta}{6} \rho_t \langle u^* \rangle \left( 1 - \frac{\langle u^* \rangle_i^2}{\langle u^* \rangle} \right) \Delta x \Delta y \quad \text{for} \quad \langle u^* \rangle > \langle u^* \rangle_i \quad (4)$$

The model coefficient $\zeta$ is set to $1.0 \times 10^{-5}$. When erosion does not occur ($\langle u^* \rangle \le \langle u^* \rangle_i$), a zero gradient condition is adopted for $\langle \Phi_{surf} \rangle$ because $E_{surf} = 0$ in such a situation. Snow deposition per unit time and per unit area is expressed as $\langle \Phi \rangle_p \langle w f \rangle$. Here, the subscript “p” represents the value at the first grid point. The deposition fluxes $D_{sky}$ and $D_{surf}$ kg/(m²·s) corresponding to $\langle \Phi_{sky} \rangle$ and $\langle \Phi_{surf} \rangle$ are given as follows:

$$D_{sky} = -\langle \Phi_{sky} \rangle_p \langle w f_{sky} \rangle \Delta \lambda \quad (5)$$

$$D_{surf} = -\langle \Phi_{surf} \rangle_p \langle w f_{surf} \rangle \Delta \lambda \quad (6)$$

Finally, the net deposition rate $M_{total}$ kg/(m²·s) is given as shown below:

$$M_{total} = D_{sky} + D_{surf} + E_{surf} \quad (7)$$

The change in snow depth per unit time $\Delta z_s$ m/s is obtained from $M_{total}$ divided by the accumulated snow density $\rho_s$ kg/m³:

$$\Delta z_s = \frac{M_{total}}{\rho_s} \quad (8)$$

Then, the change in snow depth during snowfall period was calculated by multiplying the variation of snow depth per unit time by the time during snowfall period. The period was set to an hour in each simulation as mentioned later in chapter 3. In this study, the influence of surface change on the flow field was not considered.

# 3 Prediction of snowdrift around buildings

## 3.1 Computational conditions

A field measurement was conducted on the campus of Hokkaido University of Science in Sapporo, Japan, from December 2009 to February 2010 by Tsutsumi et al. [6]. The building array consisted of 12 one-story prefabricated buildings built in three rows and four columns with intervals of 2.7 m. The dimensions of each building were 2.4 m (D) × 5.4 m (W) × 2.7 m (H) as shown in Figure 1. The depth of snow around the buildings was measured using a 3D laser scanner with a resolution of 0.5 m. Although the snow depths were measured on January 7 and February 8, 2010, we selected the snowdrift event that occurred from December 12, 2009 to January 7, 2010 as the first test case to investigate the performance of the prediction model. The meteorological data used as boundary conditions were obtained from the observational data of Yamaguchi meteorological station, which is the station nearest to the measurement site. Figure 2 shows a wind rose diagram obtained from Yamaguchi meteorological station during the measurement period. The prevailing wind was west-northwest. Figure 3 shows time histories of the measured

hourly averaged data at Yamaguchi station during the measurement period (Dec. 12, 2009–Jan. 7, 2010). The wind velocity was measured at a height of 10 m. It was assumed that snowfall occurred when precipitation was observed and the temperature did not exceed 2.0°C. The steady-state change of snow depth per unit time given by eq.(8) was computed with the hourly averaged meteorological data when snowfall was considered to occur. Then, the hourly change in snow depth was calculated by multiplying the variation in snow depth per unit time by 3600 s (= 1 h). Finally, the total change in snow depth was estimated by summing the hourly change in snow depth during snowdrift events. Inflow and top boundary conditions of $\langle \Phi_{sky} \rangle$ were given by using eq.(9):

$$\langle \Phi_{sky} \rangle = \frac{\rho_s h^*}{\langle w f_{sky} \rangle \Delta t_s} \quad (9)$$

Here, $h^*$ m is the reference snow depth defined as a snow depth at each snowfall for a period $\Delta t_s$[s] (= 3600 s) in regions unaffected by buildings. The reference snow depth was estimated by the following relation:

$$h^* = \frac{P \rho_w}{\rho_s} \quad (10)$$

where, $P$, $\rho_w$ and $\rho_s$ mean hourly measured precipitation

<figure><img src="image_37.png" /><figcaption>Figure 1. Schematic view of field measurement [6]</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 2. Wind rose at Yamaguchi station (Dec. 12, 2009–Jan. 7, 2010)</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>(1) Wind velocity</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(2) Precipitation and temperature</figcaption></figure>

<caption>**Figure 3.** Time histories of measured hourly averaged data at Yamaguchi station (Dec. 12, 2009–Jan. 7, 2010)</caption>

<table><caption>**Table 1.** Computational conditions</caption><tr><td>Grid discretization</td><td>50(x₁) × 42(x₂) × 21(x₃)</td></tr><tr><td>Inlet boundary</td><td>⟨u₁⟩, ⟨u₂⟩: z₀-type logarithmic law<br>(z₀ = 0.41 m which corresponds to the roughness length surrounding low-rise houses)<br>⟨u₃⟩ = 0<br>k = ⟨u*⟩² / C<sub>μ</sub><sup>0.5</sup><br>ε: assuming the local equilibrium of production of k and energy dissipation rate (P<sub>k</sub> = ε)</td></tr><tr><td>Outlet boundary</td><td>Zero gradient condition</td></tr><tr><td>Ground boundary</td><td>z₀-type logarithmic law<br>(z₀ = 1.0 × 10⁻³ m: roughness length on snow surface)</td></tr><tr><td>Building wall boundary</td><td>Generalized logarithmic law</td></tr><tr><td>Lateral and upper boundaries</td><td>Slip wall condition</td></tr><tr><td>Spatial derivatives for convection term</td><td>QUICK scheme</td></tr><tr><td>Turbulence model</td><td>Revised k-ε model<br>proposed by Durbin [8]</td></tr></table>

at Yamaguchi meteorological station, a density of water and a density of accumulated snow, respectively. The density of accumulated snow was assumed to 100 kg/m³ as that of fresh snow. The vertical profile of mean wind velocity at inflow boundary obeys the z₀-type logarithmic law. The roughness length z₀ was set to 0.41 m corresponding to surrounding low-rise houses around the measurement site. Other computational conditions for CFD are listed in Table 1. These conditions basically follow guidelines for practical applications of CFD to pedestrian wind environment proposed by the Architectural Institute of Japan [9]. The revised k-ε model

<table><caption>**Table 2.** Parameters of snow</caption><tr><td>Snowfall velocity<br>corresponds to ⟨Φ<sub>sky</sub>⟩: ⟨wf<sub>sky</sub>⟩</td><td>− 1.0 m/s</td></tr><tr><td>Snowfall velocity<br>corresponds to ⟨Φ<sub>surf</sub>⟩: ⟨wf<sub>surf</sub>⟩</td><td>− 0.2 m/s</td></tr><tr><td>Threshold friction velocity: ⟨u*⟩<sub>t</sub></td><td>0.20 m/s</td></tr><tr><td>Density of accumulated snow: ρ<sub>s</sub></td><td>100 kg/m³</td></tr></table>

proposed by Durbin [8] was employed as the base model to calculate eddy viscosity. Parameters related to the property of snow are summarized in Table 2. For the purpose of this study, it was assumed that other factors related to changes in snow depth, such as snowmelt, were negligible.

## 3.2 Results and discussion

The snow depth was normalized by the reference snow depth h* that was not affected by buildings.

Figure 4 shows the distribution of the measured and predicted snow depths on January 7. The eroded regions between the longer sides of the buildings as observed in the CFD results correspond to the regions as determined from the measurement. However, the large snow depth measured in front of the windward building was not reproduced by CFD. One reason for this disagreement can be that an hourly averaged wind velocity was applied to the computation. Hence, although large erosion usually occurs when the instantaneous wind velocity is greater than the threshold velocity in real situations, this effect is not included in this CFD model.

Figure 5 compares the longitudinal distributions of snow depths along sections A-A' and B-B' shown in Figure 1. In section A-A', which represents the space obstructed by buildings, the distribution of snow depth obtained by CFD shows good agreement with that obtained from the measurement, except for the region in front of windward buildings. The observed erosion between buildings and the accumulation of snow near the
building wall are well reproduced. In section B-B', which represents the open space allowing flow in the prevailing wind direction, the CFD generally underestimated snow depth, and peak values observed in the measurement results were not well reproduced by CFD.

This underestimation is also attributed to the fact that this simulation did not consider the temporal high wind speed velocity, as mentioned earlier. Therefore, the amount of snow directly transported from the windward side of buildings to the region between buildings was underestimated.

## 4 Conclusions

Snowdrift in and around building arrays was predicted using CFD. The predicted distribution of snow depth agreed well with the measured snow depth values, except for the regions in front of windward buildings and the open space allowing flow in the prevailing wind direction. The effect of temporal change in wind velocity on the accuracy of steady-state simulation will be examined in future studies.

## References

1. Okaze T., Mochida A., Tominaga Y., Ito Y., Yoshino H., CFD prediction of snowdrift around a cube using two transport equations for drifting snow density, Proceedings of the 5th International Symposium on Computational Wind Engineering (2010)

2. Oikawa S., Tomabechi T., Daily observation of snowdrifts around a model building. Proceedings of the 4th International Conference on Snow Engineering, 149-156 (2000)

3. Bang B., Nielsen A., Sundsbø P. A., Wilk T., Computer simulation of wind speed, wind pressure and snow accumulation around buildings (SNOW-SIM). Energy and Buildings, 21, 235-243 (1994)

4. Thiis K. T., Large scale studies of development of snowdrifts around buildings. Journal of Wind Engineering and Industrial Aerodynamics, 91, 829-839 (2003)

5. Beyers M., Waechter B., Modeling transient snowdrift development around complex three-dimensional structures. Journal of Wind Engineering and Industrial Aerodynamics, 96, 1603-1615 (2008)

6. Tsutsumi T., Chiba T., Tomabechi T., Snowdrefts on and around buildings based on field measurement, Proceedings of the 7th International Conference on Snow Engineering, 9-17 (2012)

7. Shao Y., Li A., Numerical modeling of saltation in the atmospheric surface layer, Boundary-Layer Meteorology, 91, 199-225 (1999)

8. Durbin, P.A., On the k-3 stagnation point anomaly. Int. J. Heat and Fluid Flow, 17, 89-90 (1996)

9. Tominaga, Y., Mochida, A., Yoshie, R., Kataoka, H., Nozu, T., Yoshikawa, M., Shirasawa T., AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings. Journal of Wind Engineering and Industrial Aerodynamics, 96(10-11), 1749-1761 (2008)

<figure><img src="image_41.png" /><figcaption>Figure 4. Horizontal distribution of snow depth on January 7, 2010</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Figure 5. Longitudinal distribution of snow depth on January 7, 2010</figcaption></figure>
# Outdoors experiments of snowdrift on typical cubes based on axial flow fan matrix in Harbin

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

**Abstract.** Snowdrifts around different model buildings were investigated in wind tunnel and outdoor environment, considering the serious problems caused by snowdrifts around buildings when formed on undesirable places. The formation of snowdrifts is highly connected to the wind speed around the buildings. Some experiments have been carried out to find the relationship between the snow deposition around the building and wind speed. The experiments were conducted with the same model, 1000×1000×1000 mm³ but under three different wind speeds. All of the experiments were based on a new type of laboratory equipment - axial flow fan matrix, which was placed outdoor. The laboratory equipment can produce stable wind field within a certain range and a speed range of 1~10.5 m/s. The laboratory equipment was placed in a narrow strip between two buildings in Harbin China. The temperature around was measured during the experiments. The snowdrift height around the model was surveyed in approximately 1000 points and contour maps produced. The results proved that different wind speeds had a significant influence on the pattern of the snowdrifts produced.

## 1 Introduction

Snowdrift formation around buildings should be considered areas prone to be snowy and windy to ensure proper accessibility and safety of buildings. Most of the researchers use numerical simulation to model snow accumulation around buildings, while this tool needs to be carefully validated before it can be used for planning purposes. Important understanding for use in such validation is gathered in outdoor full-scale experiments and wind tunnels experiment. Full-scale experiment is sparsely reported. As an appropriate profiler to simulate snowdrift problems around buildings, wind tunnel experiment methods have been explored for many years. Related research has been undertaken in several studies.

Naaim et al. [1] reviewed the conditions for a correct simulation of snowdrifts in wind tunnels, or outdoor environment, by comparing requirements proposed by four scientists (i.e., Kind, 1975, 1985; Iversen, 1979, 1980, 1981, 1982, 1984; Anno, 1984; Tabler, 1980). Delpech [2] 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. [3] carried out a series of wind tunnel measurements over a loose snow surface to investigate the characteristics of drifting snow phenomenon 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 above, 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 [4] conducted a series of experiments in wind tunnel to investigate the motion of natural snow. The results indicated that the mass flux rates of snow exponentially decreased with

height whereas the snow transport rate increased exponentially with wind velocity.

Although various studies dealing with snowdrift simulation have been carried out in recent years, due to the limitation of the experimental instrumentation and funds, such experiments suffer a lot from scaling problems and artificial snow [5, 6]. Compromises must be made based on the objective of the experiment.

To solve the problems, Harbin Institute of Technology (HIT) proposed a new kind of experimental equipment, axial flow fan matrix, which aimed at the experiments of snow drifting. The experimental equipment was composed of a 2×2 fan matrix, wind rectifier, and a snow simulator. Considering the device is sited outdoors, the experimental equipment was placed in a narrow strip between two buildings to reduce the interference of the external environment. In order to meet the experimental requirements, the laboratory equipment could produce stable wind field within a certain range and a speed range of 1~11.5 m/s.

Treating the building as a cube with flat roof, several experiments were conducted to find the relationship between wind speeds and the pattern of the snow accumulation around buildings.

## 2 Outdoors experiments of snowdrifts around cube model

### 2.1 Experiment condition and method

The experiment was performed within the campus of Harbin Institute of Technology (HIT), China (Figure 1). It was a part of a large snow engineering program. The natural wind speed in this area has been previously found to be very small during winter. This condition makes the area ideal for snowdrift measurements.

<figure><img src="image_43.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_44.png" /><figcaption>**Figure 1.** Locations of experiments</figcaption></figure>

To find the relationship between wind speeds and the pattern of the snow accumulation around buildings, three experiments using the same model but under three different wind speeds-2.5m/s, 3.5m/s, 4.5m/s, were conducted. The model was of cube shape with dimensions of 1000 × 1000 × 1000 mm3.

The experiments were conducted on February 1st, 2016, based on the axial flow fan matrix (Figure 2). The snow used during the experiments comes from real snow. During the experimental process, snow density and natural temperature were measured. Natural wind speed and direction were measured at 1.8m on a PC-4 automatic weather station (Figure 3), and the sampling interval was 1 minute.

<figure><img src="image_45.png" /><figcaption>**Figure 2.** Experiment equipment</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Figure 3.** PC-4 automatic weather station</figcaption></figure>

## 2.2 Results

Figure 5 shows the snowdrift geometry measured at afternoon on February 1st, 2016. The natural temperature is about -13°C.

During the experiments, the natural wind in the experiment site was measured (Figure 4).

<figure><img src="image_47.png" /><figcaption>**Figure 4.** Wind frequency diagram in the winter of 2015-2016</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>**Figure 5(a).** Photos of snowdrift with 2.5m/s wind speed</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>**Figure 5(b).** Photos of snowdrift with 3.5m/s wind speed</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>**Figure 5(c).** Photos of snowdrift with 4.5m/s wind speed</figcaption></figure>

The incidence angle of the wind was 0° on the windward wall. During the experiments, the model was fixed on the experimental platform with a uniform flat bottom. The distance between the model and the outlet surface of wind was 2m.

Every experiment would last for 10min. At the end of the experiment, the height of the snowdrifts was measured with a straight edge tool in approximately 1000
points around each building. Because the data were measured manually, the measurement accuracy was down to millimeter level. The points were chosen for optimal description of the curvature of the snowdrifts. The Kriging method was used to interpolate a surface between the measured points. Figure 6 showed the contour map of snowdrift.

<figure><img src="image_51.png" /><figcaption>Figure 6 (a). Contour map of snowdrift with 2.5m/s wind speed</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Figure 6 (b). Contour map of snowdrift with 3.5m/s wind speed</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>Figure 6 (c). Contour map of snowdrift with 4.5m/s wind speed</figcaption></figure>

## 3 Comparison and analysis

The temperature during the experiment ranged between -10°C and -13°C. The 1 min average natural wind speed in the period was ranging from 2.5 to 4.5 m/s.

The snowdrifts around the same model building, presented in Fig. 6, showed more or less the same basic characteristics. The accumulation pattern proved the close relationship between wind speed and snowdrift around buildings.

The results of the three experiments show some identical characteristics. The snowdrifts showed a symmetrical distribution, similar to the wind-induced vortices. Just close to the upwind, there was a short steep inclining snowdrift, but the range of the snowdrift was short. There was a sharp increasing of snowdrift in an area, which was from a certain distance to the upwind surface. With the distance away from the upwind surface becoming larger, the accumulation effect weakened. In the leeward side of the model, there was less but uniformed snow accumulation.

However, because the simulated wind speed was different, each experiment result showed unique characteristics. Comparing the three snow distribution contour maps we could conclude that:

Figure 6 (a) shows the results of the experiment under the wind speed 2.5 m/s. The height of snowdrift close to the upwind surface was the highest and the range of snowdrift was the largest. At about 0.6m away from the upwind surface, there came an area of obvious snow deposition and a maximum height of snowdrift. With the increase of distance away from upwind surface, the height of snowdrift decreased. Although there was no snowdrift near the model on both sides, at about 0.6m away from the upwind surface, there was some small snowdrift. At about 1.2 m away from downwind, there were two obvious snow depositions on both sides. The whole snowdrifts showed a typical horseshoe pattern.
Figure 6 (b) shows the results of the experiment under the wind speed 3.5 m/s. The height of snowdrift close to the upwind surface was slightly lower and the range of snowdrift was slightly smaller. At about 0.75m away from the upwind surface, there was a steep inclining snowdrift ending in a sharp edge and a maximum height of snowdrift, indicating the extension of the upwind recirculation zone. With the increase of distance away from upwind surface, the height of snowdrift decreased. At about 0.8m away from the upwind surface, there was some small snowdrift. There appeared an area with no snowdrift in the range of 0.5m~1m away from the downwind surface and no obvious area of snow deposition appeared on both sides.

Figure 6 (c) shows the results of the experiment under the wind speed 4.5 m/s. The height of snowdrift close to the upwind surface was the lowest and the range of snowdrift was smallest. At about 1.0m away from the upwind surface, there was a most obvious steep inclining snowdrift ending in a sharp edge and a maximum height of snowdrift, indicating the extension of the upwind recirculation zone. In the measuring range, there was no snowdrift on both sides of the model. A typical horseshoe pattern of snowdrift appeared at the area close to the downwind surface.

## 4 Conclusions

This paper has investigated the characteristics of snowdrift in outdoor experiments based on axial flow fan matrix. The results proved that different wind speeds had significant influences on the pattern of the snowdrifts produced. As the wind speed increased, the snowdrift just close to the upwind surface decreased rapidly, the snowdrift away from the upwind surface became stronger and steeper and the distance between the sharp snowdrift and the upwind surface became larger. Although the sphere of influence was not measured, we could predict that the range of wind influence on snowdrift became larger as the wind speed increased.

## References

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

2. Delpech Ph, Palier P, Gandemer J. Snowdrifting Simulation around Antarctic Buildings[J]. Journal of Wind Engineering and Industrial Aerodynamics, 74-76, 567-576 (1998)

3. Tsubasa Okaze, Akashi Mochida, Yoshihide Tominaga. Wind tunnel investigation of drifting snow development in a boundary layer[J]. Journal of Wind Engineering and Industrial Aerodynamics, 104-106, 532-539(2012)

4. LüXiaoHui, Huang N, Tong D. Wind tunnel experiments on natural snow drift[J]. Science China Technology Science, 55, 927-938 (2012)

5. Thomas K. Thiis, Yngvar Gjessing. Large-scale measurements of snowdrifts around flat-roofed and single-pitch-roofed buildings[J]. Cold Regions Science and Technology 30. 175-181(1999)

6. J.H.M. Beyers, T.M. Harms. Outdoors modelling of snowdrift at SANAE IV Research Station, Antarctica[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2002, 91 : 51-569. (2004)

7. A. Mecke, I. Lee, J.R. Baker jr., M.M. Banaszak Holl, B.G. Orr, Eur. Phys. J. E **14**, 7 (2004)

8. M. Ben Rabha, M.F. Boujmil, M. Saadoun, B. Bessaïs, Eur. Phys. J. Appl. Phys. (to be published)

9. Luigi T.De Luca, *Propulsion physics* (EDP Sciences, Les Ulis, 2009)

10. F. De Lillo, F. Cecconi, G. Lacorata, A. Vulpiani, EPL, **84** (2008)
# Numerical simulation of snowdrift on a membrane roof and wind-induced response analysis under coupled wind and snow loads

SUN Xiaoying¹, HE Rijin¹, WU Yue¹

¹Key Lab of Structure Dynamic Behaviour and Control of the Ministry of Education, Harbin Institute of Technology, 150090 Harbin, China

**Abstract.** Long-span structures are quite sensitive to uneven distributed snow loads, therefore it is vital to estimate the wind-induced redistribution snow loads accurately. This paper presents the results of numerical simulation of snowdrift on a long-span membrane roof at first. Based on Euler-Euler method in multi-phase flow theory, the numerical simulation adopt Mixture model by employing commercial CFD software FLUENT. Assuming the snow is evenly distributed at the initial time in simulation, and the initial snow depth is calculated by the referenced snow pressure in Chinese Load Code, combined with the snow deposition and erosion model, the result is given as the counter map of snow distribution coefficient under different wind velocities and directions, and we give the most unfavourable form of snow distribution. Next, the uneven distributed snow is set as static load on the membrane roof, and the unsteady aerodynamic force is obtained by numerical simulation either, then the wind-induced response under coupled wind and snow load is analysed, the results show that uneven distributed snow is a disadvantage mechanical condition to this structure. Finally, under the condition of laying air-blowing pipelines on the roof, the snow sweeping process is simulated by employing transient solver in FLUENT, the result is given as the time dependent residual snow mass fraction and snow mass flux which provide some useful suggestions on practical application.

## 1 Introduction

Wind causes substantial snow accumulation and unevenly redistributed on and around buildings which often leads to structural damage or even collapse. Long-span structures are quite sensitive to unevenly distributed snow loads, therefore it is vital to estimate the wind-induced redistribution snow loads accurately.

Transportations of snow particles can be classified into three process i.e. creep, saltation and suspension ([1]). Creep is a process in which snow particles moves by rolling or creeping on the snow surface under low wind speed, the height of creep is usually below 0.01m. Saltation is a process in which snow particles move with repeat leaping up or jumping on the snow surface under relatively high wind speed, the height of saltation is usually below 0.1m. If wind speed continues to increase, the snow particles are transported upwards by turbulent eddies and moved far downwind, this process is called suspension, the floating height of snow particles in suspension is about 0.1~100m. The snow transport rate in creep is so small that can be neglected, as well as previous research, the snow movement in creep is not considered.

Currently, there are mainly three methods to estimate snowdrift namely field observation, wind tunnel test and numerical simulation. Field observation is the most realistic approach to obtain the uneven snow distribution, but this approach requires precise measurements for snow scale and meteorology conditions and relies heavily on the meteorology conditions. Wind tunnel test constrained greatly by the similarity criterion and cost a lot. Numerical simulation was firstly adapted to analysis snowdrift by Uematsu at 1991 ([2]) and now widely used by researchers.

## 2 Numerical modelling

Based on one way coupled assumption, the snow particle effects on the turbulence are not considered, the simulated turbulence is influenced directly by the snow phase through the wall function.

Based on Euler-Euler approach in multiphase flow theory, there are three multiphase models supported by FLUENT, i.e. the VOF model, the Mixture model and the Eulerian model. Compared these three models, Mixture model allowing better stability and smaller computation labour, while sharing nearly the same precision with Euler model. For this reason, this paper adopt Mixture model to simulate snowdrift.

### 2.1 Governing Equations ([3])

In Mixture Model, the air and snow shares the same mass continuity equation and momentum conservation equation as following formulas.

#### 2.1.1 Mass continuity equation

$$ \frac{\partial \rho_m}{\partial t} + \nabla(\rho_m \bar{v}_m) = 0 \quad (1) $$

Where $\rho_m = \sum_{k=1}^{2} f_k \rho_k$ is the mixture density of air and snow, $f_k$ is the volume fraction of the air ($k=1$) and snow ($k=2$), $\bar{v}_m = \sum_{k=1}^{2} (f_k \rho_k \bar{v}_k) / \rho_m$ is velocity of the mixture phase.

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

<footnote>Corresponding author: sxy_hit@163.com</footnote>
### 2.1.2 Momentum Conservation Equation

$$ \frac{\partial(\rho_m \bar{v}_m)}{\partial t} + \nabla \cdot (\rho_m \bar{v}_m \bar{v}_m) = -\nabla p + \nabla \cdot \left[ \mu_m (\nabla \cdot \bar{v}_m + \nabla \cdot \bar{v}_m^T) \right] + \rho_m \bar{g} + \bar{F} + \nabla \cdot \left( \sum_{k=1}^{2} f_k \rho_k \bar{v}_{dr,k} \bar{v}_{dr,k} \right) \quad (2) $$

Where $\mu_m = \sum_{k=1}^{2} f_k \rho_k$ is the mixture dynamic viscosity, $\bar{F}$ is the body force, $\bar{v}_{dr,k} = \bar{v}_k - \bar{v}_m$ is the drift velocity for phase k.

### 2.1.3 Snow Volume Fraction Equation

$$ \frac{\partial(\rho_s f)}{\partial t} + \nabla \cdot (\rho_s f \bar{v}_m) = -\nabla \cdot (f \rho_s \bar{v}_{dr,s}) + \sum (m_{sa} - m_{as}) \quad (3) $$

Where $\rho_s$ and $f$ are the snow density and volume fraction respectively, $m_{sa}$ and $m_{as}$ are the mass converted between air and snow, in this paper $m_{sa}$ and $m_{as}$ refers to 0.

## 2.2 Deposition and erosion model

Wall shear stress is a very important parameter in snowdrift, the snow particle begins to move when the wall shear stress $\tau$ on the snow surface exceed the threshold wall shear stress $\tau_t$, otherwise the snow particle remains stationary. The threshold wall shear stress $\tau_t$ is usually determined through field observation and the wall shear stress can be derived from the logarithmic law. Here and after, we use friction velocity $u_*$ to replace the wall shear stress $\tau$ and obeys the following relation:

$$ u_* = \sqrt{\tau/\rho} = \frac{\kappa u(z)}{\ln(z/z_s)} \quad (4) $$

Where $\kappa$ refers to Von Karman constant which equal to 0.4, $z_s$ is the aerodynamic roughness height over a snow surface which is optimized to $3.0 \times 10^{-5} m$.

The snow mass flux of erosion $q_{ero}$ is simplified as a function of the excess friction velocity available above the threshold friction velocity. The deposition flux $q_{dep}$ is calculated to be proportional to snow concentration $f$ and the snow fall velocity $w_f$.

$$ q_s = \begin{cases} q_{ero} = A_{ero}(u_*^2 - u_{st}^2) & u_* > u_{st} \\ q_{dep} = \rho_s f w_f \frac{u_{st}^2 - u_*^2}{u_{st}^2} & u_* \le u_{st} \end{cases} \quad (5) $$

Where $A_{ero}$ is a proportionality coefficient representing the snow pack bounding strength which is chosen to be $7.0 \times 10^{-4}$ ([4]), $u_{st}$ is the threshold friction velocity which is chosen to be 0.2 m/s ([1]), $w_f$ is the snowfall velocity of snow particles which is chosen to be 0.2 m/s.

If the initial snow depth is assumed to be $h_0$ and the time of wind acting on snow surface is $\Delta t$, then the depth of snow can be represented as:

$$ h = h_0 + \frac{q_s}{\rho_s} \Delta t \quad (6) $$

## 2.3 Turbulence Model

The turbulence model employed in numerical simulation of snowdrift varies from scholars. Uematsu ([2]) adopt the mixing length theory in his pioneering works on CFD prediction of snowdrift, Beyers (2004, [5]) adopt the standard $k-\epsilon$ model and points out that it is necessary to critically examine the influence of inaccuracies isotropic turbulence assumption on the snow accumulation and erosion, ZHOU and LI ([6]) adopt the modified RSM model which take account of the snow particle's effect on turbulence. HONG ([7]) compared the results of different turbulence models with field observation and optimized $k\bar{d}l\omega$ model for snowdrift simulation.

In this paper, the $k\bar{d}l\omega$ model is adopted. The $k\bar{d}l\omega$ model is considered to be a three-equation eddy-viscosity type, which includes transport equations for turbulent kinetic energy ($k_T$), laminar kinetic energy ($k_L$), and the inverse turbulent time scale ($\omega$).

$$ \frac{Dk_T}{Dt} = P_{k_T} + R + R_{NAT} - wk_T - D_T + \frac{\partial}{\partial x_j} \left[ \left( v + \frac{\alpha_T}{\alpha_K} \right) \frac{\partial k_T}{\partial x_j} \right] \quad (7) $$

$$ \frac{Dk_L}{Dt} = P_{k_L} - R - R_{NAT} - D_L + \frac{\partial}{\partial x_j} \left[ \left( v + \frac{\alpha_L}{\alpha_K} \right) \frac{\partial k_L}{\partial x_j} \right] \quad (8) $$

$$ \frac{D\omega}{Dt} = C_{\omega 1} \frac{\omega}{k_T} P_{k_T} + \left( \frac{C_{\omega 2} R}{f_{\omega}} - 1 \right) \left( R + R_{NAT} \right) \frac{\omega}{k_T} - C_{\omega 2} \omega^2 + C_{\omega 3} f_{\omega} \alpha_T f_W^2 \frac{\sqrt{k_T}}{d^3} + \frac{\partial}{\partial x_j} \left[ \left( v + \frac{\alpha_T}{\alpha_K} \right) \frac{\partial k_T}{\partial x_j} \right] \quad (9) $$

## 2.4 Boundary conditions

### 2.4.1 Wind profile

The fully developed inlet wind velocity profile is determined from the power law profile:

$$ \bar{u}(z) = u_0 \left( \frac{z}{z_0} \right)^{\alpha} \quad (10) $$

The reference velocity $u_0$ and the height above the surface $z_0$ are 10 m/s and 10 m respectively. Note that the friction velocity $u_*$ at inflow boundary should exceed the threshold friction velocity $u_{st}$ so that the snowdrift could initiate.
### 2.4.2 Snow particle concentration profile

The inlet snow particle concentration profile is given by an empirical formula proposed by Pomeroy, Gray & Male ([8], [9]), in which the inflow snow concentration is set as constant within saltation layer ($z \le h_{sal}$) and decreases exponentially with the increase of height within suspension layer ($z > h_{sal}$).

$$f = \begin{cases} \frac{0.68}{u_p u_{sal} g} u_{st} (u_s^2 - u_{st}^2) & z \le h_{sal} \\ 0.8 \exp[-1.55(4.78u_s^{-0.544} - z^{-0.544})] & z > h_{sal} \end{cases} \quad (11)$$

Where $u_p$ is the average velocity of the snow particle in saltation layer which is chosen to be 2.8 $u_{st}$ suggested by Pomeroy ([8]), $h_{sal}$ is the height of the saltation layer which is determined by the following formula ([8]):

$$h_{sal} = 1.6 \frac{u_s^2}{2g} \quad (12)$$

## 3 Numerical simulation of snowdrift

### 3.1. Geometrical modelling

The membrane structure (Figure 1) is located in Changchun, northeast China, a famous heavy-snowfall region. This mast supported tension structure spans about 188m, the highest point of the membrane part is about 18m and the height of masts is about 48m. For simplicity, the masts and surrounding landscape are neglected in CFD modelling. Due to its Hyperbolic Paraboloids building shape, snow stacks easily in the central concave region, therefore the snow load is one of the control loads in structural design.

<figure><img src="image_55.png" /><figcaption>Figure 1. Long-span membrane structure.</figcaption></figure>

Based the symmetry properties of this structure, four wind directions (i.e. 0°, 30°, 60°, 90°) are considered to estimate the effect of wind direction on snow distribution. The dimension of the computational region was selected as X×Y×Z=10L×21B×7H which is large enough to neglect the influence of the computational region boundary on interior flow field, where L, B, H refers to the across-wind dimension, along-wind dimension and the height of the structure respectively.

<figure><img src="image_56.png" /><figcaption>Figure 2. Sketch of model and computational region.</figcaption></figure>

### 3.2 Results

The following results are given as counters of the Snow distribution factor $\mu_r$ which is defined as:

$$\mu_r = \frac{S_k}{S_0} \quad (13)$$

Where $S_k$ refers to standard value of snow loads which can be calculated as $S_k = \gamma_{snow} h_{snow}$, $\gamma_{snow}$ is the snow density which is chosen as 1500 $N/m^3$ by measured data, $h_{snow}$ is calculated by Eq. (6). $S_0$ is the reference snow pressure which is 0.45 $kN/m^3$ according to Chinese load code for the design of building structures (GB 5009-2012, [10]).

<figure><img src="image_57.png" /><figcaption>Figure 3. Snow distribution factor under 0° wind direction</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>Figure 4. Snow distribution factor under 30° wind direction</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>Figure 5. Snow distribution factor under 60° wind direction</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>Figure 6. Snow distribution factor under 90° wind direction</figcaption></figure>
Under 0° wind direction, the whole membrane roof is represented as erosion region. In the middle part of the membrane roof, there is nearly no snow exists ($\mu_r$ =0) because of high wind velocity, but there is a strip area in the windward side where snow accumulates ($\mu_r$ =0.1~0.7) because of the flow separation. In both sides of the membrane roof, snow accumulates due to low wind speed, and snow accumulates heavier in leeward side of the membrane roof than that in windward side, this may attributes to the snow particles in windward side were blown away and deposits in leeward side of the roof.

Under 30° and 60° wind direction, the whole membrane roof is represented as erosion region too, and their snow distribution patterns are very similar with each other. Snow mainly accumulates in the left part of the membrane roof and the 60° case accumulates much more snow than that of 30° case, this non-symmetrical load may bring about disadvantage effect on the membrane structure.

Under 90° wind direction, the whole membrane roof is mainly represented as deposition region which is very different from other three cases. Most part of the roof has a snow distribution factor $\mu_r$ > 1 and the maximum snow distribution factor occurs in the leeward part excess 1.2. There are two major reasons: one is the profile of the membrane roof is very close to streamline in 90° case, this results in week flow separation and relatively low wall shear stress on the snow surface, so most part of the roof has a snow distribution factor $\mu_r$ > 1. The other reason is that the structure has a dimension of 188m along windward direction, it is long enough that the snow particles blown away in the windward part deposits substantially in the leeward part.

Comprehensively analysing these four cases, the membrane roof is divided into 64 zones (Figure 7) to give the most disadvantageous snow distribution factor of each partition for design use, the partition scheme is in accordance with the PTFE and ETFE zone. Due to the symmetry, only 16 zones were presented to give the most disadvantageous snow distribution factor.

<figure><img src="image_61.png" /><figcaption>**Figure 7.** Sketch of partition schemes of the membrane roof</figcaption></figure>

Table 1 presents the most disadvantageous snow distribution factors, it is the maximum value of snow distribution factor in each zone among these four wind direction cases. From table 1 we can find that the most disadvantageous snow distribution factors in concave area (zone 10~15) excess 1.0 entirely, the maximum value come up to 1.23. Note that the results in this paper is obtained under the condition of the wind velocity at 10m height is 10m/s, the more disadvantage situation may

occur at lower wind speed but longer acting time.

<table><caption>**Table 1.** Most disadvantageous Snow distribution factors</caption><thead><tr><th>Zone Number</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><th>μ<sub>r</sub></th><td>0.95</td><td>0.94</td><td>0.94</td><td>1.01</td><td>1.05</td><td>1.07</td><td>0.93</td><td>0.77</td></tr><tr><th>Zone Number</th><td>9</td><td>10</td><td>11</td><td>12</td><td>13</td><td>14</td><td>15</td><td>16</td></tr><tr><th>μ<sub>r</sub></th><td>0.99</td><td>1.01</td><td>1.09</td><td>1.19</td><td>1.23</td><td>1.22</td><td>1.10</td><td>0.95</td></tr></tbody></table>

# 4 Numerical simulation of snow-sweeping

In northeast China, the long winter causes durable snow load. What's more, the snow melts due to indoor heating and then frozen if the temperature goes down, this melt-frozen process causes the snow density increases with time which increases the roof load. For this reason, it is necessary to sweep the snow away especially for long-span roofs.

## 4.1. Geometrical modelling and parameters

In order to sweep away snow on the membrane roof, we designed an air blowing pipeline system see Figure 8. The pipelines are installed along with the boundary of the PTFE and ETFE membrane zones, each pipeline has several air-outlets. The pipeline has a circular cross section with a diameter of 10cm; the air-outlet has a rectangle cross section with a dimension of 12cm×6cm.

<figure><img src="image_62.png" /><figcaption>**Figure 8.** Sketch of air blowing pipeline system (a) layout plan, (b) three-dimensional effect graph</figcaption></figure>

The initial snow depth is set as 40cm which is about 1.23h₀, the air blowing velocity is selected as 2m/s, 3m/s, 4m/s, 5m/s respectively to estimate the effect of air blowing velocity on snow sweeping results.

The flow field changes rapidly and usually the snow
sweeping process finished in a short time, this process cannot be simulated by steady solver. So the transient solver is adopted which can obtain the residual snow mass concentration. The multiphase model is selected as Mixture model either, and time step in transient simulation is select as small as 10⁻³s for convergence requirement.

## 4.2. Results

The results are given as residual snow mass fraction and dimensionless snow mass flux defined blew.

The residual snow mass fraction is defined as (14), where *m* refers to the residual snow mass at time *T*, and *m*₀ refers to the snow mass at the initial moment.

$$f = \frac{m}{m_0} \qquad (14)$$

In order to make a comparison of snow sweeping efficiency in different cases, the dimensionless snow mass flux is defined as (15), where *Q* is the summation flux on all the outlets of the computational region, and the *Q*<sub>max</sub> is the maximum snow mass flux in all these 4 cases which occurs in the initial moment of the 5m/s case, *Q*<sub>max</sub>=16870 kg·s⁻¹.

$$q = \frac{Q}{Q_{\max}} \qquad (15)$$

<figure><img src="image_63.png" /><figcaption>Figure 9. Residual snow mass fraction</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>Figure 10. Snow volume fraction contour map of the 5m/s case</figcaption></figure>

Figure 9 compared the residual snow mass fraction varying pattern under four air blowing velocity cases, it is clear that the higher air blowing velocity is, the faster snow sweeping process is, thus the snow sweeping efficiency is positively related to the air blowing velocity which verifies our ordinary experience. Through further analysis we can see that the snow sweeping efficiency decreases rapidly with time, about 70% of snow is blew away during the first 4s.

If the residual snow mass fraction equals to 5% can be set as the competition criterion of the snow sweeping process, the completion time of 2m/s, 3m/s, 4m/s, 5m/s case is 14s, 12s, 10s, 8.5s respectively, which is approximately linearly related to the air blowing velocity. Figure 10 gives the snow volume fraction contour map at some key point in time of the 5m/s case.
<figure><img src="image_65.png" /><figcaption>Figure 11. Dimensionless snow mass flux</figcaption></figure>

Figure 11 gives the dimensionless snow mass flux. During $t=0\sim2$s, the dimensionless snow mass flux is positively related to the air blowing velocity too. But the varying tendency show some difference from the residual snow mass fraction after 4s, the curve show a negative correlation between the dimensionless snow mass flux and the air blowing velocity, this is because most the snow is blew away during the initial stage in the high air blowing velocity case, so the snow mass flux decreases rapidly and the higher air blowing velocity have a lower snow mass flux, but the difference is so small that doesn't affect the final completion order.

## 5 Conclusions

Based on the Mixture multiphase model, a three dimensional numerical simulation of snowdrift and snow sweeping were presented, the simulation obtained the most disadvantageous snow distribution factors of the membrane roof which indicates concave area (zone 10~15, see Figure 7) accumulates the largest snow loads, with the maximum snow distribution factor come up to 1.23. The snow sweeping simulation indicates 70% of snow were swept away during the first few seconds and the total time of the sweeping process is within 14s. In our practice operation, we can use a high velocity in the initial stage which has a high efficiency, after the majority of the snow were swept away, we can use a relatively low velocity to save energy.

## References

1. Tominaga Y, Okaze T, Mochida A. CFD modeling of snowdrift around a building: An overview of models and evaluation of a new approach. *BUILD ENVIRON*, 46(4), 899-910 (2011).

2. Uematsu T, Nakata T, Takeuchi K. Three dimensional numerical simulation of snowdrift. *COLD REG SCI TECHNOL*, 1(20), 65-73 (1991).

3. Inc. A. ANSYS FLUENT Theory Guide. (Ed.2011)

4. Beyers M, Waechter B. Modeling transient snowdrift development around complex three-dimensional

structures. *J WIND ENG IND AEROD*, 96(10-11), 1603-1615 (2008).

5. Beyers JHM, Sundsbø PA, Harms TM. Numerical simulation of three-dimensional, transient snow drifting around a cube. *J WIND ENG IND AEROD*, 92(9), 725-747 (2004).

6. Li X.F, Zhou X.Y, Gu M. Modified RSM turbulence model in snowdrift simulation. *14th National Wind Engineering Symposium*. Beijing, China. (2009)

7. Hong C.B. Numerical simulation of snowdrifting on typical form of long-span roofs. MSc thesis. Harbin Institute of Technology, (2012)

8. Pomeroy JW. Saltation of Snow. *Water resource research*, (1990).

9. Pomeroy JW. Steady-state suspension of snow. *J HYDROL*, (1992).

10. Chinese load code for the design of building structures (GB 50009). (Ed.2012)
# Developing experimental method for investigating snow deposition around buildings using snow substitutes

Jennifer Fiebig¹,ª, Hans Holger Hundborg Koss¹

¹Technical University of Denmark (DTU), Dep. of Civil Engineering, Building 118, 2800 Kgs. Lyngby, Denmark

**Abstract.** The phenomenon of wind-driven snow in the arctic and other snow-prone regions has been a significant problem for the built environment for centuries, but the integration of snow deposition in structural design and urban planning is still insufficient and challenging. A traditional way of studying snow deposition on and around buildings is wind-tunnel testing on a reduced scale using a substitute material for snow. This paper gives an overview of experimental approaches from the literature and presents the development of a test set-up in which the emphasis is on the physical properties and aerodynamic performance of the substitute material. The assessment of its fluid mechanical similarity is based on phenomenological comparison with Benchmark observations in the nature. The experimental study was carried out in the small boundary-layer wind tunnel at DTU Civil Engineering and full-scale observations are conducted on research sites in Greenland. The fall velocity of the substitution in comparison to the prototype snow could not determine a velocity scale for the experiment. On the other hand, the phenomenological effect of snowdrift is sufficiently accurate to replicate the deposition characteristic and the ground-area distribution.

## 1 Introduction

Adaption to extreme climatic conditions has been a challenge since the first human settlements in arctic and other snow prone regions. The evolution of traditional buildings and shelters was based on experience and led to optimised designs with respect to construction resources and environmental impact. However, the transition from settlements to cities with higher development density has led to modern building designs that deviate significantly from the traditional housing concepts and fail to take into account many aspects of environmental impact. One important aspect is the accumulation of wind-driven snow on and around buildings. The observed climatic trend of increasing temperature adds to the severity of this problem. Higher temperatures cause higher snow densities and an increasing tendency for combined snow and rain events resulting in wet and heavy snow. As a consequence, the loads from massive snow accumulation frequently exceed the carrying capacity of roof structures leading to local damages or to the total collapse of the supporting structure. Melt water penetrates the building causing deterioration of both the material structure and the indoor climate. Snow depositions around and between buildings block roads and building entrances and require frequent clearing.

Distribution and magnitude of these phenomena depend to a large extent on the air motion around buildings. Wind tunnel testing is a proven method for simulating the airflow interaction of building structures and provides a suitable platform for investigating snow transportation and deposition. However, experimental simulations of snowdrift in a boundary-layer wind tunnel struggle with scaling issues. A review of research on modelling snowdrift showed a parametric variability problem when to satisfy a large number of requirements [1]. A considerable amount of authors proposed accurate parameters for snowdrift in wind tunnel testing without

agreeing on the most important parameters to fulfil when simulating snow drift [2,3].

Additionally, the model of snowdrift differs from nature due to thermodynamic effects in the multiphase flow [4]. The mass is transported by the turbulent fluxes and because of the rate of phase change, ice to water vapour occurs. Additional simulation difficulties derive from the behaviour of snow consisting of different modes: a surface creep layer, a saltation layer, and suspension layer. [5,6] An observation in-situ shows the effect of the 2-phase flow like accumulation formation on facades, indicated in Figure 1 and the snow trace phenomenon in Figure 2. The simulation methods are still under development to fulfil the similitude laws within these layers and recent studies focus on particle density of the substitute in association with time duration of snow events and wind velocity hence large scale observation showed a correlation [2], [7–9].

<figure><img src="image_66.png" /><figcaption>**Figure 1** Snow deposition phenomenon, Hospital in Nuuk, Jennifer Fiebig 2015.</figcaption></figure>

General observation in snowdrift testing is the relevance of the accuracy in flow scale, geometric scale, particle saltation length, snowdrift rate and time scale to obtain representative results of the physical snowdrift

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

<footnote>ª Corresponding author: jenfi@byg.dtu.dk</footnote>
modelling [1, 2, 3, 4, 7, 8, 14, 15, 16]. The aim of the study is to establish a well-documented experimental setup suitable for an extensive parametric study of snow deposition and accumulation and sufficiently accurate to conclude not only on the ground-area distribution but also on the accumulation height. The setup will be used as main tool in a research project investigating the interaction between architectural designs, snow loading and ground-level accumulation between buildings. Furthermore, the documentation of the used material shall allow for benchmarking with numerical methods of snow drift simulation.

<figure><img src="image_68.png" /><figcaption>Figure 2 The tracing effect of snowdrift on a roof structure, Jennifer Fiebig 2015.</figcaption></figure>

## 2 Full-scale and experiment

### 2.1 Observation in situ

<figure><img src="image_69.png" /><figcaption>Figure 3 Snowflakes on millimetre paper with division 1 x 1mm Nuuk 2015.</figcaption></figure>

First in-situ observations of snow drift and accumulation were conducted on DTU premises in Sisimiut, Greenland, operated by the Center of Arctic Technology and Sustainable Solutions (ARTEK). Paralleling to the photometric documentation of snow accumulation on and

around buildings weather data was recorded using a medium level compact wireless weather station. The recoded data together with the in-situ observation provided information about the phenomenon in comparison to a model in the wind tunnel simulation. Based on the limited prototype data in-situ of snow deposition on a building over certain duration, an additional comparison from literature was needed. The density $ρ_p$ of prototype snow in the Antarctic has been measured by Mitsuhashi in 1982 with values from 500-900 kg/m³. The particle diameter for prototype snow in the Antarctic had a range of 100-150 µm [14]. On both sides in Greenland, Nuuk and Sisimiut, a millimetre paper method was selected in order to identify the size of snowflakes. Figure 3 illustrates the fallen snowflakes in Nuuk on a millimetre paper with 1x1mm division. Weather data was recorded on 15/12/2015 with an air temperature of -7°C, a humidity of 76 %, atmospheric pressure of 995.5 hPa and a wind direction from north. The classification of prototype snow in comparison to model snow validated the findings in the experimental development in respect to the similarity law for scaling.

### 2.2 Wind tunnel preparation

The study was carried out in the wind tunnel laboratory at the Technical University of Denmark (DTU), Department of Civil Engineering. The wind tunnel has a closed air circulation (Göttinger type) with a test section that can be operated in a fully closed, open and in blockage-tolerant configuration.

<figure><img src="image_70.png" /><figcaption>Figure 4. Overall layout of the small boundary-layer wind tunnel at DTU Civil Engineering.</figcaption></figure>

Figure 4 shows the overall layout and dimension of the wind tunnel. The test section has a cross-flow dimension of 0.48m in width and 0.48m in height and a length of 1.5m. The turbulent flow reflecting the main properties of the atmospheric boundary-layer (ABL) is generated in the flow processing unit upstream the test section. Generally the boundary-layer wind profile is defined by the logarithmic equation:

$$u(z) = \frac{u_s}{k} \ln \left( \frac{z}{z_0} \right) \quad (1)$$

Deviations from the logarithmic profile can occur if complex terrain or katabatic winds and see-breeze define the conditions [10], [11]. The roughness height $z_0$ for snow covered surfaces is varying with friction velocity and is to define as following:
$$z_0 = c \frac{u_*^2}{g} \qquad (2)$$

The value c is the constant for the different snow types with values of 0.015 found by Joffre (1982) and 0.018 by Schmidt (1982) for fresh cold snow[11]. Literature provides different measured roughness heights from a flat snow surface with $z_0 = 0.0001$ to a snow covered surface $z_0 = 0.002$ [10–12]. Additionally, the snow flux influences the wind profile shape and the wind velocity which was found by Föhn [13]. For the experimental set-up the mean wind velocity was measured without model snow (Figure 5). Eventually some airborne particles would fly into the sensitive openings of the measurement advice. Alternatively the simulation with model snow was measured with an Anemometer. The flow properties of the simulated atmospheric boundary-layer (ABL) have been measured analogically with a pitot tube and digitally with the cobra probe. In this experiment the exponent model (Eq.3) from the Eurocode for the normalization of the measured wind profile was applied:

$$\frac{u(z)}{u(z_{ref})} = \left(\frac{z-d_0}{z_{ref}}\right)^{\alpha} \qquad (3)$$

The factor α= 0.17 fitted a category between II and III for suburban and villages. This category matched with findings in literature for snowdrift simulation and also describes the rocky landscape and small villages in Greenland [14]. As earlier mentioned by Föhn the wind profile changes when snow particles fall into the flow. The effect was also considered in a snowdrift simulation by D.H. Kim, K.C.S. Kwok, and H.F. Rohde. The authors measured the generated wind profile without and with model snow and measured that the curve of the wind profile changed from category II to I and reflected thereby the Antarctic characteristic [14], [15]. In this experiment we assumed that the generated category between III and II did not have major changes and defined a snow covered surface in category II.

The wind tunnel simulation shall be capable of reflecting the ground-near transportation but also the accumulation on roof areas. Especially for the latter process the substitution material need to be airborne at higher levels to enter the flow field around roof tops and has to be able to accumulate on flat and on inclined surfaces. For this purpose a special seeding mechanism was installed to have control over the snow rate and simultaneously be able to test both, the suspension layer and the surface near saltation process.

At the beginning of the test section a set-up for the seeding technique with snow substitute was mounted on top of the wind tunnel (Figure 6). The seeding process is generated by a motor witch sets a wooden box with a mesh at the bottom, in motion. Through this technique the supplied particles of the substitute fall into the generated ABL flow and start drifting within the test section. For the snowdrift simulation a blockage-tolerant configuration was installed within the test section.

<figure><img src="image_71.png" /><figcaption>Figure 5 generated atmospheric boundary-layer (ABL)</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>Figure 6 Sketch of snowdrift configuration.</figcaption></figure>

## 2.3 Model and prototype snow Materials

Modelling snow with substitution material implies inevitably a simplification of the physical process regarding thermodynamics, wind-driven transportation and compacting. As a consequence, the ambition of the simulation does not aim on establishing a full physical similarity but rather a phenomenological similarity to observations in nature (chapter 2.1). Furthermore, it is not the author's intention to develop a simulation method from scratch but to draw from experience of earlier studies reported in literature. Here, all considerations take their origin in the understanding and/or description of the process in nature.

Snow is constantly in a sublimation process meaning that it melts and re-freezes again or even vaporises. Prototype snow is deforming into a round shape like grains undergoing a metamorphism process when it deposits on the surface [11]. This behaviour can be simulated by erosion/ saltation tests (e.g. Bagnold's sand erosion tests) using granular substitute materials suitable
for modelling of the phenomenon. The granular material is less than ideal to accumulate e.g. on a test steep roof model or other specific constructions. Further studies have been initiated to test the performance with powdery materials [14] in order to improve the representation of snow drift and snow deposition.

Following granular-type (saltation layer) materials have been specified so far: 1. Two types of Semolina (hard grains of flour), and 2. wheat bran and saw wood ash. As for the powdery substitute (suspension layer) the following materials have been looked at: baking powder, wheat flour and potato flour. The characterisation of the substitution material focused on following aspects:

- Geometric shape and size of the material particles

- Particle density

- Aerodynamic specification of fall velocity

Under a scanning electron microscope (SEM) different samples of grains and powder materials were investigated. Through this method, we were able to determine the average size of the particles and describe the structure and surface. With a SEM analysis (scanning electron microscope), zoomed images of the substitute samples are displayed by scanning it with a focused beam of electrons. Following samples have been scanned and analysed regard their size and shape:

Substitute as powders

- Baking powder

- Potato powder

- Wheat flour

The powders are fine particles and have been zoomed into 100, 500 and 1000 times. For the determination of the particle sizes the scans of 500 times zoom were chosen for a suitable estimation of the structure and size.

Substitute as grains

- Wheat bran

- Semolina

- Fine Semolina / Groats

- Saw dust

Grains are coarse and larger than the powdery materials, therefore the scans were zoomed in with 100 and 500 times. A photometric analysis determined the length and the radios of the particles. For the powders in 500 zoom, the division of the reference lineal was 100 µm, which is equal to 400px. The grains in 100 times zoomed had a reference lineal of 500 µm and equal to 400px. The pixel displays the optical reference of the different sizes in the picture. Figure 7 and 8 illustrate the general method of the measurement, both the grains and the powder materials.

<figure><img src="image_73.png" /><figcaption>Figure 7 Semolina under the microscope. Division of the reference lineal is 500 µm.</figcaption></figure>

After specifying structure and size of the particles a pycnometer method was accomplished to determine also the density of the substitute in g/m³:

$$ \rho = \frac{m}{v} \qquad (4) $$

For the determination the DS standard DS/CEN ISO/TS 17892-3:2004 was applied. For each substitute material 4g was supplied to the pycnometer. Eventually powder material would dissolve in water therefore ethanol was used to calculate the density of the substitutes.

<figure><img src="image_74.png" /><figcaption>Figure 8 Potato powder under the microscope. Division of the reference lineal is 100 µm.</figcaption></figure>

An aerodynamic specification of the fall velocity was determined as a third classification for the substitutes. The particle fall velocity was measured photometric and visualised with a laser beam. Following configurations with the Camera (Canon EOS 700D) have been adjusted: ISO 1600; Exposure 1/20 sec (0.05sec); f/9. In numerus pictures the mean fall velocity has been read from a reference lineal in the layer of the laser beam. The fall
velocity shall give an estimation of the velocity scale for the simulation:

$$L_v = \frac{v_m}{v_p} \qquad (5)$$

Particles have been supplied through the seeding mechanism and have been photographed while falling through a laser beam. Particles which fall into the laser beam layer got visible. The reference lineal indicates how long the particle falls within the exposure time of the camera. Figure 9 shows the length of falling potato powder particles in an exposure time of 0.05sec.

<figure><img src="image_75.png" /><figcaption>Figure 9 Determination of fall velocity m/s of Potato powder in 0.05sec.</figcaption></figure>

## 2.4 Test performance with snowdrift

After the specification of the substitution material and the preparation of the wind tunnel a test performance with substitution was simulated. Earlier researcher (Kind [6] and D.H. Kim et al. [15]) found the best agreement of similarity with high density materials. Sodium Bicarbonate with 1000 kg/m³ and a size of 50 µm was found to reproduce the most similar snow accumulation shape compared to the nature, while nine other different materials have been tested [15].

Based on the findings in literature that the particle size has to be considered to not be smaller than 30 µm, rather bigger we selected potato powder for the simulation [15]. The test model was a black concrete rectangular of 7 x 7 cm for a better contrast between the

white substitution and the model. The model was placed in the middle of the test section. A tachometer controlled the rotation speed of the seeding mechanism which was filled with the substitute. The supplied material of potato powder was 2,28kg in the seeding mechanism. With velocities from 1,3m/s to 3m/s snow accumulation in the middle of the wind tunnel section occurred. In order to control at which wind-speed accumulation in the middle of the section starts to build up and at the same time being able to measure while particles drift in the section, an anemometer was used. After the duration of 8:20 minutes a mass of 2,28kg fell rough the seeding process.

## 3 Results

As mentioned previously, the aim of the simulation is to develop an experimental method which can reproduce the phenomenological similarity of snow effects observed in the nature. In this section, we evaluated the specification methods before the test performance with substitution was carried out and compare it with data from literature.

The photographic measurement with millimetre paper in Nuuk shows an average diameter of snowflakes from 1000-3000 µm (Figure 3). Similar results were found in Sisimiut from 500-3000µm (data not shown). This result is significantly different from the measurements found in literature in the Antarctic, shown earlier in section 2.1. Because of the phase change, the snowflakes in-situ formed into a bigger formation. Nevertheless, the results do not show a correlation to the actual wind-driven snow events in Greenland and deposition of snow on and around structures. But they are parameters for similarity laws and define the scale in simulating the nature phenomenon. Therefore the investigation of different substitutes in size, density and the terminal fall velocity are specified. The results are summarised in Table 1.

<table><caption>Table 1 Specification of substitute size *r* or *l*, density *ρ*<sub>m</sub>, fall velocity *v*<sub>m</sub></caption><thead><tr><th>Substitute</th><th>d [µm]</th><th>ρ [g/cm3]</th><th>v [m/s]</th></tr></thead><tbody><tr><td>Baking powder</td><td>2.25-37.25</td><td>2.0717</td><td>0.98</td></tr><tr><td>Flour (Wheat)</td><td>15.0-75.0</td><td>1.4564</td><td>0.72</td></tr><tr><td>Fine semolina (mannagryn)</td><td>71.0-500</td><td>1.4199</td><td>1.28</td></tr><tr><td>Potato powder</td><td>10.0-62,5</td><td>1.4934</td><td>0.76</td></tr><tr><td>Saw dust</td><td>88-500</td><td>1.3517</td><td>0.33</td></tr><tr><td>Semolina</td><td>675</td><td>1.3128</td><td>0.95</td></tr><tr><td>Wheat bran</td><td>152.5-562.5</td><td>1.3176</td><td>0.75</td></tr></tbody></table>

<table><caption>Table 2 Overview of specification of snow</caption><thead><tr><th>Prototype<br>snow</th><th>d [µm]</th><th>ρ [g/cm3]</th><th>v [m/s]</th></tr></thead><tbody><tr><td></td><td>500-3000</td><td>0.500-0.900<sup>[14]</sup></td><td>0.3-2.0<sup>[15]</sup></td></tr></tbody></table>
As a comparison, the determination of the terminal fall velocity of snow with 0.3-2.0m/s depending on the snow type measured by Mellor 1965[15] is shown in Table 2. The results in Table 1 show that the sizes of the substitution are smaller compared to the prototype, but the density is higher and the terminal fall velocity is similar. That means a length $L_v$ for the velocity scale can't be defined by the aerodynamic specification of fall velocity and the particle size for model and prototype. In fact the significance and relevance of similarity has been discussed by other authors while reproducing the snow accumulation characteristic without fulfilling the physical requirements correctly, but experiments have shown similar snow formations which have been observed in the nature[1, 2, 3, 4, 7, 8, 14, 15, 16]

<figure><img src="image_76.png" /><figcaption>**Figure 10** Snowdrift with potato powder every 60sec with 3 m/s.</figcaption></figure>

From this starting point the snowdrift simulation was carried out without fulfilling the similarity law between model and prototype particle. Mainly the snow effects in suspension and saltation shall be reflected on a simple model shape like a cube. In agreement with the literature the wind profile should be generated correctly to transport the particles in their aerodynamic physical nature. The accumulation and deposition on the model in the test performance over the duration of 8.20 minutes is presented in Figure 11. First particles start to accumulate on the sides walls of the rectangular and shape pattern in a circle upwards. The edges of the rectangular show an aerodynamic flow pattern as in erosion tests (e.g. Bagnold's sand erosion tests). Deposition on the front

side accumulates and the side walls show a snow-flow pattern on the rectangular. The particles attach to each other and start to build a formation on the walls.

<figure><img src="image_77.png" /><figcaption>**Figure 11** Snow drift pattern after 8.20 minutes on a rectangular model.</figcaption></figure>

In general, aerodynamic snow-effects illustrated in Figure 11 are reproduced with the experimental method. Picture 6 in Figure 10 in the time lapse shows the trace of saltation process behind the rooftop. Likewise, the phenomenon of deposition of particles in front of the model is displayed. After the duration of 8:20 minutes the accumulation shape in the front is 4cm high and the side shape 2cm. On the backside of the rectangular the aerodynamic flow pattern is clearly printed.

## 3 Discussion

The integration of snow deposition in structural design and urban planning through an experimental method has been investigated mainly in the theoretical improvement of the similarity requirements. Furthermore, it is still insufficient and challenging to fulfil the physical laws of a 2-phase with snowdrift.

Literature and the research so far could not agree on a methodology to establish a correct physical similarity, as reported in the previous sections. So far the aspect of the phenomenological improvement have been largely ignored or wasn't in focus. The development of the documented experimental setup is suitable for an extensive parametric study of snow deposition and accumulation and sufficiently accurate to conclude not only on the ground-area suspension but also on roof structures in the suspension layer. The set-up will be the main tool for parametric simulations to investigate the interaction between architectural designs, snow loading and ground-level accumulation between buildings. The length $L_v$ was not found therefore the simulation was a dimensionless study with aspect on the characteristic of snow deposition on a structure. Here the significance of the relevance of snowdrift similitude parameters occurs, as other researcher discussed before [16]. The development of the experimental set-up allows a flexibility of reproducing not only the saltation process and the accumulation on the near surface ground.
Furthermore, the documentation of the used material and the parametric possibility allows benchmark, with numerical methods of snow drift simulation, to collect a wide range of data.

Architectural case studies with a parametric investigation in snow deposition, accumulation or even the effect of a trace line behind a roof structure can be simulated for data collection and optimizing the structural design of buildings. A similar tracing effect in literature [16], here called the reverse riming formation has been found and simulated from a reference case. Unlike this experiment, the snow particles started to shape a formation based on the long duration. The research in this paper and the comparison to literature shows that not all similarity laws are needed when simulating the phenomenon of snow accumulation and deposition on structural design.

## 4 Conclusions

Main significance in the future is the phenomenological accuracy of different shape formation which has to be studied in a parametric study with regular verification to the nature. An accurate generated atmospheric boundary-layer for simulating snowdrift is needed for modelling the nature phenomenon of wind-driven snow. Main focus is the phenomenological similarity of the snow deposition / formation without fulfilling all similarity laws for the substitution in regard to the prototype. This will forward the experimental development and the design process of urban planning in cold regions.

## 5 Acknowledgements

Financial support for the research project is partly provided by Laust Løgstrup, Director of Qeqqata Kommunia in Greenland, and fonds from the Technical University of Denmark DTU.

## References

1. F. Naaim-Bouvet and M. Naaim, Ann. Glaciol., **26**, 212–216 (1998)

2. W. Wang, H. Liao, M. Li, and H. Huang, Open J. Civ. Eng., **3**, 13–17 (2013)

3. F. Naaim-Bouvet, Surv. Geophys, **16**, 711–727 (1995)

4. R. A. Schmidt, Rev. Geophys. Sp. Phys., 20, 39–44 (1982)

5. R. J. Kind, J. Wind Eng. Ind. Aerodyn., 36, 855–866, (1990)

6. R. J. Kind, Cold Reg. Sci. Technol., 12, 217–228 (1986)

7. C.-E. Teleman, Constr. Arhit., **60**, 147–163 (2014)

8. X. Zhou, J. Hu, and M. Gu, Nat. Hazards, **74**, 1629–1648 (2014)

9. T. K. Thiis, J. Wind Eng. Ind. Aerodyn., **91**, 829–839 (2003)

10. S. Thordarson, PhD thesis, Norwegian University of Science and Technology NTNU, (2002)

11. C. Jaedicke, PhD thesis, University of Bergen, (2001)

12. J. W. Pomeroy and D. H. Male, Seasonal Snowcovers: Physics, Chemistry, Hydrology, Springer Netherlands (1987)

13. P. M. B. Föhn, J. Glaciol., **26**, 469–480 (1980)

14. Smedley, D. J., Kwok, K. C S, Kim, D. H., J.Wind Eng. Ind. Aerodyn., **50**, 153-162, 1993

15. D.H. Kim, K.C.S. Kwok, H.F. Rohde, Report, Similitude Requirements of Snowdrift Modelling for Antarctic Environment, University of Sydney, (1991)

16. K.C.S. Kwok, D.H. Kim, D. J. Smedley, H.F. Rohde, J.Wind Eng. Ind. Aerodyn., **41-44**, 2797-2808, (1992)
# Use of numerical simulations of snow drift in planning of infrastructure – A case study from Northern Norway

Thomas K. Thiis<fnref n="1" />, <fnref n="2" />

**Abstract.** A large industrial development is planned in Honningsvåg in Northern Norway. The area is characterized by subarctic climate with long cold winters. The vegetation on the site is sparse, and the wind exposure gives problems with snow drifting. This paper is a study of the snow drifting at the site and around some of the planned buildings. Based on numerical simulations and meteorological observations, snow depth maps of the site are produced for four different meteorological conditions.

## 1 Introduction

Generally, an assessment of snow drifts around buildings and in terrain should be performed in “scale levels”. The first important level is to determine which zones in the terrain that are naturally subjected to snow accumulation and large snow depth. When the most critical zones are identified and avoided, the second level would be an assessment of the position and orientation of the buildings on the site. The last level would be to identify any possible small scale problems including position of doors, sunshade and ventilation openings. This report mainly deals with the first level, snow distribution in the terrain.

## 2 Field site

The project in Honningsvåg involves a large development of the surrounding terrain as well as introduction of several new buildings of different sizes. The planning of the Veidnes project has not yet reached the detail phase, and the study is therefore concentrated on how the terrain influences the snow accumulation. In addition, the effect of some of the larger buildings has been included. Figure 1 shows a model of the planned project together with a wind rose showing wind statistics during the winter months. Figure 2 shows a panorama of the site taken the 2. March 2014 from the hill in the south-east part of the area. This is before the snow melting season has started.

The snow depth maps were made for four different meteorological conditions. These are the maximum short

2013/2014 in which there was a field inspection, the year representing the mean yearly cold precipitation and the extreme year, represented by the year with largest sum of winter precipitation since 1957 at the closest meteorological station. The wind conditions used in the simulations are the same for all the meteorological years and is explained later. Table 1 shows the meteorological conditions from the three selected years and the short-term maximum. Table 2 shows the distribution of joint precipitation and wind in the area since 1957. Most of the precipitation incidents occur with a precipitation intensity of less than 5 mm/24h and a wind speed between 5-10 m/s. The climate at the site is mild enough to prevent ice-covering of the fiord. This limits the fetch at the site to the open water body.

<figure><img src="image_78.png" /><figcaption>**Figure 1.** 3D model of study site. Terrain cutting in red.</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>**Figure 2.** Panorama view of site, seen from south-east.</figcaption></figure>

term continuous snowfall, which will be the design criterion for the snow clearing utility service, the year

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

<footnote n="1">¹Norwegian University of Life Sciences, Ås, Norway</footnote>

<footnote n="2">²Multiconsult ASA, Tromsø, Norway</footnote>

<footnote>ª Corresponding author: thomas.thiis@nmbu.no</footnote>
<table><caption>**Table 1.** Cold precipitation at Repvåg, close to Veidnes.</caption><thead><tr><th>Period</th><th>mm cold precipitation (-1 C)</th><th>Hour of precipitation at 5mm/12h</th></tr></thead><tbody><tr><td>Max. short term snowfall</td><td>36</td><td>86.4</td></tr><tr><td>2013/14</td><td>121</td><td>290.4</td></tr><tr><td>Mean</td><td>176</td><td>422.4</td></tr><tr><td>Extreme</td><td>382</td><td>916.8</td></tr></tbody></table>

Using Tabler's [1] guideline for determining yearly snow transport, the yearly mean snow transport is calculated to 24 tons/meter. The same guidelines classifies the snow transport to be controlled by the snowfall. Tabler describes 24 tons/meter as "light-to-moderate" snow drifting. To simplify the simulations, only the meteorological situation with cold precipitation are considered. In addition, the wind conditions during precipitation are divided in three speed classes and six directional classes, given in table 3.

<table><caption>**Table 2.** Joint distribution of cold precipitation (Air temperature<-1°C) and average wind speed during precipitation time and 12 hours after, from Slettnes and Repvåg (1957-2012), cold precipitation 31% of the time</caption><caption>Joint distribution of cold precipitation and average wind speed [% of cold precipitation time]</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="6">Sum of cold precipitation last 24 hours [mm water equivalent]</th><th rowspan="2">All amounts</th></tr><tr><th>0-5</th><th>5-10</th><th>10-15</th><th>15-20</th><th>20-25</th><th>&gt;25</th></tr></thead><tbody><tr><td rowspan="4">Average 10 min wind speed during and 12 hours after precipitation [m/s]</td><td>0-5</td><td>23.10</td><td>2.26</td><td>0.64</td><td>0.35</td><td>0.13</td><td>0.19</td><td>26.67</td></tr><tr><td>5-10</td><td>52.61</td><td>8.22</td><td>2.45</td><td>0.83</td><td>0.22</td><td>0.38</td><td>64.72</td></tr><tr><td>10-15</td><td>3.98</td><td>2.36</td><td>0.99</td><td>0.13</td><td>0.10</td><td>0.00</td><td>7.55</td></tr><tr><td>&gt;15</td><td>0.25</td><td>0.48</td><td>0.22</td><td>0.03</td><td>0.06</td><td>0.00</td><td>1.05</td></tr><tr><td></td><td>All speeds</td><td>79.96</td><td>13.32</td><td>4.30</td><td>1.34</td><td>0.51</td><td>0.57</td><td>100.00</td></tr></tbody></table>

<table><caption>**Table 3.** Wind climate at Veidnes during precipitation at Repvåg.</caption><thead><tr><th>Direction[deg]</th><th>0m/s</th><th>5m/s</th><th>10m/s</th><th>SUM</th></tr></thead><tbody><tr><td>15</td><td>2.2%</td><td>16.9%</td><td>2.6%</td><td>21.6%</td></tr><tr><td>75</td><td>2.6%</td><td>7.4%</td><td>0.8%</td><td>10.7%</td></tr><tr><td>135</td><td>1.4%</td><td>1.7%</td><td>0.2%</td><td>3.2%</td></tr><tr><td>195</td><td>4.6%</td><td>5.5%</td><td>2.1%</td><td>12.2%</td></tr><tr><td>255</td><td>4.4%</td><td>14.6%</td><td>5.3%</td><td>24.3%</td></tr><tr><td>315</td><td>2.0%</td><td>19.0%</td><td>7.1%</td><td>28.0%</td></tr><tr><td>SUM</td><td>17.1%</td><td>65.0%</td><td>17.9%</td><td>100.0%</td></tr></tbody></table>

## 3 Method

To simulate the spatial wind field and snow accumulation in the terrain, the general purpose finite volume CFD code ANSYS CFX 14 was applied. Twelve different characteristic weather situations were simulated in addition to the snow accumulation during calm weather. The simulations results in accumulation flux maps were snow accumulation [kg/m²s] is simulated for each grid cell

near the surface. The thirteen characteristic situations, (twelve situations of joint wind and precipitation and one with precipitation in calm weather) were then weighted according to the measured climate at the site (table 3) and combined to form maps showing the snow depth distribution at the end of years with different weather.

The CFD code solves the incompressible, time averaged Navier-Stokes equations, using the $k-\epsilon$ RNG turbulence model to close the set of equations. It is widely used in simulations of wind around buildings and it is capable of producing realistic results. The two phase fluid flow is solved with the Euler-Euler approach. Both the air and the snow phase possess its own flow field and the fluids interact via the drag force calculated using the Schiller-Naumann drag model. The simulation boundary conditions are assumed to be blowing snow with a precipitation intensity of 5 [mm/12h] which is equivalent to 1.157E-04 [mm/s] or a flux of 1.157E-04 [kg/m²s]. With an estimated terminal fall velocity of 0.5 [m/s] (Mason, 1962), the snow mass concentration in the air becomes 2.3E-4 [kg/m³]. In the present analysis, the density of snow particles when suspended in the air is set to 150 kg/m³ and the density of the snowpack is set to 300 kg/m³. The radius of the snow particles is set to 0.15 mm, a parameter that is mainly important for the determination of the particle/air interaction. The vertical wind velocity distribution, at the inlet boundary, follows the logarithmic expression in Eq. (1).

$$u(z) = \frac{u_{*}}{K} \ln \left( \frac{z}{z_{0}} \right) \quad (1)$$

where $u(z)$ is the wind speed at a given height $z$, $z_0$ is the aerodynamically height of the roughness elements, set to 0.0001m, $K$ is the von Karman's constant (equal to 0.41), and $u_*$ is the friction velocity set according to Table 4. Two wind speeds are simulated for each wind direction; these are 5 m/s and 10 m/s.

<table><caption>**Table 4.** Inlet velocity coefficients of the two simulation cases, according to Eq. (1)</caption><thead><tr><th></th><th>u*</th><th>u(10m)</th></tr></thead><tbody><tr><td>U<sub>5</sub></td><td>0.18</td><td>5</td></tr><tr><td>U<sub>10</sub></td><td>0.37</td><td>10</td></tr></tbody></table>

### 3.1 Simulation of snow erosion and accumulation

To simulate the snow depth distribution, thirteen characteristic situations are considered. For the twelve characteristic situations where wind is a participating mechanism, the erosion ($q_{ero}$) and deposition ($q_{dep}$) fluxes are calculated according to Naaim et al. [2]:

$$q_{ero} = B \rho (u_{*}^{2} - u_{t}^{2}) \quad \text{if} \quad u_{*} \ge u_{t} \quad (2)$$

$$q_{dep} = C w_{f} \frac{u_{*}^{2} - u_{*}^{2}}{u_{*}^{2}} \quad \text{if} \quad u_{*} < u_{t} \quad (3)$$

where $B$ is a coefficient representing the intergranular bonding in the surface layer, $\rho$ is the density of the air and
$u_*$ is the air friction velocity. By definition, $u_{*t}$ is the threshold friction velocity of the snow below which aeolian snow transport ceases; it ranges from 0.07 to 0.25 m/s for fresh and dry snow and from 0.25 to 1 m/s for old, wind hardened snow, and is here set to 0.3 m/s. This value, which is marginally higher than “fresh and dry snow”, allows for slightly higher deposition fluxes in areas of $u_* < u_{*t}$ and is meant to compensate for the exclusion of episodes of wind without precipitation. C is the snow concentration in the air with units [kg/m³] and $w_f$ is the terminal fall velocity of a snow particle. In the present simulation the value of $B\rho$ is set to $7*10^{-4}$ and $w_f$ is set equal to 0.5 m/s. The resulting net flux is then determined by

$$q_{erodep} = q_{dep} - q_{ero} \quad (4)$$

In this form, the effect of impinging particles on the erosion of the surface is ignored. Also the effect of the change of geometric shape of the snow surface is ignored. Naaim et al. [2] produced realistic simulations of snowdrifts around snow fences with this model. However, they state that, in leeward accumulation zones, the model is very sensitive to the ratio $u_*/u_{*t}$, implying that the quality of the air flow simulation near the surface is important. The method has previously been used to calculate snow drift on roofs [3]

## 4 Results and discussion

Based on numerical simulations and meteorological observations, snow depth maps of the site are prepared for four different meteorological periods. An example of the characteristic snow flux deposition map for 5 m/s and wind direction 315° is shown in figure 3.

<figure><img src="image_81.png" /><figcaption>Figure 3. Snow accumulation flux at Veidnes, wind direction 315°, wind speed 5 m/s.</figcaption></figure>

Figure 4 shows snow depths in 25 points measured with a snow sonde in March 2014, before the snow cover has started to melt. The simulations representing the winter of 2013/14 shows a variation of snow depth from close to zero to deeper than 2.5 meters, figure 5. This is confirmed

by in-situ measurements from the same winter given in figure 4.

<figure><img src="image_82.png" /><figcaption>Figure 4. Snow depth measurements [cm] at some point of interest 2. March 2014.</figcaption></figure>

<figure><img src="image_83.png" /><figcaption>Figure 5. Simulated snow depth distribution at the site, winter 2013/14</figcaption></figure>

The snow depth distribution at the site in a year with mean cold precipitation is shown in figure 6. The differences from the 2013/14 winter is small in the areas exposed to wind. However, in the terrain cutting close to the center of the site will experience considerably more snow. This is because the terrain cutting creates a sheltered area where blowing snow particles will accumulate. This is even more pronounced in the year of extreme cold precipitation, shown in figure 7.
<figure><img src="image_84.png" /><figcaption>**Figure 6.** Simulated snow depth distribution at the site, year with mean precipitation.</figcaption></figure>

The extreme conditions indicate that large parts of the planned terrain cuts will be covered in snowdrifts deeper than 2.5 meters. In such conditions, also the buildings at the site, located in the western part of the site, will be largely affected by the snow. The snowdrifts surrounding the buildings will exceed 2.5 m height. However, due to the large size of the numerical grid cells compared to the building structure in this area, the results are indicative. More accurate could be obtained with grid refinement in this area. The simulated snow depth in the area situated leeward of the hill, south-east in the simulation domain, is relatively deep compared to the rest of the domain. The terrain slope in this area is steep, some places more than 45°, and snow sliding is to be expected. The simulations do not account for snow sliding, thus the simulated snow depth in this area is presumably too deep. However, the accumulation potential indicates that this area might be subjected to avalanche hazard.

The last simulation concerns the snow accumulation during a short-term extreme snowstorm. This will serve as an estimate of the on-site capacity for snow clearing during operation. The meteorological measurements shows that the most severe snowing event since 1957 was 36 mm water equivalents in 24 hours. The wind direction of extreme precipitation events is usually north-west, thus 315° is chosen as the wind direction for this case. The wind speed distribution of the case is based on the mean wind climate during precipitation given in table 3, ie. 17% of the time 0 m/s, 65% of the time with wind speed 5 m/s and 18% of the time with wind speed 10 m/s. Note that the snow density of the accumulated snow is set to 200 kg/m³ in this case, which will influence the calculated snow depths. Figure 8 shows the snow depth distribution after such snowstorm. The results show that also a short-term extreme event will create snowdrifts which are more than 2.5 meter deep in some areas. Figure 9 shows a close up of the snowdrifts around the buildings in the west part of the site. The simulations indicate the position and size of the small-scale snow drifts around the buildings.

<figure><img src="image_85.png" /><figcaption>**Figure 7.** Simulated snow depth distribution at the site, year with extreme precipitation.</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>**Figure 8.** Simulated snow depth distribution at the site, short-term snowstorm.</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>**Figure 9.** Simulated snow depth distribution around buildings in the western part of the site after short-term snowstorm.</figcaption></figure>
## 5 Conclusions

Based on numerical simulations and meteorological observations, snow depth maps of Veidnes are produced for four different meteorological conditions. The snow drifting severity at Veidnes is classified as “light to moderate” according to literature. The simulations demonstrate how numerical simulations can be used as a planning tool in larger industrial sites where snow drifting is expected. Even if the in-situ measurements does not allow for direct comparison due to the low number of observations, the numerical method reproduces the main characteristics of the observed snow cover of the winter of 2013/14. This gives credibility to the use of the results for planning of the site.

The simulations representing the winter of 2013/14 shows a variation of snow depth from close to zero to deeper than 2.5 meters. The extreme conditions indicate that large parts of the planned terrain cuts will be covered in snowdrifts deeper than 2.5 meters. In such conditions, buildings at Veidnes will be affected by the snow.

It is recommended to avoid placing buildings and infrastructure close to terrain cuts. It is further recommended to install sufficient snow fencing to reduce snow accumulation around buildings and in the terrain.

## Acknowledgement

The author wants to thank STATOIL ASA for permitting the use of the case-study in this paper. Thanks to Juliane Borge and Juni Vaardal-Lunde for the climate data analysis.

## References

1. Tabler R.D. (2003) Controlling blowing and drifting snow with snow fences and road design, NCHRP project 20-7 (147)

2. Naaim, M., Naaim-Bouvet, F., Martinez, H. (1998) Numerical simulation of drifting snow: erosion and deposition models, Annals of Glaciology 26: 191-196.

3. Thiis, T.K., Potac, J., Ramberg, J. F. 3D numerical simulations and full scale measurements of snow depositions on a curved roof.: *Proceedings from 5th European & African Conference on Wind Engineering*. : Firenze University Press 2009 ISBN 9788864530383. s. 371-374
Session 2

Structural loading I
# Probabilistic concepts in snow engineering - from observations to the specification of consistent design values including climate change

Michael Kasperski

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

**Abstract.** Consistent probabilistic concepts are required to specify the design snow loads for structures. The paper presents some basics on probability theory and discusses the fundamental reliability demands. For complex problems simulation techniques can be applied. The paper explains how a bad quality of the generated random numbers can ruin a simulation. For the snow load on the ground, extreme value statistics can be performed based on annual extremes or based on independent peaks over a threshold. The paper points out that both approaches lead to the same probability information. In Europe, exceptional snow loads have been considered as design scenario. The respective events have been identified as outliers in the ensembles. The paper reveals that the strategy used in the scope of the European snow research project fails and may lead to considerable underestimations of the design snow loads. A consistent strategy is presented how to deal with outliers. The snow loads on roofs are influenced by the wind, thus making the shape factor a random variable. A probabilistic model is available for gable roofs. The paper presents corresponding results for the design values. Furthermore, the combination of wind and snow is discussed. Finally, the question of possible long term trends in the snow climate is studied.

## 1 Introduction

Probability theory is a branch of mathematics which is required in structural engineering as soon as it comes to the specification of design values for action effects or the resistance. The long way from observations or experimental results to the final design values requires a lot of choices which are difficult to evaluate in their consequences since humans have no natural understanding of the concept of probability.

This unfortunate shortcoming in the evolution of human thinking can be illustrated by the following example which refers to a famous TV-show in the USA and is known as the Monty-Hall dilemma. After several challenges the final candidate finds him- or herself in front of three doors. Behind one of the doors, there is the fantastic first price; the room behind the two other doors is empty. The candidate is asked by the hosting show-master Monty Hall to choose a door. Then the show-master - knowing where the first price is hiding - opens one of the remaining doors which with no surprise gives view to an empty room. Finally, the show-master asks the candidate if he or she remains with the initial choice or if the candidate wants to change the door. The puzzling question arises which is the best strategy to optimize or maximize the probability of winning the first price: Staying faithfully with the initial choice or fickly make a new choice? Most people have difficulties to find the answer for this problem of conditional probabilities; and it is easy to convince a layman by simply claiming that it doesn't matter which door is chosen since when two doors are still closed, the probability is 50:50 for getting the first price anyway.

The second section of this paper therefore briefly summarizes some basics in probability theory and presents the fundamental reliability demands which are applied in structural design.

When mathematical skills are missing or when problems become too complex for an analytical approach, researchers and engineers today can apply the modern tool of simulation. The core of simulation techniques is an efficient and reliable pseudorandom number generator. Basic demands will be discussed in section 3.

Section 4 presents the probabilistic models for the snow loads on the ground. Two approaches are discussed: the extreme value statistics based on annual extremes and the so-called peaks-over-threshold approach. The latter introduces as further variable the number of events per year. If this variable is modelled by the Poisson distribution both approaches are identical.

In the European snow load codification, exceptional snow loads have been included. These loads have been identified as statistical outliers. Section 5 discusses the statistical stability of identifying outliers in a confined ensemble and shows that with the wrong strategy estimated design snow loads become too small.

Snow and wind belong together like wind and waves. The accumulation and distribution of snow is influenced by the wind. Furthermore, in the design, combinations of these two climatological actions have to be considered. Section 6 presents the influence of wind on the design values of snow on gable roofs and discusses the derivation of combination factors. A consistent approach is presented which specifies a single combined wind-snow load case.

Over the recent decades, winters in Germany have been different to what they used to be in the past. Missing snow falls torture the tourist industry, and in many tourist areas snow canons are used to obtain at least small stripes of snow to allow for skiing. The question arises, if it is possible to identify long term trends in the snow climate which justify reducing the design snow loads. Section 7 presents a probabilistic approach to this problem.

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

<footnote>a Corresponding author: michael.kasperski@rub.de</footnote>
## 2 Basics in probability theory and fundamental reliability demands

The classical definition of probability is based on a chance experiment which can result in N mutually exclusive and equally likely outcomes. If $N_A$ of these outcomes result in the occurrence of the event A, then the probability of A is simply obtained as the ratio of the number of successful outcomes to the number of all possible outcomes:

$$p(A) = \frac{N_A}{N} \qquad (1)$$

A simple example is given with the chance experiment of throwing a die. The set of all possible outcomes, called the sample space, is obtained as the integer numbers from 1 to 6. If the event A is defined, for instance, as getting a result larger than 4, two of the six outcomes are successful. So, the probability of A becomes $2/6 = 1/3$.

The generalized definition of probability which can be applied to any random process is obtained as the relative frequency of the event A which is the ratio of $n_A$ occurrences in N observations for N going to infinity:

$$\lim_{N \to \infty} p(A) = \frac{n_A}{N} \qquad (2)$$

$n_A$ - occurrences of event A

For a discrete process, the event A can be defined as a set of specific outcomes; for a continuous process, the event A has to be specified as a set of specific ranges of outcomes. Basically, the probability lies in the range from 0 to 1. The probability 0 is obtained for the impossible event; the probability 1 corresponds to the certain event.

The complement of an event 'A' is the event 'not A'. The event A and its complement are mutually exclusive. The probability of the complementary event is obtained as:

$$p(\text{not } A) = 1 - p(A) \qquad (3)$$

The distribution of the probability for a random process x over the range of all possible outcomes is called the probability density function f(x). For a continuous process, the probability is obtained as the area under the probability density in the ranges of the specified outcomes. The integral over the probability function from $-\infty$ to x is called the cumulative probability distribution or often simply the probability distribution F(x):

$$F(x) = \int_{-\infty}^{x} f(u) du \qquad (4)$$

The basic relation between the probability density and the cumulative probability distribution is shown in figure 1. The area under the probability density function f(x) from $-\infty$ to a specific value x occurs as the amplitude of

the cumulative probability distribution F(x). The cumulative probability distribution F(x) therefore describes the non-exceedance probability of a specific value x. Basically, the probability density has no negative values, and F(x) is monotonically increasing from 0 to 1.

If a chance experiment, having the probability $p_1$ for a specific event A, is repeated several times, the probability that the event A is occurring increases with each further trial. From everyday experience, it is obvious that the corresponding accumulation of probability has to be non-linear. The number of trials until event A finally occurs is smaller than infinity; however, it may by a large number.

For L repetitions the accumulated probability $p_L$ for the event A occurring at least once in the L trials is obtained as:

$$p_L(A) = 1 - (1-p_1(A))^L \qquad (5)$$

Figure 2 shows the basic accumulation behaviour on the example of a perfect coin, a perfect die and the once-in-fifty-years event. The respective events are: occurrence of head, occurrence of six and exceedance of the once-in-fifty-years event.

<figure><img src="image_89.png" /><figcaption>probability density function f(x)</figcaption></figure>

<figure><img src="image_90.png" /><figcaption>cumulative probability distribution F(x)</figcaption></figure>

**Figure 1:** Relation between the probability density function f(x) and the cumulative probability distribution F(x)
<figure><img src="image_91.png" /><figcaption>**Figure 2:** Accumulation of the probability of event A occurring at least once in L trials</figcaption></figure>

Equation (5) specifies the average probability of event A occurring at least once in L trials. The number of occurrences of event A in L trials forms a random sequence. The corresponding probability density is given with the Bernoulli distribution:

$$f(N) = \binom{N}{L} \cdot p^N \cdot q^{L-N} \quad (6)$$

L - number of trials

N - number of successful trials

p - probability of success

q - probability of failure = 1 - p

For the event of getting no head for throwing the coin twice, the sample space is no head, one head and two heads. The corresponding probabilities are 0.25, 0.50 and 0.25. For the die-experiment the sample space and the corresponding probabilities are summarized in table 1 assuming a series of six throws.

<table><caption>**Table 1:** Probabilities of the random number of successful trials in 6 throws for the event getting a six for a perfect die</caption><tr><td>N<sub>6</sub></td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td></tr><tr><td>p</td><td>0.34</td><td>0.40</td><td>0.20</td><td>0.054</td><td>0.0080</td><td>0.00064</td><td>1/6<sup>6</sup></td></tr></table>

N<sub>6</sub> - number of successful trials in 6 throws

The number of trials k needed to get one success follows the geometric distribution:

$$f(k) = (1-p)^{k-1} \cdot p \quad (7)$$

For experiments with a small probability of success p, i.e. for rare events, and a great number of repetitions L, the Bernoulli distribution can be simplified and becomes in the limit for p → 0 and N → ∞ the Poisson distribution:

$$f(N) = \frac{\lambda^N}{N!} \cdot \exp(-\lambda) \quad (8)$$

λ - average number of occurrence during the given time interval

For the once-in-50-years event, N<sub>av</sub> becomes 1 for 50 years observation. The probability of getting no exceedance of the once-in-50-years event in 50 years observation is 0.368. This approximated value differs slightly from the true value of 0.98<sup>50</sup> = 0.364. One single exceedance in 50 years is obtained with the probability of 0.368, two exceedances occur with a probability of 0.184, three exceedances with 0.061, four with 0.015, five with 0.003 and so on.

The inter-arrival times of a Poisson process follow the exponential distribution which is given as follows:

$$f(T) = \frac{1}{\lambda} \cdot \exp\left(-\frac{T}{\lambda}\right) \quad (9)$$

$$F(T) = 1 - \exp\left(-\frac{T}{\lambda}\right)$$

with 1/λ - average arrival rate

The average arrival rate is 1/50 for the once-in-50-years event. The corresponding probability density and cumulative probability distribution are shown in figure 3.

<figure><img src="image_92.png" /><figcaption>probability density f(T)</figcaption></figure>

<figure><img src="image_93.png" /><figcaption>cumulative probability distribution F(T)</figcaption></figure>

<caption>**Figure 3:** Inter-arrival time of exceedance of the once-in-50-years event</caption>