Introduction

Reference Uematsu, Kaneda, Takeuchi, Nakata and YukumiUematsu and others (1989) tried to research drifting snow phenomena ranging from air flows to snowdrifts. Their models were two-dimensional except in the case of the suspension phenomenon. In this research the authors were not able to reproduce a snowdrift in a complicated topography, which points out the need for consideration of the suspension phenomenon. Recently there have been a few papers on this topic, but none of much significance. No numerical simulations have considered the air flows which lead to drifting snow as well as both the saltation and the suspended motion of snow particles.

Numerical Model

Basic equations

Air flow

Snowdrift development is mainly controlled by the friction velocity u_*, which can be obtained from a velocity profile.

A numerical simulation method developed by Reference ArisawaArisawa (1987) is applied to obtain a three-dimensional velocity field. Arisawa’s method is a so-called control volume method based on Patanker’s algorithm (Reference PatankerPatanker, 1981). The Reynolds equations with boundary conditions are numerically solved using this method.

Snow transport

Snow transport is possible by saltation, suspension or creep. As mentioned above, only saltation and suspension are taken into consideration in this simulation.

Up to now, many empirical formulae describing snowdrift transport rate by saltation have been issued (e.g. Reference KobayashiKobayashi, 1972). However, the application of these empirical formulae are limited strictly to the place and time of the observation because they contain empirical constants. Therefore, in view of future applications, a new theoretical equation should be derived from the particle trajectory.

Assuming that the saltation phenomenon occurs continuously, the following four equations for saltation of snow particles are derived. The parameters used in the derivation of the trajectory of saltating particles are shown in Figure 1.

The shear stress T ₀ at the snow surface is expressed by the following formula

where TS 0 is the shear stress caused by collision of saltating particles with the surface and T _t is the shear stress due to wind. Here, T _t is assumed to be the critical shear stress, which is consistent with the assumption of Reference OwenOwen (1964).

TS 0 can be expressed as

where M is the mass flux per unit area and time, U ₂ is the horizontal component of the particle velocity before collision with the surface, and U ₁ is the horizontal component of the particle velocity after collision. The condition for snow particles undergoing constant saltation can be given by the equation

Here, m is the mass of a single particle given by (β/6)πd ³, d is the diameter of the particle and β is the density of ice. Ν is the number of particles passing per unit width per unit time, and L_s is the saltation length. The lefthand side of Equation (1.3) gives the number of particles passing per unit area per unit time, and the righthand side gives the number of particles landing on the surface.

Fig.1 . Sketch for the definition of trajectory parameters.

Finally, the snowdrift transport rate qB is defined as

Using Equations (1.1), (1.2), (1.3) and (1.4), the following theoretical equation for particle transport rate can be derived:

Here, H _s is the saltation height, α is the fly-out angle, e is the restitution coefficient, S = /ρ — 1, g is the acceleration due to gravity, and u _*t is the threshold friction velocity. Although there are many variables in Equation (1.5), each has a physical meaning.

Equation (1.5) is checked in section 3 by comparison with observations. The values for L _s, H _s, d and e are given as follows: e = 0.8, obtained from the collision experiment on the ice particles against the ice surface (Reference Araoka and MaenoAraoka and Maeno, 1978); the other three are constant values since L _s is a function of a friction velocity, judging from the work of Reference KikuchiKikuchi (1981). As a result, Equation (1.5) is reduced to

Equation (1.6) is compared with the data collected by Reference KobayashiKobayashi (1972). Figure 2 shows that this semi-theoretical equation compares well with measured values.

Fig.2 . Snowdrift transport rate against friction velocity by Kobayashi (1972) and the simulated curve of the semi-theoretical formula.

Suspension is modeled using the diffusion equation.

The deposition (or erosion) rate can be calculated by simulating the drift density in the saltation layer. The deposition rate, d, can be expressed as

Here, W _f is the snowfall velocity. In Equation (1.7) the maximum erosion rate is assumed to be equal to the maximum deposition rate.

The maximum erosion rate, e _r, can be expressed as

Finally, we can calculate the snowdrift rate, s, as

where s is defined as the mass of snow accumulation on a unit horizontal area per unit time.

2. Numerical aspects

Various methods for numerical calculation in fluid mechanics are found in Reference RoacheRoache’s textbook (1978). Among them, the “control volume method” is presumed to be the most physically precise. The SIMPLER (Semi-Implicit Method for Pressure-Linked Equations Revised) model, which Patanker (1981) proposed, is used here.

3. Initial and boundary conditions

This simulation is roughly divided into two parts, flow and snowdrift density. Each initial and boundary condition is applied here. To determine the snowdrift density, the diffusion equation must be solved. In doing so, setting the boundary and initial conditions becomes a problem.

Air flow

