Introduction
Drifting snow is an important topic in both the engineering and climatological fields. Although recent research has revealed many basic properties of drifting snow, the time needed to reach steady-state conditions, including not only the saltation but also suspension, remains controversial. This is due to a lack of observational data, and the difficulty of expressing the transition between saltation and suspension using the majority of the current mathematical models.
In this study, the time evolution of drifting snow under a steady wind with turbulence was estimated using a numerical model. In the model, Lagrangian stochastic theory is used to incorporate the effect of turbulence on the motion of drifting-snow particles (Reference Sato, Uematsu and KanedaSato and others, 1997; Reference Sundsbø and HansenSundsbø and Hansen, 1997;Reference Taylor, Li and WilsonTaylor and others, 2002).This method enables us to discuss both the saltation and the suspension process. Aerodynamic entrainment, grain/bed collision (splash process), wind modification, and particle size distribution, according to the self-regulatory saltation model of sand (Reference Anderson and HaffAnderson and Haff, 1988, Reference Anderson and Haff1991;Reference McEwan and WillettsMcEwan and Willetts, 1991, Reference McEwan and Willetts1993), are also taken into account. The calculated time-scale results will be discussed in conjunction with the development of a spatial length scale of drifting snow, and will be compared with results reported in other studies on the spatial development of drifting snow.
Numerical Model
The model used in this study is one-dimensional and time-dependent. It calculates each particle trajectory under a steady wind with turbulent fluctuation over a flat, uniform snow surface. Although this approach requires large computational loads and cannot be applied to large areas (Reference GauerGauer, 2001), it is physically more sound than other methods and is useful for exploring the basic dynamics of drifting snow.
Furthermore, the model includes four distinct sub-processes (aerodynamic entrainment, the trajectory of the grains, splash process and the wind velocity modification) which are essential for aeolian sediment transport.
Particle trajectories are calculated within a small domain over a flat surface with periodic boundary conditions, under a horizontally uniform wind (Fig. 1). In the domain, snow particles start to move from random positions on the bed, owing to the aerodynamic entrainment. The number of entrained grains is calculated at each time-step.When a particle collides with the surface, the number, the velocities and the size of ejected grains are determined by splash functions. Thermodynamic effects (sublimation of the snow particles, humidity, temperature, etc.) are not taken into consideration; specific attention is paid to the grain dynamics in this study. The particles passing the end of the grid at a given height were counted in order to estimate the vertical distribution of the streamwise mass flux. The wind modification owing to the momentum exchange between the grains and the wind is also calculated at each vertical grid.
Fig. 1. Schematic diagram of the calculation configuration. A periodic boundary condition was applied for the trajectory calculation (dashed line).
Within the atmospheric surface layer over a plane, the mean wind velocity is typically large with respect to the horizontal wind fluctuation, and the mean vertical wind is zero. Therefore, in the two-dimensional case, the wind velocity vector u can be approximated by u = (U;w′), where U is the mean wind speed and w′ the vertical wind fluctuation. Themeanwind field follows the Reynolds-averaged Navier–Stokes equation for a steady state (@/@t = 0) and horizontal homogeneity (@/@x = 0). Using Prandtl’s mixing length theory, we obtain the horizontal momentum equation for air
where ρf is the density of air, and κ the von Kármán constant. The coupling force Fx corresponds to the rate of the time change of particle momentum as
where n is the total number of particles per unit volume at height z;mthe particle mass, and (dup(z)/dt) the horizontal acceleration of a particle at height z above the surface.
w′ is determined by the random-flight model (e.g.Reference Wilson and SawfordWilson and Sawford, 1996). For a fluid parcel i, w′i at time t þ Δt is expressed as
where Δt is a time increment, σw is the standard deviation of the vertical wind, ξ is a Gaussian random variable with zero mean and unit variance and TL is the Lagrangian time-scale. In the surface layer, TL can be expressed in terms of friction velocity u* and height above the surface z as TL = 0.4z/u* (Reference Wilson and SawfordWilson and Sawford, 1996). Lagrangian tracking techniques are suitable for turbulent particle dispersion, and a similar method has been applied to the simulation of blowing- and drifting-snow particles emitted from a fixed point, by Reference Sato, Uematsu and KanedaSato and others (1997), Reference Sundsbø and HansenSundsbø and Hansen (1997) and Reference Taylor, Li and WilsonTaylor and others (2002).
In the calculation of particle trajectories, we should obtain the vertical wind fluctuation w0pi at each particle’s position, which is strictly different from w′i because of the particle inertia and gravitational settling (crossing-trajectories effect). Reference Hunt and NalpanisHunt and Nalpanis (1985) modeled the vertical fluid velocity seen by a solid particle using Equation (3), but with a Lagrangian time-scale,
where A1 is a constant of O(1). In this study, we use Equation (4) as the Lagrangian time-scale. A1 is set to 0.5 according to Reference AndersonAnderson (1987).
The grain-borne shear stress τg is calculated from the particle trajectories as
where n↓ and n↑ show the particle number per unit area in unit time (the vertical number flux density of particles) on moving upward and downward, respectively. upi↑ and upi↓ are the horizontal speed of the ith particle moving upward and downward, respectively.
Motion of the particles is mainly determined by two forces: gravitational force and drag force due to the wind (Reference Nishimura and HuntNishimura and Hunt, 2000; Reference GauerGauer, 2001). Assuming the particle shape is spherical, the equation of motion for a snow particle is
where g is the gravitational acceleration, d the particle diameter, and up and u are particle and wind velocity, respectively. Cd is the drag coefficient, which is determined according to Reference Morsi and AlexanderMorsi and Alexander (1972). VR is the relative velocity between a particle and the wind, defined as
Particle trajectories are determined by Equation (6) provided the initial conditions and the wind field are known. To provide initial conditions for the particles, aerodynamic entrainment and splash processes are considered. The process of aerodynamic entrainment is modeled based on the excess shear stress rule, which was first proposed by Anderson and Haff (Reference Anderson and Haff1988, Reference Anderson and Haff1991). In regard to the splash process, experimentally determined splash functions expressing the probability distributions of the vertical restitution coefficient, horizontal restitution coefficient and ejection number, formulated by Reference Sugiura and MaenoSugiura and Maeno (2000), are used. Their experiment was conducted in a cold wind tunnel (–15°C) with natural compacted snow. The snow particles on the surface were not sintered. Hence, our model simulations are also restricted over the loose snow. In addition, functions are obtained at relatively low wind speed. Since the splash process can be drastically modified by the strong turbulence at higher wind speeds (Reference Nishimura and HuntNishimura and Hunt, 2000), Sugiura and Maeno’s splash functions would be applicable only at moderate wind (say 5–10ms–1 at height 1 m).The particle diameters at the lower boundary were given with a gamma random number distribution, according to the previous observations (e.g. Reference BuddBudd, 1966).
Table 1 and Table 2 show the list of the model parameters and the initial and lower boundary conditions for the wind field, respectively. In Table 1, the threshold friction velocity of drifting snow u*t and the roughness length over the snow surface z0 are determined by a wind-tunnel experiment (Reference Nemoto and NishimuraNemoto and Nishimura, 2001). In the simulation, we use two time-steps, 10–2s and 10–4 s; the latter is for the grain trajectory calculation and is two orders of magnitude smaller than the former for wind calculation, in order to adequately resolve each grain trajectory. The vertical dimension is discretized into slices, and within each slice the wind speed is calculated by the balance of the body force due to the grains and the wind shear stress (Reference McEwan and WillettsMcEwan and Willetts, 1991, Reference McEwan and Willetts1993). Detailed descriptions of the model and calculation procedures have been provided in Reference NemotoNemoto (2002) and will be published separately.
Table 1. List of the model parameters
Table 2. Initial and lower boundary conditions for the wind field
Results
Overall features of time evolution
First, calculation was conducted for an initial friction velocity u*i of 0.32 ms–1. This corresponds to a wind speed of 8.7 ms–1 at height 1 m, which is almost the same as Reference TakeuchiTakeuchi’s (1980) field observation, making intercomparison easier. Figure 2 shows the time evolution of the total mass flux over 10 s from the onset of drifting snow. During the first few seconds, the total flux Q increases rapidly, then it reaches a steady state (dQ/dt ~ 0), slightly decaying. These features were found in the wind-tunnel experiment of sand saltation by Reference ButterfieldButterfield (1991) and were also obtained by self-regulatory saltation models (Reference Anderson and HaffAnderson and Haff, 1991; Reference McEwan and WillettsMcEwan and Willetts, 1991, Reference McEwan and Willetts1993).
Fig. 2. Time evolution of the total mass flux. The initial friction velocity was 0.32 m s–1.
Figure 3 shows the changes in friction velocity, the effective roughness length and the grain-borne shear stress over time. The friction velocity and the effective roughness length were calculated by assuming a logarithmic relation for the wind profile above the saltation layer (~5 cm). As the total flux increased, the friction velocity and the effective roughness length also increased, reaching a steady state after around 5 s. The increase in the friction velocity and the effective roughness length corresponds to a change in the wind structure. This is caused by the momentum exchange between the grains and the wind near the surface, and can also be recognized from the change in the wind profile above the saltation layer. The grain-borne shear stress also increases during the first seconds and then reaches a steady state. Its trend is similar to the change in the total flux in Figure 2.
Fig. 3. Time evolution of the friction velocity (a), roughness length (b) and the grain-borne shear stress (c). The initial friction velocity was 0.32 m s–1.
Time evolution of mass-flux profile
Figure 4 shows the time evolution of the horizontal mass flux. At t = 0.6 s, the mass flux extends only from the surface to the height of z = 0.1m. The mass flux shows an increase over time. For t > 1.5 s, the mass flux below z = 0.1m does not change much. This means that the mass-flux profile near the surface has reached a steady state, which is consistent with the results of the steady-state saltation model of sand (Reference Anderson and HaffAnderson and Haff, 1991; Reference McEwan and WillettsMcEwan and Willetts, 1991, Reference McEwan and Willetts1993). In contrast, the mass flux above z = 0.1m shows a gradual increase during the calculation. It appears that at a height of several meters, it takes 50 s for the mass flux to converge into a steady profile; drifting-snow particles reached a height of 10 m. The shape of the calculated vertical profile was compared with the measurements not only in a cold wind tunnel by Reference Nishimura, Sugiura, Nemoto and MaenoNishimura and others (1998) and Reference Sugiura, Nishimura, Maeno and KimuraSugiura and others (1998) but also at Mizuho station, Antarctica, and its quantitative accuracy was verified (e.g. Reference NemotoNemoto, 2002).
Fig. 4. Time evolution of the horizontalmass-flux profile. The initial friction velocity was 0.32 m s–1.
Figure 5 shows the time evolution of the horizontal mass flux over 60 s, at four different heights. Steady state of the mass flux at each height is defined by ∂q(t; z)/∂t = 0, where q is the mass flux of drifting-snow particles at height z. At z < 0.2m, mass flux reaches a maximum value within a few seconds, after which it decreases gradually to a saturation value. However, at z >1m, the increase in mass flux is much slower. At z = 1.25m, a noticeable mass flux was observed after around 5 s, and it reached a steady state after approximately 30 s. The time to attain steady mass flux is much longer, at z = 5 m (approximately 55 s as shown in Fig. 5). These results indicate that it takes about 50 s for the flux profile of both the saltation and the suspension layer to reach a steady state in our case.
Fig. 5. Time evolution of the horizontal mass flux at four heights: (a) 5 m; (b) 1.25 m; (c) 0.16 m; and (d) 0.04 m. The initial friction velocity was 0.32 m s–1.
It should also be noted that the mass flux at z = 0.16m shows a strongly skewed variation compared to other levels (Fig. 5). Reference NemotoNemoto (2002) revealed that the particle transport mode changes from saltation to suspension around height 0.1m, which is also implied by the change in mass-flux profile in Figure 4. Above 0.1m, small particles 5100 μm in diameter are dominant, and the mass flux is several orders of magnitude smaller than it is near the surface. The height around 0.1m (0.16m in Fig. 5) corresponds to the transition zone, and larger particles 250–300 μm in diameter frequently jump into it. These particles are two orders of magnitude larger in mass and make a large contribution to the mass flux there. Thus a number of spikes appear in Figure 5c.
Figure 6 shows the 0.5 s running mean of horizontal mass flux over 90 s at five heights. It clearly indicates the vertical profile of the time for the mass flux to achieve a steady state. Although the highest case (z = 11m) still shows a slight increase even at the end of the calculation (90 s), we can reasonably say ∂q(t; z)/∂t = 0 is satisfied. In general, it takes longer to reach equilibrium with increasing height, but ∂q/∂t approaches zero at t > 50 s.
Fig. 6. Running mean horizontal mass flux (0.5 s) at five heights. The initial friction velocity was 0.32 m s–1.
The time for the mass flux to reach a steady state depends on the wind speed. Figure 7a and b show the time evolution of the horizontal mass flux with initial friction velocities of u*i = 0.23 and 0.39 ms–1, respectively. In the former case, the mass flux achieves equilibrium within 10 s, while 50–60 s are needed in the latter case. At u*i = 0.23ms–1, the wind strength is near the threshold and the suspension layer does not develop much. Thus 10 s is enough to reach equilibrium. On the contrary, at u*i = 0.39ms–1, suspension due to the turbulent diffusion is significant and the equilibrium time becomes 60 s, which is around 10 s longer than the time at u*i = 0.32ms–1.
Fig. 7. Time evolution of the horizontal mass-flux profile. The initial friction velocity was 0.23 m s–1 (a); 0.39 m s–1 (b).
Time evolution of wind speed
Figure 8 shows the time evolution of the wind profile. A characteristic modification for aeolian sediment transport processes is clearly shown: increase in effective roughness length, leading to an increase in friction velocity (Reference OwenOwen, 1964). Figure 9 shows the time evolution of the wind speed over 60 s at four different heights. The wind speed below 0.2 mdecreases immediately after the onset of drifting snow. This decrease is caused by the momentum exchange between the drifting-snow particles and the wind. In Figure 5, the mass flux near the surface shows an overshooting feature (rapid increase to a maximum value before a decrease to a steady mass flux). The overshoot phenomenon can be explained as follows. In the early stage, drifting snow develops rapidly by the increase of entrained and splashed grains. As the number of drifting-snow particles increases, part of the momentum of the wind is transferred to the particles, leading to a deceleration of the wind near the surface. This reduces the ejection rate of grains at the surface, and less momentum is transferred from the wind to the particles. This feedback mechanism leads finally to a steady state of the mass flux and the wind speed.
At a higher position, the change in wind speed is slow and small, as it takes time for the wind profile to adapt to the change in the apparent roughness length. In our case, it takes approximately 50 s to reach a steady wind profile, about as long as it takes for the mass flux at z > 1m (Fig. 5.) However, the change in the wind speed is small for heights above 1m, in particular for the wind at z = 5 m, where the mass flux is about four orders of magnitude smaller than it is near the surface.
Figure 10 a and b show the time evolution of the wind profile for different initial friction velocities. When u*i is 0.23 ms–1, the wind-profile modification is fairly small, because the wind speed is near the threshold (0.21ms–1; Table 1) and the particle concentration is rather small. The wind modification is more significant in Figure 10b than it is at u*i = 0.32ms–1 (Fig. 8), but the time needed for the wind profile to reach a steady state is approximately the same in both cases (50 s).
Fig. 8. Time evolution of the wind profile. The initial friction velocity was 0.32 m s –1.
Fig. 9. Time evolution of the wind speed at four heights: (a) 5 m; (b) 1.25 m; (c) 0.16 m; and (d) 0.04 m. The initial friction velocity was 0.32 m s–1.
Fig. 10. Time evolution of the wind profile. The initial friction velocity was 0.23 m s–1 (a); 0.39 m s–1 (b).
Discussion and Conclusions
In this study, we have shown results of time-series data obtained by a new numerical model of drifting snow. The basic concept of the model is the same as that of the self-regulating saltation model for sand (Reference Anderson and HaffAnderson and Haff, 1988, Reference Anderson and Haff1991; Reference McEwan and WillettsMcEwan and Willetts, 1991, Reference McEwan and Willetts1993). In addition, the effect of vertical turbulent wind on the particle trajectories is incorporated. Hence, the model can represent the mass-flux development within both the saltation layer and the suspension layer.
In the cases examined (u*i = 0.32 and 0.39ms–1), the time it takes for the total mass flux to reach a steady state appears to be 3–5 s. Vertical profiles of horizontal mass flux near the surface (i.e. in the saltation layer) show a similar tendency. In contrast, for the wind structure and the whole mass-flux profile, it takes >50 s to reach a steady state. This longer duration is due to the relatively long time needed to reach the final mass-flux profile, which extends to an order of a few meters high due to turbulent diffusion. In the surface layer, the eddy diffusivity coefficient Km can be expressed in terms of friction velocity and height above the surface as Km = κu*z. From this relation, we define a diffusion time-scale T(s) as
If we consider the case u* = 0.32ms–1, the thickness of the drifting-snowlayer z is approximately 10 mas shown in Figure 4. Therefore, from Equation (8), we obtain T ~ 78 s, which is of the same order as the calculated time. This consistency indicates that the turbulence has a significant effect on the dispersion and development of drifting snow.
Although the horizontal mass flux near the surface reaches saturation in the early stages of drifting snow, horizontal mass flux at z > 1m continues to increase at this stage and takes much longer to reach saturation. This result seems to contradict the fact that the total mass flux attains a steady state in a relatively short time. However, as shown in Figure 4, the mass flux in the suspension layer is more than three orders ofmagnitude smaller than the mass flux near the surface. Thus, no change is observed in the mass flux near the surface.
Time-series data of the wind speeds at u*i = 0.32ms–1 (Fig. 9) show that the wind structure also changes noticeably during approximately 50 s, particularly at high positions. Reference McEwan and WillettsMcEwan andWilletts (1993) reported that in their simulations it took >40 s for the boundary layer to respond to the development of sand saltation at u*i = 0.32ms–1. They considered that the increase in effective roughness, accompanied by the sand saltation, gave rise to the change in the wind structure at the surface and propagated upwards from there. Our results for the response time of the wind structure (~50 s) are similar to theirs, although the transition from saltation to suspension is included in our model. In the suspension layer, the momentum exchange between grains and the wind is small because the particle number flux is much smaller than it is in the saltation layer, and the particles have small mass (i.e. small diameter) and the relative velocity is small. The effect of particle inertia is presumably negligible in the suspension layer. Therefore, we can conclude that the wind modification obtained in our model is mainly induced by the saltating particles near the surface, which have large diameter and more inertia.
In this study, the time taken for the saltation and suspension layer of drifting and blowing snow to reach a steady state is found to be several tens of seconds, when the initial friction velocity u* = 0.32 and 0.39ms–1. However, some numerical models based on continuous theory (Reference Déry, Taylor and XiaoDéry and others, 1998; Reference MannMann, 1998; Reference BintanjaBintanja, 2000; Reference Xiao and TaylorXiao and Taylor, 2002) suggest a much longer time is needed, such as several tens of minutes (Reference Xiao, Bintanja, Déry, Mann and TaylorXiao and others, 2000). The discrepancy can be explained as follows. Our model incorporates the effect of turbulence on the motion of drifting-snow particles and describes both saltation and suspension processes. Thus the wind-speed modification (reduction) due to the momentum exchange between grains and the wind near the bed is also taken into account. It is important because it causes a sudden decrease in the wind near the surface (Fig. 9) and leads to the reduction of the ejection rate on the bed. On the other hand, much attention has been paid to suspension in other models; the saltation layer was given as a lower boundary condition. The wind modification above the saltation layer can be slower than it is in the saltation layer, which leads to a slow evolution of the system. The grain/bed collision process (splash function), which is omitted from other models, may affect the time-scale. It is also possible that the thermodynamic effect, which is not considered in our model (e.g. sublimation of the drifting-snow particles), lengthens the time required to achieve a steady state. As a next step, we will include the thermodynamics in our model to see how it affects the time-scale concerned.
The spatial scale of the development of the saltation and suspension layer has remained a controversial issue because there are few examples of observation; Reference TakeuchiTakeuchi’s (1980) field observation may be the only one. The model used in this study calculates the time development of vertical distribution of the mass flux and the wind; fundamentally, it is a one-dimensional model. Although the approach is physically more sound than that of other models, it is difficult to apply a vast multidimensional area of calculation with this type of model because the number of particles in the calculation is limited, as already mentioned. However, despite this restriction, applying Taylor’s frozen turbulence approximation we can indirectly estimate the spatial scale of drifting-snow development. Considering a homogeneous field in a uniform mean wind, spatial changes in a turbulence variable and its temporal changes at a fixed point can be related to each other as ∂/∂t = -U∂/∂x, where U is the mean wind velocity. The time development of the horizontal mass-flux profile is shown in Figure 4. Applying Taylor’s hypothesis, we can closely relate the time development to the spatial development of the horizontal mass-flux profile in the mean wind direction. Given that the mean wind speed at 1 m height was approximately 8 ms–1, the profiles in Figure 4 (t =0.6, 1.5, 2.5, 10, 30 and 50 s) can be considered to correspond to x = 4.8, 12.0, 20.0, 80.0, 240.0 and 400 m, respectively. Thus, we can reasonably say that the streamwise distance reaching the drifting-snow steady state is approximately 400 m.Reference TakeuchiTakeuchi (1980) observed the horizontal distribution of the total mass flux in drifting snow. His observation indicated that the total mass flux increased rapidly from the onset to 100–200m, and reached a steady state after around 200–400m for a wind speed at 1m height of 7– 8 ms–1. The estimated length scale fromTaylor’s hypothesis shows good agreement. The length scale obviously depends on the wind speed. In fact, several observations have showed that drifting snow achieves a steady state within a shorter length scale when saltation is dominant (Table 3). This feature is successfully reproduced by our model. However, for the larger wind (say u* > 1ms–1), much improvement of the model will be necessary. In particular, the splash process, which can be a function not only of the friction velocity but also of the temperature and snow properties, is important for the development of drifting snow (see Reference Naaim-Bouvet and NaaimNaaim-Bouvet and Naaim, 1998). More investigation and model improvement will increase our understanding of the drifting- and blowing-snow process and will make it possible to estimate the transport rate under different conditions.
In this study, the horizontal wind U changes in time but is assumed to be uniform within a small domain. In reality, however, drifting snow develops in the mean wind direction, and the horizontal advection term (∂/∂x) would not be zero; U is a function of both time and space. Reference Shao and LiShao and Li (1999) reported a two-dimensional saltation model for sand. Their results indicate U(x) decreases with x from the starting point of saltation, and as far as the saltation layer is concerned, equilibrium saltation is achieved at a distance around 15 m, which corresponds to t < 1 s. The reduction of U in the mean wind direction would shorten the equilibrium time. The extension of our model to two dimensions is the logical next step to investigate precisely the time and length scale of drifting snow.
Table 3. Observed distance for saltation layer to reach equilibrium
Acknowledgements
The authors gratefully acknowledge helpful discussions withT. Ohata, N. Maeno, T.Yamada and T. Sone from the Institute of LowTemperature Science, Hokkaido University. Thanks are also due to T. Sato of the Nagaoka Institute of Snow and Ice Studies, Shinjo Branch for valuable comments and suggestions. We also gratefully acknowledge the assistance of I. K. McEwan of the University of Aberdeen, Scotland, with the numerical model.
