Notation
1. Introduction
In Alpine terrain, redistribution of snow by wind is a major factor in snow pack evolution as well as in the formation and effectiveness of avalanches. Thus, it is of great interest to predict snow redistribution and its influence on avalanche danger with regard to avalanche forecasting and landuse planning.
Avalanche warning on a regional scale requires information on the occurrence of snowdrift and on the (new) snow distribution on a mesoscale (several 10 km^{2}). Scarcity of land resources is a major restricting factor in the economic development of many mountain villages. In these villages, avalanche danger often causes restrictions on land use. As great interest exists in using all available land resources, efforts are made to develop numerical avalanche models to reliably determine endangered zones. These models require highly resolved input data on the snowmass distribution in avalanche release zones (microscale; several 100 m^{2}) which depends strongly on snowdrift. In both cases, one is interested in suitable tools to assess, analyze or forecast the effect of blowing and drifting snow in Alpine topography. Corresponding to World Meteorological Organization (WMO) standard observations, the term drifting snow is used to describe the nearsurface wind transport, whereas the expression blowing snow is reserved for situations where particles rise to a height of 1.8 m and above. However, in the following the term drifting snow or wind transport will be used to refer synonymously to both phenomena.
The snowdrift model presented here is based on many previous field studies of snow transport as well as on physical and numerical simulations for aeolian particle transport. The field studies show that snow transport is basically given by the transport in saltation and in suspension (Reference KobayashiKobayashi, 1972; Reference RadokRadok, 1977; Reference TakeuchiTakeuchi, 1980; Reference SchmidtSchmidt, 1982, Reference Schmidt1986; Reference Mellor and FellersMellor and Fellers, 1986). These studies were carried out on plains and aimed at determining the transport rate for steadystate conditions. The attained empirical formulations show that transport rate is proportional to a power of the wind speed, where the power varies typically between 2 and 4. Only a few attempts have been made to measure snowdrift in mountain areas (e.g. Reference FohnFöhn, 1980; Reference Schmidt, Meister and GublerSchmidt and others, 1984; Reference MeisterMeister, 1987) where drift flux measurements were made on a crest. Reference Fohn and MeisterFöhn and Meister (1983) also measured the spatial distribution of snow in the area surrounding a crest. Similar measurements were made by Reference CastelleCastelle (1995) for a mountain pass. Physical simulations have been used to study the deposition patterns around obstacles (e.g. Reference IversenIverson, 1980; Reference AnnoAnno, 1985) or have dealt with particular features of the drift process (e.g. the particle–bed impact as studied by Reference Maeno, Araoka, Nishimura and KanedaMaeno and others (1979, Reference Maeno1985, Reference Maeno, Nishimura, Sugiura and Kosugi1995), Reference Willets, Rice and BarndorffNielsenWillets and Rice (1985), Reference Rice, Willets and McEwanRice and others (1995, Reference Rice, Willets and McEwan1996) and Reference Nishimura, Sugiura, Nemoto and MaenoNishimura and others (1998). Numerical models of snow or sand drift differ widely in scope and focus. Specific processes have been studied, for example, simulating the grain impact (Reference Werner and HaffWerner and Haff, 1988) or particle transport in saltation (Reference Anderson and HaffAnderson and Haff, 1988, Reference Anderson and Haff1991; Reference WernerWerner, 1990; Reference McEwan and WilletsMcEwan and Willets, 1991; Reference SorensenSørensen, 1991). Also models based on heuristic rules for avalanche forecasting have been developed, like ELSA (Reference Mases, Buisson, Frey and SivardièreMases and others, 1995).
Most transport models for snow can be classified into two categories: Eulerian–Lagrangian (e.g. Reference Sato, Uematsu, Kaneda, Izumi, Nakamura and SackSato and others, 1997; Reference Sundsbø, Hansen, Izumi, Nakamura and SackSundsbo and Hansen, 1997) or Eulerian–Eulerian (Reference Uematsu, Kaneda, Takeuchi, Nakata and YukumiUematsu and others, 1989; Reference Liston, Brown and DentListon and others, 1993). However, Reference Masselot and ChopardMasselot and Chopard (1998) follow a different approach, using a lattice Boltzmann model. Particle tracking, which is done in the Eulerian–Lagrangian approach, can require large computational effort and restrict the application of those models to small areas (e.g. the estimated number of particles in the saltation layer is about 1–5 × 10^{6} m^{−2}). All the models discussed above are twodimensional. More recent models address the requirements of threedimensional transport modeling (Reference Liston and SturmListon and Sturm, 1998; Reference Naaim, NaaimBouvet and MartinezNaaim and others, 1998). These models are Eulerian–Eulerian. Common to most transport models is the separation between saltation and suspension. Either the transport rate for suspension is defined by a height integration of the flux with given profiles for the wind speed and concentration (Reference PomeroyPomeroy, 1989; Reference Liston and SturmListon and Sturm, 1998), or a separate balance equation for the particle concentration is solved (Reference UematsuUematsu, 1993; Reference Naaim, NaaimBouvet and MartinezNaaim and others, 1998). The assumption that the wind profile develops fully into a logarithmic profile over the total height of the suspension layer, which Reference PomeroyPomeroy (1989) and Reference Liston and SturmListon and Sturm (1998) take for granted, is critical in Alpine terrain. For example, measurements (Reference FohnFöhn, 1980) show that the wind profile on a crest is very different from a logarithmic profile. This difference also affects the concentration profile. Almost all descriptions of the transport in saltation (e.g. Reference PomeroyPomeroy, 1989; Reference UematsuUematsu, 1993; Reference Liston and SturmListon and Sturm, 1998; Reference Naaim, NaaimBouvet and MartinezNaaim and others, 1998) rely on Reference PomeroyPomeroy’s (1989) empirical formulation, which is defined for steadystate conditions. Reference Pomeroy and GrayPomeroy and Gray (1990) assume that it takes about 300–500 m for a boundarylayer flow that is approximately 3 m in depth to reach steady conditions. In Alpine terrain, with its high variability, 300–500 m is a long distance. Reference CastelleCastelle (1995) pointed out that the drift transport seldom if ever reaches steady state in mountainous terrain. Recently, attempts have been made to parameterize the deviation from steady state (Reference Liston and SturmListon and Sturm, 1998; Reference Naaim, NaaimBouvet and MartinezNaaim and others, 1998).
Due to the limited range of applicability because of the simplifications that were used, the above models cannot accurately represent complex topographies like an Alpine terrain and nonsteady wind fields. To overcome some of the weaknesses of the previous models, the model presented here uses a different approach for the description of saltation. Based on particle motion, a continuum formulation is derived for the saltation layer. This formulation includes deposition and erosion, distinguishing between aerodynamic entrainment and particles ejected due to impact, and does not require steadystate conditions.
2. Numerical model
Blowing and drifting snow is strongly related to the atmospheric boundary layer (ABL). Erosion, transport and deposition of snow depend on the wind field and the turbulence within the lowest 10–100 m of the boundary layer, the socalled surface layer (SL). It is widely accepted that the primary transport mechanisms are saltation and suspension. The transport due to saltation usually starts at light wind speeds, M _{10} = 5–8 m^{−2}, where M _{10} is the wind speed at 10 m (Reference SchmidtSchmidt, 1980; Reference CastelleCastelle, 1995; Reference Li and PomeroyLi and Pomeroy, 1997). Typical transport rates (integrated over the height) range from 0.001 to 0.03 kg m^{−1} s^{−1} (Reference TakeuchiTakeuchi, 1980; Reference CastelleCastelle, 1995). Saltation is usually restricted to a vertical extension of about 0.01–0.1 m (Reference KobayashiKobayashi, 1972), and drift densities range from 0.1 to 1 kgm^{−3} (Reference Mellor and FellersMellor and Fellers, 1986; Reference Pomeroy and GrayPomeroy and Gray, 1990; Reference CastelleCastelle, 1995). For comparison, even if the transition from saltation to suspension is more or less continuous, noticeable suspension starts at moderate wind speeds (M _{10} = 7–11 m^{−2}), depending on the particle size and flow turbulence. An estimation of the required shear velocities is given in Reference GauerGauer (1999). The typical mean grainsize in drifting snow is about 150–200 µm (Reference SchmidtSchmidt, 1984; Reference PopePope, unpublished). Transport rates for suspension range from 0.01 to 0.1 kg m^{−1} s^{−1} and drift density is <1 kg m^{−3} (Reference Mellor and FellersMellor and Fellers, 1986; Reference CastelleCastelle, 1995). The vertical extension can reach several tens of meters, though the drift density usually decreases significantly with increasing height above the surface.
For a review of current conceptual understanding of aeolian particle transport, the reader is referred to Reference Anderson, Sørensen and WilletsAnderson and others (1991) or Reference McEwan, Willets and PyeMcEwan and Willets (1993). In a greatly simplified form, one can imagine the transport due to blowing and drifting snow as follows (see also Fig. 1): If the wind blowing over a snow surface becomes sufficiently strong, and wind shear exceeds a certain critical value, the socalled threshold, some grains are set in motion by the wind. Initially, only a few of these will be lifted off the surface and accelerated by the wind. Some of them will gain enough energy to rebound and/or eject other grains on impact. At the onset, the number of grains resulting from an impact, the socalled mean replacement capacity, is more than one on average. This results in an exponential increase of grains in the saltation layer, which follow more or less ballistic trajectories determined by the timeaveraged wind profile. Grains within the saltation layer are not — or only weakly — influenced by turbulence of the wind. As more and more grains saltate, the vertical wind profile changes due to the considerable extraction of momentum from the airflow by the grains in motion. Now the grains gain less energy, and the probability of a grain rebounding and/or ejecting other grains on impact decreases until the mean replacement capacity reaches the equilibrium value of one. At this stage, the number of saltating grains fluctuates around a certain value, sometimes called the saturation value, which depends on the drivingwind field and certainly on the properties of the snow surface. The surface properties determine the fluid threshold as well as rebound and dislodgements effected by collisions of grains with the surface, and thus play an important role in saltation. In turbulent and gusty winds, a certain number of particles will be caught by eddies and will travel a significantly greater distance without surface contact. The ratio between mass transport by suspension and mass transport by saltation increases with increasing wind speeds and turbulence. In all cases, gravity counteracts the wind forces on the grains in saltation and suspension.
In short, wind transport can be regarded as the result of five closely/mutually linked processes:

aerodynamic entrainment,

grain trajectories,

grain–snowpack impacts,

modification of the wind field,

transport due to turbulent suspension.
In the following, a numerical model is proposed that uses a continuum mechanical approach to describe blowing and drifting snow, and makes use of a splitting of the two dominant transport processes into two mutually coupled layers. This procedure is chosen to reduce the computational effort, which would drastically increase if the saltation layer was to be fully resolved. The transition mode between saltation and suspension, the socalled modified saltation, is not separately considered and it is assumed that those particles are either in the saltation layer or in suspension. The complete derivation for the model is given in Reference GauerGauer (1999).
2.1. Windfield and suspension modeling
As outlined above, the driving force for winddrifting snow is the airflow within the surface layer. This flow can be described using the equation of state (1), and the conservation equations for mass, momentum (3), moisture (4) and heat (5). For atmospheric flows smaller than mesoγscale (≪ 12 km), the conservation equation for mass usually reduces to the incompressibility approximation (2). Note that Einstein summation notation will be used in the following. In Equations (1–5), all unprimed quantities are to be understood as Reynolds averaged quantities.
In Equation (4), q _{T} denotes the total specific humidity of air, which can be split into a vapor and a nonvapor part, using q _{T} = q + q _{L}. The heat balance (Equation (5)) is expressed in terms of the potential temperature, . The virtual temperature, T _{v}, is defined as the temperature at which dry air must be in order to have the same density as moist air at the temperature T. It is given by
where r is the mixing ratio, and r _{sat} and r _{L} are the saturation and solid/liquid water mixing ratios, respectively. By substituting Equation (6) into Equation (1), effects of mass concentration on the wind flow can be taken into account. For a complete derivation of these equations, the reader is referred to the literature (Reference PichlerPichler, 1986, Reference StullStull, 1997). In the following, it will be assumed that temperature effects, sublimation and effects due to the Coriolis term, f _{c} ∊_{ij} _{3} u_{j} , are negligible for the present purpose of the model. It is assumed that snow grains travel with a velocity , where W _{f} is the absolute value of the freefall velocity of a grain. This assumption implies that all grains have the same size. On the conditions above, Equation (5) is negligible, a fixed reference temperature, T _{ref}, can be used, and Equation (4) can be replaced by the balance equation for only the solid part,
to describe precipitation as well as the suspension transport. The solid part of the moisture can be written as q _{L} = cρ _{p} /ρ _{a}, where c is the volume fraction of snow and ρ _{p} is the intrinsic density of ice. On the righthand side of Equation (7), a sink/source term for the phase change due to sublimation is neglected. The sublimation rate of an ice particle is determined by the difference between the saturation vapor pressure over ice, as a function of the particle radius, and the watervapor pressure of the ambient air. At the present stage, the model will be used to represent the redistribution of snow during precipitation events with a relative humidity of approximately 100%. In this case, it is reasonable to assume that the difference in the vapor pressures is small and hence the sublimation loss negligible. Otherwise, the sublimation rate must be parameterized (as in Reference Liston and SturmListon and Sturm, 1998), or in a consistent manner the phase exchange, −E/ρ _{a}, must be included on the righthand side of Equation (7) and an additional equation for the vapor part is needed. As phase exchange depends on the grainsize, the grainsize distribution should be considered. Also the latentheat transfer should be taken into account, so that Equation (5) is no longer negligible. The particle temperature, which determines the saturation vapor pressure, is influenced by the absorbed solar radiation. This fact as well as Equation (5) require an approximation of the incoming solar radiation.
For turbulent closure, the wellknown e–∊ model is used (see, e.g., Reference RodiRodi, 1980).This was first used in a snowdrift model by Reference Liston, Brown and DentListon and others (1993) in their twodimensional model, and has been used to simulate boundarylayer evolution, flow changes in roughness and topography, and seabreeze fronts (Reference StullStull, 1997). It is a oneandahalforder closure and is used as a compromise between accuracy and effort. Here, Reynolds stress and Reynolds flux are approximated by
and
where the turbulent viscosity is parameterized by
Separate balance equations are formulated for the turbulent kinetic energy, e, and the dissipation rate, ∊, i.e.
and
The empirical constants are set to c_{µ} = 0.09, c _{1} _{∊} = 1.44, c _{2} _{∊} = 1.92, σ _{T} = 0.9, (σ _{e} = 1.0 and σ_{∊} = 1.3 (Reference RodiRodi, 1980).
The description of the suspension layer above is similar to simple mixture formulations used in the description of passive tracers and in airpollution modeling. As in those models, momentum transfer between the different phases is assumed negligible. However, mass change due to the suspension of snow is included here. Reference Naaim, NaaimBouvet and MartinezNaaim and others (1998) used a slightly modified e–∊ for their snowdrift model to account for the effect that particles have on turbulence.
2.2. Saltationlayer modeling
Within the saltation layer, the assumption of a negligible momentum transfer breaks down, as the acceleration of massive grains exerts an additional force on the wind. Hence, within the saltation layer, the mass balance and the momentum balance are formulated for each of the two phases, dry air and snow. Scaling analysis (Reference GauerGauer, 1999) shows that mass conservation within the saltation layer is primarily determined by the conservation of the snow mass
and that the reduced mixture momentum equation describes the balance of the force necessary to accelerate the saltating particles and the driving forces, represented by the Reynolds/ shear stress and gravity. Therefore,
where ρ _{d} is the density of dry air. In the numerical model, the following discretized volumeaveraged versions of Equations (13) and (14) are used
and
where sub/superscripts 〈1〉, 〈2〉, 〈3〉, 〈4〉, 〈5〉 and 〈6〉 stand for the direction of the control volume in x, y, z, −x, −y and −z, respectively, and are the approximated slopeparallel acceleration acting on a rebounding or ejected grain, respectively, during its jump, and c _{R} and c _{E} are the corresponding parts of the volume fraction of snow, c, within the saltation layer. Let us assume that the ratio of rebounding and ejected grains in a control volume will remain unchanged due to advection during a timestep. In this case,
and
where and are the calculated particle concentrations of rebounding and ejected particles, respectively, from the previous timestep. The advection coefficients are given by
where , denotes the vector of the mean particle velocity and are the area vectors of the sides in the direction k of the control volume. The volume V is approximated by A^{〈6〉} h _{s}, where A ^{〈6〉} is the surface area and h _{s} is the saltation height. A parameterization of ,h _{s}, , etc., will be defined later on.
2.2.1. Concentration exchange between saltation and suspension layer
The settling due to gravity, J _{cs}, and the turbulent entrainment, J _{ct}, determine the boundary conditions of the concentration exchange between the saltation layer and the suspension layer. Grains settle out of the suspension layer through the boundary area, A^{〈3〉}, at velocity . Turbulent fluctuations will also influence the settling. On the other hand, grains in the saltation layer can be caught by turbulent eddies and lifted off, and then transported by the mean flow. For the description of the mass exchange the following statement is used:
where A _{〈3〉} is set to
and c _{〈3〉} is given by
and the diffusion coefficient describing the turbulent dispersion is set to
n_{i} is the normal vector to the boundary area between the saltation layer and the suspension layer, is the area vector. , is the wind velocity at the saltationlayer height and c _{s} is the volume fraction of snow in the suspension layer above. The value is obtained by similarity theory for the neutral boundary layer and corresponds to the vertical fluctuation velocity (Stull, 1987). The statement for the mass exchange is adapted from a hybrid differencing scheme in which central differencing is used if the mesh Peclet number (C/D) is < 2, and upwind differencing, neglecting diffusion, is used if the mesh Peclet number is > 2.
2.2.2. Erosion and deposition
Although Reference BagnoldBagnold (1954), in one of the first thorough examinations of the drift processes, recognized the difference between fluid threshold and impact threshold, most models do not distinguish between aerodynamic entrainment and ejecta due to impact. Reference Anderson, Sørensen and WilletsAnderson and others (1991) as well as Reference McEwan, Willets and PyeMcEwan and Willets (1993) emphasized the importance of the impact process for particle entrainment, and it is observed that threshold wind speeds are significantly reduced if there are already snow grains in the air.
The mass exchange between the snowpack and the saltation layer is determined by erosion and deposition of snow grains. As mentioned above, two different mechanisms are important for erosion: aerodynamic entrainment and particle impact. At present, little is known about the exact nature of either mechanism in the case of snow. In the absence of better knowledge on aerodynamic entrainment, the statement proposed by Reference Anderson and HaffAnderson and Haff (1991) is adopted, i.e. a linear relation between the number of entrained grains per unit area and unit time, N _{ae}, and the excess shear stress is used. Thus,
where τ _{a} is the socalled airborne shear stress and is assumed to be equal to the shear stress of the airflow at the bed, τ _{c} is the critical fluid shear stress for entrainment and ζ is an empirical constant. As no value is known for the case of snow, ζ = 10^{5} grains N^{−1} s^{−1} is used, which is the value used by Anderson and Haff for sand. Actually, the choice of ζ is not that important, as fluid entrainment essentially acts during the initial development of wind transport. Once there are some particles in transport, they act as seeding agents for further dislodgements. The impacts between these particles and the surface eject more particles, leading to a cascade as mentioned above. This cascade drives the transition from an aerodynamically controlled to an impactcontrolled system. More important is the choice of τ _{c}, which determines the onset of drift, at least during periods without additional precipitation. τ _{c} is highly dependent on the snowpack properties and meteorological conditions such as particle size and shape, bonding, temperature and humidity. Reference SchmidtSchmidt (1980) tried to estimate τ _{c} for wind conditions over a flat snow surface. According to him, values for τ _{c} range from 0.04 to 9000 Pa. In comparison, typical values for the wind shear stress/Reynolds stress, τ _{a}, during snowdrift vary from 0.06 to 1 Pa, indicating that aerodynamical entrainment is only possible within a small range. The lower end of the estimation above is in agreement with observations based on the onset of drift (see, e.g., Reference CastelleCastelle, 1995; Reference Li and PomeroyLi and Pomeroy, 1997). Little is known about the influence of snowpack properties and meteorological conditions on τ _{c}. Reference Li and PomeroyLi and Pomeroy (1997) show a relation between air temperature and the mean threshold wind speed. An attempt to relate the threshold to snowpack properties, such as particle size, sphericity and dendricity, can be found in Reference Guyomarc’h and MérindolGuyomarc’h and Mérindol (1998).
The erosion due to particle impacts is determined by the number of particle impacts per unit time on the snow surface and the number of ejected particles per impact. The impact rate is given by
where V _{p} is the particle volume and t _{J} is the average duration of a particle jump. To estimate the number of particles ejected per impact, N _{EpI}, it is reasonable to assume that this number is a function of the difference of the kinetic energy of a particle before the impact, E _{I}, and after it, E _{R}, of the total kinetic energy of the ejected particles, N _{EpI} E _{E}, of the energy dissipated into the bed during the impact, E _{D}, and, if cohesion is included, of the total energy of the bonds, N _{EpI} E _{B}, which had to be broken. From the energy balance, we obtain
All these energies can be expected to depend on the mechanical properties of the snowpack, but the exact relationship is as yet unknown. To estimate a value for N _{EpI}, some assumptions are made. The assessment for E _{D} was based on the numerical simulations of Reference Anderson and HaffAnderson and Haff (1991) who found that, in the case of sand (E _{B} = 0), the mean number of ejecta increases approximately linearly with the impact speed, where the ejection rate N _{EpI} ≈ 1.5 m^{−1} s U _{I}. For the assessment of E _{I}, E _{R} and E _{E} see Table 1. Figure 2 shows an approximation of the ejection rate with E _{B} as parameter. The range of E _{B} from 10^{−10} to 10^{−8} J seems to be reasonable regarding different snow types (Reference SchmidtSchmidt, 1980; Reference GublerGubler, 1982). Using Equations (26) and (27), the ejection rate can be expressed as
Here and in Equations (29) and (30), there is a distinction between rebounding particles and those particles which were just ejected because their respective impact speed U _{I} and jump time t_{J} are significantly different. Both kinds are marked by subscript ()_{R} and ()_{E}, respectively.
Snow grains in saltation tend to rebound at the surface. On the other hand, if a snow grain does not rebound it will be deposited at the surface. Hence, if P _{R} is the probability that a snow grain will rebound on impact, 1−P _{R} is the probability that it will be deposited. Thus, the deposition rate can be described by
Physical simulations (Reference Kosugi, Nishimura and MaenoKosugi and others, 1995; Reference Rice, Willets and McEwanRice and others, 1995, Reference Rice, Willets and McEwan1996) show that the probability of a particle rebounding is approximately 0.6–0.9 for the investigated impact speeds. In the absence of better data, the relation
is used, which Reference Anderson and HaffAnderson and Haff (1991) proposed for sand. Assuming that the rebound probability for a snow grain at the surface of a snowpack is smaller than for a sand particle on a bed of sand, the coefficient ξ is set to 1.5 m^{−1} s instead of the original value, 2 m^{−1} s. Actually, the rebound probability is expected to depend on the snowpack properties, but little is known about the relationship as yet.
2.2.3. Modification of the wind field
Reference BagnoldBagnold (1954) was the first to note the feedback between particle movement and the driving wind. Thereafter, it became common to parameterize the interaction between the particles in the saltation layer (extraction of momentum) and the wind above, using an increased roughness length, z _{0}, in the equation for the logarithmic wind profile
Reference OwenOwen (1964) assumed that the saltationlayer height scales as and that the roughness length can be given by
However, there are two weak points. First, it is unreasonable to assume that the roughness length depends only on the saltationlayer height and that there is no dependency on the number of particles within the saltation layer. This oversimplification might explain why there is a large difference in the reported values for the constant k. Owen proposed k = 0.02, whereas Reference Rasmussen, Sørensen, Willets and BarndorffNielsenRasmussen and others (1985) reported values of 0.14 and 0.18, and Reference Pomeroy and GrayPomeroy and Gray (1990) set k to 0.1203. Second, Equation (31) cannot represent the rapid decrease in the wind speed just above the top of the saltation layer (h _{s} ≈ 1–2 cm), an effect which can be observed, for example, in Bagnold’s classical windvelocity plot (Reference BagnoldBagnold, 1954, fig. 17).
The primary effect that the particles have on the airflow is to provide a spatially distributed momentum sink, with sink strength dτ _{a}(z)/dz ≅ −dτ _{g}(z)/dz per unit volume of fluid. τ _{a}(z) is the flux of fluid momentum across the level z, and τ _{g}(z) represents the net flux of slopeparallel grain momentum across that level. dτ _{g}(z)/dz corresponds to the interaction (drag) forces on the wind, due to the slopeparallel acceleration of the grains by the air. This acceleration mainly takes place within the saltation layer. The total shear stress is given by
where is the stress at the top and above the saltation layer. Vegetation canopies have a similar effect on flow. In boundarylayer meteorology, it is customary to describe the effect of vegetation on the flow above by a generalized logarithmic law (Reference ThomThom, 1971; Reference JacksonJackson, 1981; Reference WieringaWieringa, 1993). In such a case,
where d is the displacement height and z _{0} is the roughness length. According to Reference ThomThom (1971) and Reference JacksonJackson (1981), z _{0} expresses the magnitude of forces acting on the surface, whereas d is related to the distribution in z of these forces. Hence, for a known distribution of the interaction (drag) forces within the saltation layer, F _{DW}(z), the displacement height can be calculated according to the statement of Reference Shaw and PereiraShaw and Pereira (1982):
where and τ _{a} are the airborne stresses at the top of the saltation layer and at the surface of the snowpack, respectively. Assuming that the particles do not influence each other, the interaction force can be written as
where N(z) is the number of particles in a height z with a slope parallel to the relative velocity . C _{D} is the drag coefficient and d _{P} is the grain diameter. In Equation (37), only the acceleration of the grains due to the wind is represented. The reverse transfer of momentum due to descending grains from faster flow levels causes an enhancement of the air velocity close to the surface. This departure from the logarithmic wind profile within the saltation layer is observed, for example, by Reference BagnoldBagnold (1954) and is also reported by Reference Maeno, Araoka, Nishimura and KanedaMaeno and others (1979). For the parameterization, this departure is taken into account, as discussed in the following subsection.
For the description of the roughness length, the following statement is used:
where z _{0s} is the aerodynamic roughness length for a snow cover without saltation, typically in the range 10^{−4} to 10^{−3} m (Reference StullStull, 1997, p. 380), here set to 3 × 10^{−4} m, and, in the absence of data for snow transport, λ is assumed to be 0.1, a value used for vegetation. Equation (34) is used as a boundary condition for the wind speed, at the interface between the saltation and suspension layers (see section 2.4), assuming that, at least in the lowest part of the suspension layer, a logarithmic profile exists.
2.3. Parameterization
To close the system of equations for the saltation model, h _{s}, d, , , etc., must be determined.
The transport within the saltation layer is regulated by the motion of saltating grains, bouncing off the snow surface and following ballistic trajectories. It is widely accepted that the saltation trajectories are governed by the following relation:
where the terms on the righthand side denote the gravitational force (I), the drag force (II), the lift force (III), Magnus force (IV), added mass force (V), pressuregradient term (VI), electrostatic force (VII) and the Basset history term (VIII) (for reference see, e.g., Reference Clift, Grace and WeberClift and others, 1978, ch.11; Reference Schmidt, Dent and SchmidtSchmidt and others, 1998). With the aid of an appropriate scale analysis (Reference GauerGauer, 1999), it can be shown that, in a first approximation, Equation (38) reduces to a balance of inertial, drag and gravity forces, such that:
where U_{ϱi} is the relative velocity between the particle and the airflow (U_{ϱi} = U_{Pi } – u_{i} ). These equations form a system of differential equations which can be solved for given initial conditions, U_{P} _{0i }, and a known wind field, u_{i} . They build the basis for parameterization of the saltation layer.
Physical simulations with sand and snow (Reference Kosugi, Nishimura and MaenoKosugi and others, 1995; Reference Rice, Willets and McEwanRice and others, 1995, Reference Rice, Willets and McEwan1996; Reference Nishimura, Sugiura, Nemoto and MaenoNishimura and others, 1998), as well as numerical simulations (Reference Werner and HaffWerner and Haff, 1988; Reference Anderson and HaffAnderson and Haff, 1991) of the particle–bed impact, have been carried out. All these studies confirm that impacting particles behave in similar ways, in that they have a certain probability to rebound, and that the mean impact, launch angles, and ratio of impact to rebound velocity vary only over a small range. Ejected particles also behave in a specific manner. Typical values of the quantities describing the particle–bed impact are shown in Table 1 (for notation see Fig. 3).
ς _{R}, ς _{E}, α _{R}, and α _{R} describe the particle–snowpack interaction on impact (e.g. ς _{R} characterizes the energy loss of a rebounding grain on impact). All these parameters are expected to be affected by the mechanical properties of the snowpack, but little is known about the relationship as yet. For our parameterization, the following values are chosen: ς _{R} = 0.5, ς _{E} = 0.1, α _{R} = 40° and α _{E} = 65°. Furthermore, within the saltation layer the similarity assumptions are made, i.e. the airflow velocity, u_{i} , and the concentration, c, maintain approximately similar profiles in the z direction (perpendicular to the surface) as they change in time or in the x and y directions. In particular, it is assumed that
where ξ_{u} = [1 + ς(1 − z/h _{s})]^{−2} with ς = 0.35, and ξ _{c} = 1 is assumed. Based on the assumptions above and on Equation (39) as well as on values given in Table 1, it is possible to calculate the trajectory of a single grain. To be consistent, the determination of the wind field within the saltation and the particle trajectories should be carried out simultaneously. The use of the profile function (Equation (40)) is a simplification, where it is assumed that the momentum extraction and the modification of u_{i} , due to the particles is similar to the case of vegetation. Figure 4 shows an example of such a calculation. Here, the normalized horizontal velocity of a rebounding particle is presented as a function of the normalized height for the ascending as well as the descending branch of the trajectory. The mean horizontal velocity, , and the assumed wind profile within the saltation layer are also shown. The mean particle velocity is approximately . Calculations with varying show similar results. In the saltation model, is set to 0.58 . The other missing parameters of the particle trajectories (e.g. h _{s}, d, , etc.,) can also be determined in this way. Figure 5 presents those parameters as functions of for a rebounding grain. Similar calculations were used to determine the corresponding quantities for ejected grains.
2.4. Implementation
The model is implemented in the commercial flow solver CFX4.1 from AEA Technology, England (AEA Technology, 1995). This program is based on the finitevolume technique and uses a bodyfitted grid which can be modified during the calculation. Thus the modification of the topography due to deposition and erosion of snow can be modeled. Socalled Fortran user routines provide the programming interface for problemspecific modifications of the code and have proved flexible enough to accommodate all the extra computations for the simulation of blowing and drifting snow. Hence, the program serves as a framework, providing the solver for the wind field and for the suspension mode, whereas the calculations of the saltation mode were embedded in Fortran routines. Although the wind field, suspension and saltation form a mutually coupled system, here a certain decoupling was desired. The wind field and the suspension mode are solved together and provide the boundary condition for the saltationlayer calculation, which is carried out between two timesteps. On the other hand, the boundary conditions for the next timestep were obtained from this calculation. In order to evaluate the snowpack, the grid is adjusted to the newsnow depth at regular time intervals. For that purpose, at every timestep ΔHS _{i} is calculated according to
where ΔHS_{e} is the remaining erodible snow depth and c _{λ} is the volume fraction of snow in the snowpack. After a certain number of timesteps, n, the grid is adjusted to the newsnow depth:
3. Field measurements and simulation results
To test the model, two methods were pursued. Firstly, the model was compared with experimental results described in the literature. Secondly, it was tested against data from the Gaudergrat experimental site in the Swiss Alps.
Reference KobayashiKobayashi’s (1972) are one example of measurements reported in the literature. Figure 6 shows a comparison between his field experiments and drift simulations. The mass flux, Q, is plotted against the wind speed at 1 m. Kobayashi noted that due to the scattering in the drift rate, especially at low wind speeds, no distinct threshold wind speed for the occurrence of snowdrift was found, but it was in the range 4–6 m s^{−1}. For the simulation, τ _{c} was set to 0.05 Pa, and E _{B} was assumed to be 10^{−9} J. Taking into account the uncertainties in the snowpack properties, the simulation and the measured data are in good agreement.
A logarithmic plot of the simulated steady mass flux vs M _{0}._{25} is given in Figure 7. At low reference velocities, a rapid increase in the mass flux can be observed. Then the rate of increase diminishes with increasing M _{0.25}, but is still nonlinear. The power decreases from approximately 4 to 2 with increasing wind speed. This decrease indicates that with increasing wind speed a saturation of transport occurs. The power range agrees quite well with those of empirical formulations given in the literature (Reference KobayashiKobayashi, 1972; Reference TakeuchiTakeuchi, 1980; Reference SchmidtSchmidt, 1982, Reference Schmidt1986; Reference Nishimura, Sugiura, Nemoto and MaenoNishimura and others, 1998). The fit in Figure 7 corresponds to a power of approximately 3, an oftfound value in empirical formulations. A separate comparison between measured and simulated drift rates in saltation shows that the model reproduces these rates quite well (Reference GauerGauer, 1999).
In the next step, the model was applied to the Gaudergrat experimental site for a fullscale simulation of an Alpine crest. Gaudergrat ridge, 2 km north of Weissfluhjoch/Davos, has a rather sharp crest — the slope angles range from 28° to 38° — and might be regarded as prototypical of Swiss Alpine topography. The prevailing wind direction during strong precipitation periods is more or less northwest and thus perpendicular to the crest line. Up to now, field measurements of drifting snow in complex terrain are rather scarce, and, although the wind as driving force is of major importance, no example is known where the wind field and the snow distribution have been simultaneously measured. In winter 1996/97, the ridge was equipped with five wind masts in the surrounding area (M78, M76, M75, M74/73 and M72) to obtain a good impression of the wind field around the crest during interesting snowdrift periods. In addition, meteorological and snowpack parameters were measured on both sides of the ridge (M77 and M74/73). Figure 8 shows the locations of the installed masts and the slope angles at the ridge. A complete description of the experimental site and the measurements can be found in Reference GauerGauer (1999).
For comparison with the simulations, the spatial snow redistribution was measured by soundings of the snow depth before and after a drift episode along five equidistant lines across the crest, roughly 8.5 m apart and > 200 m long. The soundings were taken at 4 m intervals along the lines; in the neighborhood of the crest this distance was reduced. Thus, every field campaign resulted in around 350 data points. All measurement points were marked by thin bamboo stakes during the first campaign, so later measurements could be taken at the same points (with an estimated horizontal deviation of < 0.1 m). In this manner, uncertainties in the depth measurements due to the smallscale topography were minimized. The settling of the snowpack, different for south and north aspects, between two successive soundings was taken into consideration. For this purpose, small pits were dug at several points of both slope sides, and the snow density was measured, and all newsnow depths are recalculated for the average newsnow density (110 kg m^{−3}). For the comparison with the simulation these newsnow depths are normalized with HN_{ref}, where HN_{ref} is the measured newsnow depth at the nearby Swiss Federal Institute for Snow and Avalanche Research study plot Weissfluhjoch (VF) for the corresponding periods. Correspondingly, in the simulation, HN_{ref} is the simulated newsnow depth unaffected by the wind. During winter 1996/97, six strong precipitation periods with snowdrift were investigated.
For the fullscale simulation of the Gaudergrat site, a horizontal area of 1000 × 600 m^{−2} was chosen around the crest line (see Fig. 8). The grid consists of 120 × 37 × 25 cells for the flow domain and 120 × 37 cells each for the saltation layer and for the snowpack. The grid spacing Δx ranges from 15 m at the inflow and the outflow to approximately 2.5 m close to the crest line, and Δy ranges from 20 m at the sides to approximately 8 m in the middle. In the flow domain, the spacing in the z direction varies between approximately 2 m at the bottom and 35 m at the top. The timestep was set to 0.25 s.
The upwind boundary conditions were derived from wind measurements at mast M78, and the precipitation rate was determined by
where RS(t) is the measured 10 min precipitation sum at VF. All given boundary conditions were allowed to vary with the time. The initial snow depth was set to 0.2 m of erodible snow. During the simulation, the evolution of the snowpack was followed so that it was possible to decide if new snow was deposited or old snow eroded. For the simulation, the critical shear stress of the old snowpack was set to τ _{c} = 0.5 Pa and E_{B} = 10^{−8} J, whereas for deposited new snow, values of τ _{c} = 0.05 Pa and E _{B} = 10^{−9} J were chosen. These estimations of the critical shear stress are based on observations in the field during drift episodes. A description of these observations and a complete description of the grid and the boundary settings is given in Reference GauerGauer (1999).
As an example of the windfield measurements, Figure 9 depicts a temporal plot of the wind speed, M, and the direction, dd, approximately 2 m above the snowpack, for drift period 19–21 March at M76 and M73/74. The results of the corresponding numerical simulations are also shown. Generally, the simulations agree quite well with the measurements. A conspicuous feature is the turn in the wind direction from northwest on the windward side to more or less northeast to south on the leeward side of the crest. This phenomenon might be caused by a channeling effect of the topography (Fig. 8) and is reproduced also by the simulation (see also Fig. 13, shown later). While the simulated wind speeds tend to be somewhat too high windward of the crest, they are in good agreement in the lee. The scattering in the measured wind direction at mast M73/74 indicates the highly turbulent character of the wind field just leeward of the crest. It is still impractical to fully resolve this turbulent flow in the natural environment. Nevertheless, for of most of the time the simulation also reflects the wind direction at that mast quite well. Considering that windflow modeling for the severe complex topography of Alpine environments is still a topic of research, the simulated wind field gives reasonable results. Comparisons at the other masts and different heights, not shown here, confirm this.
So far, there is no suitable method available for measuring the snow depth over an extended area with reasonable effort and good accuracy. In our case, the required accuracy is better than ±0.1 m, which is of the same order as the new snow depth during a typical drift period in winter 1996/97. Hence, the simulation and field measurements were compared for a small area (about 200 × 40 m^{2}) around the crest.
Figure 10a shows the crosssection of profile line 3, which connects masts M76, M75 and M74/73 as indicated. Depicted is the topography used in the grid and, for comparison, the topography given by the digital map (DHM). Although the grid represents the crest quite well, it is obviously smoother than the DHM and even the real topography. Figure 10b compares the measured and simulated newsnow depth for the precipitation period 19–21 March. One can see a high variability in the measurements, partly influenced by the topography, partly due to small moving dunes. The simulation shows reasonable agreement with the measurements. The largest deviation was found on the windward slope just behind the small depression. A possible explanation for this deviation is given below.
An example of a spatial comparison is shown in Figure 11. To reduce the variability due to varying snowpack properties and the turbulent character of the wind, present in the individual measurements, and to highlight the more characteristic deposition features, an average over all six measured storm episodes was calculated. For the simulation, a “model” period was chosen. The comparison between both appears to be justified since it is claimed that the simulated model episode is prototypical of the storm episodes of winter 1996/97.
The measurements still show some variability, but also some trends can be discerned: one can see the erosion area in the windward slope near the crest line, the formation of a cornice at the ridge, and an area just leeward of the crest where the snow depth is less than the reference depth. Downslope areas where HN is twice HN_{ref} alternate with depletion areas such as the small terrace on the right. It should be recognized that the smallscale terrain features can cause large differences in the erosion and deposition patterns. The simulation does not show the same variability, but nevertheless reproduces some characteristic features such as the windward erosion area close to the crest, alternation between depleted areas and areas of enhanced deposition. As well as the uncertainties due to the poorly known snowpack properties (e.g. the critical shear stress), another reason for the discrepancy between the measurements and the simulation can be found in the history of the snowpack evolution. Thus, snowdrift tends to fill in hollows and to even out the topography in the course of a winter until an equilibrium level of the snow surface is reached where erosion balances deposition. This phenomenon could explain the discrepancy in the depression throughout the windward slope, particularly if the initially measured snow depth at this location was already close to its equilibrium level, whereas for the simulation a uniform initial snow depth was assumed. Thus the simulation might tend to deposit more snow as expected until the equilibrium level is reached. Nevertheless, the results of the simulation give reasonable images.
Another kind of comparison is shown in Figure 12b, which depicts the measured temporal evolution of the new snow at points M77 and M74/73 (windward and leeward of the crest; both masts were equipped with ultrasonic snowdepth sensors) during a drift episode, and the corresponding simulation. This comparison shows that, despite the rough estimates of the snowpack properties, the model can reproduce the temporal evolution of the snow depth quite well. Figure 12a depicts the corresponding wind speed at the top of the crest, and the precipitation rate.
Finally, Figure 13 depicts an example of the simulated redistribution pattern of the newsnow layer for the whole simulated area. The new snow is picked up in acceleration regions such as the area close to the crest line, and at small humps and brows. At the crest line, erosion of the old snowpack is observed as well. Deposition occurs in the deceleration regions leeward of the saddle, in small gullies and in hollows. In this simulation, depositions twice the reference height and erosion depth of about 0.1 m of the old snowpack are found at the crest line. The simulation gives a reasonable image of the redistribution pattern, which tallies with the observations in the field (see, e.g., Fig. 14).
4. Concluding Remarks
The present model is a step toward accurate numerical simulation of blowing and drifting snow in Alpine terrain. The primary emphasis was on a physically based description of snowdrift suitable for complex Alpine terrain. The model includes: a fully threedimensional nonsteadystate modeling; the modeling of the two main transport modes (saltation and suspension); the dynamical modeling of deposition and erosion, distinguishing between aerodynamic entrainment and ejecta due to particle impacts; the backreaction of the particles in saltation on the flow (twoway coupling); and the adaptation of the grid to changing snow heights. The comparison between simulations and field measurements is encouraging and indicates that numerical simulations of snowdrift may develop into a powerful tool for landuse planning and avalanche forecasting.
On the other hand, the first numerical simulations also demonstrate the problems with models involving complex terrain. For example, the grid resolution must match the length scale of the terrain. In an early simulation, the simulated wind field differed significantly from the measurements, and the 25 m resolution of the available digital terrain model of the Gaudergrat experimental site proved to be insufficient to resolve key topographic features, like the sharp crest line, in a reasonable way. Special effort had been undertaken to improve the terrain model and thus the numerical grid in the region of interest close to the crest line. A critical point in windtransport simulation is the nonlinear increase of the computational effort with increase of the area of interest and/or spatial resolution of the area. At that point, a compromise must be made between a high spatial resolution and an acceptable computational effort.
Poorly known boundary conditions are a common problem for windfield simulations in Alpine terrain. To overcome this problem, one might consider using the output of a mesoscale weather model. Some of these models are under development, but are not yet fully adapted for use in mountain areas.
How mechanical properties of the snowpack affect snowdrift rates is still unknown. It is reasonable to expect parameters like erodibility, the rebound probability of particles, the energy loss during impact, the rebound angle, etc., to differ for different snow types, but no systematic investigation has been carried out to study these relationships as yet. The present parameterization of erosion and deposition used in the model is mainly based on windtunnel studies and numerical simulations for aeolian particle transport. However, it is possible to include the dependency on varying mechanical snowpack properties into the model. For that, some research will be necessary to determine these mechanical properties for the varying conditions occurring in nature. This research should include coldroom windtunnel studies, similar to those done for aeolian particle transport, using different snow types.
Another still poorly known process is the mutual influence between particle transport and turbulence structure of the airflow within the boundary layer (surface layer) (e.g. the modification of saltation, transition to suspension, or damping of the turbulence). To this end, it would be of interest to do instantaneous profile measurements of the mass flux and the wind speed close to the surface with high temporal resolutions. With improving sensor technology this becomes more and more possible.
Acknowledgements
This research was funded by the Swiss Federal Office of Education and Science under the contract OFES 93.0082 as part of the European “Human Capital and Mobility” project “Contribution to avalanche dynamics with the aim of mapping hazardous areas”, contract CE CHRXCT 93–0307. I thank D. Issler and P. Föhn for many enjoyable discussions, and B. Brown for his comments. I am also indebted to W. Ammann, C. Fierz, M. Hiller, G. Klausegger, G. Krüsi, R. Meister and P. Weilenmann, as well as the staff of the mechanical workshop, for help in preparing the field measurements. Similarly, I am grateful to C. Camponovo, S. Harvey, C. Marty, R. Mullins, M. Phillips and K. Plattner for their assistance in the field, to A. Stoffel for his geographic information system support and to B. Gauderon for his support on computerrelated questions. I also thank the editors, two unknown reviewers and L. McKittrick for their comments.