The boundary conditions are now stated.

At the lower boundary

At the top of the domain

and

where u ₀ and v ₀ are the initial wind speeds at the top of the domain and (ρ)t is the initial air density.

At the lateral boundaries, the radiation condition of Reference OrlanskiOrlanski (1976) is adopted.

The logarithmic distribution is adopted as the initial condition over the entire domain.

Snow transport

The lower boundary of the suspension layer is in contact with the saltation layer; therefore, the mean snowdrift density of the saltation layer is considered to be the lower boundary condition.

Strictly speaking, the height of the saltation layer varies with the degree of development of the snowdrift; however, it is assumed to be approximately constant at all points.

The snowfall rate is used as the upper boundary condition for suspension. The upper boundary of suspension is the height where the snowdrift density becomes zero. The upper boundary condition is expressed as

where P _re is the snowfall rate.

The snowdrift density is used as the initial condition over the entire domain.

The other variables

The five elements shown in Table 1 are input values for the other variables in the numerical simulation. The wind speed and snowfall rate are explained as the initial condition (described above). The roughness height is used to give the initial condition for air-flow calculation and at the same time to acquire the friction velocity on the snow surface from the wind-speed distribution obtained. The threshold friction velocity is used to acquire the drifting-snow transport rate of the saltation layer. Furthermore, the falling velocity is used to convert snowfall rate and the drifting-snow transport rate into the snowdrift density.

4. Flowchart of simulation

Figure 3 outlines the flow of the calculation process. First, geographical features are taken as a boundary for the calculation of air flow. The calculation of air flow is repeated until the wind velocity given by the logarithmic distribution converges.

Fig.3 . Flowchart of simulation.

Next, the snowdrift density is calculated using the previously obtained air flow. The calculations are repeated until a converged value is obtained using the snowdrift density obtained from snowfall rate for the initial condition. As a result, the snowdrift density and the snowdrift rate are obtained. A snowdrift is simulated by integrating the snowdrift rate. A snowdrift at an arbitrary time can be obtained by repeating the above calculation for a necessary period of time.

Comparison of Simulation with Observation

The above method enables us to simulate a variety of snowdrift cases. In this section, the results of the simulation are compared with the observations of the past. The result of drifting snow simulation on a plain is described because it has the simplest topographical setting among those analyzed, and much research of drifting snow has been conducted for this case.

Figure 4 shows the mesh used for the calculation, the snowdrift density (g m⁻³), and the snowdrift rate (g m⁻²s⁻¹). The horizontal axis shows the distance from the origin of the snowdrift, and the vertical axis shows the height above the snow surface. Figure 4a is the case for a particle diameter of 0.3 mm, and Figure 4b for the particle diameter of 0.1 mm. The area of oblique lines in the drift-density graph shows where the snowdrift density is greater than the snowfall density. Figure 4a and 4b shows that the snowdrift density becomes larger leeward from the origin of the snowdrift and decreases gradually as the depth increases. The depth of the high-drift density layer is 0.7 m for a particle diameter of 0.3 mm, and 4.25 m for that of 0.1 mm. There is almost no change in the height and density of the high-drift density layer. The windward lateral boundary is 40 m from the origin of the snowdrift for the particle diameter of 0.3 mm, and 110 m for that of 0.1 mm. For both particle diameters, the leeward snowdrift rate increases gradually from zero to the snowfall rate towards leeward. The snowdrift rate reaches the snowfall rate at approximately the same point that the drift density becomes constant.

Fig.4 . Calculated drift density and snowdrift rate in a plain with wind velocity of 10 m s⁻¹ at height of 10 m above snow surface, (a) Snow particle diameter is 0.3 mm. (b) Snow particle diameter is 0.1 mm.

Figure 5 shows the drift-density profile, which was obtained at steady state. The solid circles show the simulation result and the other symbols show observations from the past. The drift density decreases nearly linearly up to a height of 50 cm, and becomes constant above 50 cm with the logarithmic vertical axis. This profile is compared with the observations of Kobayashi (1979), and Reference Kaneda and UematsuKaneda and Uematsu (1982) (Fig. 6). The simulated profile falls between the density profiles for wind speeds of 9.8 m s⁻¹ and 5.0 m s⁻¹ at the height of 1 m above the snow surface (Kobayashi, 1979). Since the simulation was done for a wind speed of 7.5 m s⁻¹, the result is considered adequate.

Fig.5 . Comparison between observed and simulated drift density.

Figure 6 shows the variation of snowdrift transport rate relating to the distance from the origin of the snowdrift. The snowdrift transport rate increases in the leeward direction and becomes constant for 30 m and beyond. This tendency is consistent with the observations of Kobayashi (1972), which were obtained at wind speeds of 4.5–13.3 m s⁻¹ at 1 m above snow surface with particle diameters of 0.1−0.5 mm. Thus, this newly developed model agrees well with past observations. Reference TakeuchiTakeuchi (1980a) observed that snow transport rate becomes constant 350–400 m from the origin of the snowdrift while Kobayashi (1972) observed that point to be 30–60 m from the origin. This contradiction may be explained by the difference in particle diameter, as shown above.

Fig.6 . Comparison between observed drift transport rate (Kobayashi (1972)) and simulated rate.

Discussion

The theoretical formula of the snowdrift transport rate by saltation is examined first. Equation (1.5) for the snowdrift transport rate by saltation contains the following coefficients and constants: saltation distance; saltation height; angle of departure; restitution coefficient; diameter of snow particle; air density; density of snow particle; friction velocity; and thereshold friction velocity at which the snow particles begin to move. It is important to know which coefficients or constants make a main contribution to the snowdrift transport rate.

First, the ratio of the saltation height to the saltation length ranges from 0.1 to 0.4 as shown in Figure 7; therefore, the snowdrift transport rate ranges from 1.6 to 3.1. Next, the angle of departure is roughly in the range from 20° to 80°, and the transport rate varies from 0.4 to 5.6. Moreover, the diameter of the snow particle is in the range of 0.1–0.5 mm, reported by Kobayashi (1972), and the snowdrift transport rate varies from 0.01 to 0.02. The coefficient of restitution for the snow particle and the ice surface is roughly in the range of 0.2–0.9 by Araoka and Maeno (1978), therefore, the snowdrift transport rate varies from 0.25 to 9.0. The threshold friction velocity calculated with the diameter of the particle by White’s formula (1940) is in the range of 0.2–0.4 m s⁻¹, and the snowdrift transport rate ranges from 0.04 to 0.16. The friction velocity is in the range of 0.2–1.0 m s⁻¹, and the snowdrift transport rate varies from 0.04 to 1.0. When the maximum value of the contribution is divided by the minimum value, the resultant value indicates the degree of change in the snowdrift transport rate. The resulting ratios are: 1.9 for the ratio of saltation height to saltation length; 14 for the angle of departure; 2 for the diameter of the particle; 4.5 for the restitution coefficient; 4 for the threshold friction velocity; and 25 for the friction velocity. These values for the friction velocity, the angle of departure, restitution coefficient, and the threshold friction velocity are larger than 4; however, the angle of departure as well as the restitution coefficient are distributed near the average values as shown in Figure 7 and do not influence the snowdrift transport rate. On the other hand, the range of the friction velocity is very common, as reported by many researchers. Its accuracy must influence the presumed accuracy of the snowdrift transport rate. Therefore, it is important to compute the friction velocity precisely in this simulation model.

Fig.7 . Saltation length, angle of departure and the ratio of saltation height to saltation length against wind speed, obtained by the wind tunnel experiment (Kikuchi, 1981).

Concluding Remarks

In this study, a numerical simulation model was applied to the drifting snow and snowdrift phenomena. The following results were obtained:

A new theoretical expression for the saltation movement of snow particles is derived from the ballistic study, and thereby, the amount of drifting snow transport in the saltation layer is calculated. The computational results show good agreement with an observed result.

A three-dimensional numerical simulation model is developed for drifting snow and snowdrifts. The simulated flow is determined as follows: (1) Air flow is obtained by solving Reynolds equation and the continuity equations. (2) Using the airflow obtained, the diffusion equation is solved by introducing a semi-theoretical expression of saltation for a lower boundary condition. (3) From the computed solution of the diffusion equation, the accumulation of snow is defined as those snow particles that are not transported by saltation.

The developed numerical model was applied to the drifting snow phenomena on level ground.

Various observed snowdrifts were reproduced using this model. In the case where the snowdrift itself changes physical conditions of the air flow, the model is run repeatedly, giving the newly formed drifting snow surface as the lower boundary.

Acknowledgements

I wish to express my thanks to Drs D. Kobayashi, N. Maeno, N. Ishikawa, R. Narue, K. Nishimura, Y. Kodama and Y. Ishii of the Institute of Low Temperature Science, Hokkaido University, and Dr Κ. Chikita, Department of Geophysics, Hokkaido University, for their helpful discussions and encouragement throughout this study.

I am also grateful to Mr Y. Arisawa, Mr K. Takeuchi, Mr T. Nakata and Mr T. Kaneda of Japan Weather Association for their helpful assistance. I am indebted to Ms Y. Yamada, Japan Weather Association, for her assistance in preparing the manuscript. Finally, I also sincerely thank Mr H. Denpo, Mr M. Yamazaki and Mr M. Satake, Japan Weather Association, for their encouragement.

This paper is the abstract of a Ph.D. thesis, Hokkaido University, 1992.