Computational fluid dynamic (CFD) simulation of snowdrift in alpine environments, including a local weather model, for operational avalanche warning

Abstract

A new continuum approach to snowdrift modelling is introduced. In addition, numerical studies are carried out to identify the influence of time-varying wind conditions on snowdrift simulations. We compare the snowdrift patterns at Grimming mountain, Austria, derived using a time-averaged wind field and a time-varying wind field obtained from the numerical weather prediction model INCA. The results show significant differences in the deposition patterns and snow depth even after a 6 hour drift period. Using time-averaged boundary conditions leads to an underprediction of the resulting snow depth caused by averaging the wind speed, which lets gusts of wind disappear while snow transport is a non-linear function of the wind speed. Using numerical weather prediction models for snowdrift simulation therefore provides enhanced knowledge of the snow depth for local avalanche warning services.

Introduction

The forecasting of avalanche danger has become more important as a result of the increasing development of alpine environments for transport and tourism. The protection of civil facilities and human lives is a major aim of avalanche forecasting. Primarily, the influence of blowing and drifting snow on avalanche formation has to be emphasized. The snow depth on mountainsides and the dead load of snow slabs are important factors that determine avalanche occurrence. In addition to freshly fallen snow, snowdrift contributes significantly to the total snow on mountainsides. Mainly leeward, the snow depth is determined by the amount of drifting snow. Even on mountainsides with a slope angle as low as 30 ^˚ , very high snow depth significantly increases avalanche danger. Finally, drifting snow forms large snow cornices at mountain ridges. If those cornices break off they may trigger avalanches. Hence, wind-deposited snow forms can cause avalanche release.

It is very difficult to predict avalanches. Meteorologists have to rely on selective measurements of the actual weather situation and selective snow profiles. The avalanche hazard is determined from this sparse information (for a definition of what constitutes a hazard, see, e.g., Fell and others, 2005). A knowledge of zones where wind deposits snow and forms cornices as a function of the overall weather situation is of vital importance for avalanche warning.

The Physical Concept

Since Bagnold (1941) it is common to distinguish between three different transport modes of snow. Here, we use the more formalized definition of Anderson and others (1991) (Fig. 1).

Fig. 1. Transport modes of blowing and drifting snow: creeping/ reptation, saltation and suspension, indicating the different volume fractions in the different modes. Snowpack indicates the structural snow cover. Open circles designate grains bound on the snowpack and full circles designate grains in motion (from Schneiderbauer, 2006).

Surface creeping/reptation The transport modes that are located nearest to the ground, up to a height of the order of the surface roughness, are referred to as surface creeping and reptation. The transport of grains is not affected by the wind but depends on the impact of other particles. Nevertheless, this transport mode contributes weakly to snow transport. Typical grain sizes are of the order of 1 mm.

Saltation The grains are partially influenced by wind turbulence. If the trajectories of the grains are not influenced by turbulence the motion is pure saltation, otherwise modified saltation. Saltating grains are able to rebound and to eject other grains. The typical grain diameters in the saltation mode are 0.05–0.5mm. Most of the snow transport occurs in this mode, whereas the volume fraction of the snow particles is in the range 10 ^{− 5} to 10 ^{− 3}. Furthermore, for high wind speeds the height of the saltation layer can reach about 0.1 m.

Suspension In the suspension transport mode, located above the saltation layer, the snow particles are lifted far away from the surface and can reach a height of several hundred metres. This mode is referred to as long-range transport. The typical volume fractions are lower than 10 ^{− 4} and the grains are usually smaller than in the saltation layer.

Saltation starts at wind speeds u ₁₀ (wind speed measured at 10m height) of about 5–8ms ^{− 1}, whereas the different transport modes contribute to the snowdrift as listed in Table 1 (e.g. Gauer, 1999). As previously mentioned, saltation contributes most to the mass transport.

Table 1. Distribution of total mass transport; the values vary with wind speed

Snow transport can be described as follows. If the wind shear exceeds a certain threshold, the so-called fluid threshold (Bagnold, 1941), grains will be entrained and set in motion. They are easily accelerated by the wind because of their small mass and diameter. Already entrained grains contribute to the total shear, i.e. the wind shear required to entrain grains is lower in the presence of particles due to particle impacts. In addition, the heights of the transport modes scale with the wind speed. If wind speed is above the fluid threshold, a higher amount of snow can be transported. Due to the interaction between grains and wind, the wind field is modified. However, if the total shear is below a second threshold, the impact threshold, snow particles will be accumulated. Grains the motion of which is directed towards the snowpack are deposited or rebound with a certain probability. Only a certain number of grains can be transported by the wind. This leads to deposition when saltation exceeds saturation.

The transport processes discussed above lead to an evolution of the snowpack. The snow cover is growing in deposition zones, and the new shape therefore influences the local velocity field. The snow depth is reduced in erosion zones due to particle entrainment. Hence, the modification of the wind field causes time-dependent spreads of these zones.

State Of The Art

Two different approaches can be taken to numerically simulate avalanche hazards.

1. Simulation of the dimensions and destruction of avalanches for given boundary conditions such as snow depth, temperature, sun radiation, humidity, etc.

2. Snowdrift simulations or other snowpack simulations to determine the initial and boundary conditions for avalanche activities.

Recapitulating, the simulation of avalanche hazard zoning is focused on protecting threatened zones from avalanches; not on avoiding avalanches. This work is focused on snowdrift processes, and the following overview of existing work concentrates on snowdrift modelling.

Snowdrift models

In the following, previous work on snowdrifts is outlined in chronological order. Bagnold (1941) laid the foundation for our current understanding of aeolian snow and sand transport. In his wind-tunnel and field measurements he focused on the saltation modes, assuming that these contribute most to the mass transport. In addition, he derived an empirical relation between the fluid-induced shear stress and a logarithmic wind profile (Prandtl, 1965):

(1)

where u_* is the friction or drag velocity, z ₀ the surface roughness and κ the von Kármán constant (= 0.4). The quantity u_* is a mathematical abbreviation for the expression u_* = where τ is the fluid-induced shear stress and ρ _f denotes the density of the fluid. From measurements, Prandtl (1965) obtained some erosion and deposition criteria which

depend on a critical shear stress. For the erosion and accumulation processes, he found

(2)

where u_{* t,e} is the threshold friction velocity for the erosion process, u_{* t,i} is the threshold friction velocity for the impact process, A _e,i denote empirical parameters, ρ _p is the particle density, d _p is particle diameter and g is the acceleration due to gravity.

Anderson and others (1991) gave a more formalized classification of snow and sand transport. A more accurate description can be found in their section on the physical concept. Anderson and Haff (1991) estimated a linear dependence between the number of entrained grains per unit time and per unit area and the shear stress induced by the fluid. Additional work, based on work by Bagnold (1941) and Owen (1964) on the saltation transport mode, was done by Pomeroy and Male (1992), who derived semi-empirical functions of the height of the saltation layer, the mass flux in the saltation mode and the saturation concentration in saltation. For example, the height of the saltation layer is given by

As an early application of computational fluid dynamics, Waechter and others (1997) compared the velocity field of a flow around an arctic building with observed snow depositions and estimated that, in regions of low velocities, snow deposition is possible. Nevertheless, the work did not include important physical concepts such as modelling the snow phase or the evolution of the snowpack due to erosion and deposition processes.

Liston and others (1998) presented a three-dimensional (3-D) snow transport model to simulate the snow-depth evolution over topographically variable terrain. It is also important to note that the influence of vegetation on drifting snow has to be emphasized.

Based on the work of others (Bagnold, 1941; Clift and others, 1978; Anderson and Haff, 1991; Anderson and others, 1991), Gauer (1999) introduced a 3-D non-steady-state model in which the turbulent suspension is modelled by an incompressible fluid (air) and determined by computational fluid dynamics. The additional mass transport in the saltation mode is obtained by typical grain trajectories. Furthermore, ejection of grains due to impacts with other snow particles is included.

Two important simplifying approximations were made: (1) the modified saltation is not considered, and (2) the wind profile is assumed to be logarithmic within the saltation layer.

Thus, the formation of snow cornices on mountain ridges is not incorporated within such a modelling approach. This two-way coupling model includes the feedback caused by the evolution of the snowpack due to erosion and deposition processes. Precipitation is also included in the model. The numerical simulation results were reinforced by the related field measurements. An important restriction of the model is the high computational demand. An improved model, in which a simplified deformation of the wind field within the saltation layer is considered, was developed by Doorschot (2002).

Lehning and others (2000) combined the snow-cover model SNOWPACK (e.g. Lehning and others, 1999) with a snowdrift simulation, simulating a stationary wind field profile for a chosen time interval (about 1 hour) to avoid unacceptable computational expense.

Bagnold’s formulae (Equation (2)) were verified by Nishimura and Hunt (2000) inwind-tunnelmeasurements.Naaim-Bouvet and others (2004) compared field and wind-tunnel measurements with numerical simulations. On the basis of Pomeroy and Male’s (1992) semi-empirical relationships for the saltation layer, they introduced a simple two-dimensional (2-D) numerical model. Furthermore, they investigated several snow fences with different porosity to determine their influence.

The development of desert dunes or snow cornices is also an important consequence of sand or snowdrift. Many villages in desert regions are threatened by dunes. Sauermann (2001) introduced a model, fundamentally based on Bagnold (1941), to simulate the development of such formations.

In the models discussed above, snowdrift is considered as a continuum. The mass fluxes are determined by typical particle trajectories within the saltation layer (Gauer, 1999; Lehning and others, 2000; Nishimura and Hunt, 2000; Doorschot, 2002). To depict the saltation transport mode within a Lagrangian frame, one would need 10⁶ particles m ^{− 2} (Gauer, 1999),which would lead to an unacceptable computational effort in a non-steady-state simulation. Trenker and others (2005) investigated a Lagrangian particle model in a stationary wind field. The particles followed the wind field, because ’two-way coupling’ was neglected. Bagnold’s criteria were used to identify deposition zones.

To summarize, one can perceive the following:

• 1. The continuum saltation models are based on the properties of typical grain trajectories.

• 2. The migration of snow cornices on mountain ridges is not modelled satisfactorily.

• 3. The modelling of every particle trajectory in the saltation layer is nearly impossible (10⁶ particles m ^{− 2}).

• 4. It is necessary to simulate in three dimensions to obtain relevant and reliable predictions in practice.

• 5. A sophisticated effort is required using 3-D models.

• 6. The underlying snowdrift mechanism is not fully understood.

Schneiderbauer (2006) introduced a fully 3-D, non-steady-state snowdrift model with several improvements compared with the existing models. The determination of the accumulation and erosion zones is managed by Bagnold’s criteria (Bagnold, 1941), but, contrary to the work discussed above, the whole snow transport is mapped by a continuum approach (Schneiderbauer, 2006). Key features of this snowdrift model include the following:

• 1. Snowdrift patterns are predicted, on the basis of Bagnold’s erosion and deposition criteria. The erosion and accumulation criteria by Bagnold are supplemented by gravityinduced shear stresses depending on the slope angle.

• 2. The snow cover is modelled dynamically, coupling changes in geometry to the flow. Deformation of the snow cover is mapped in the node-normal direction (Schneiderbauer, 2006, p. 66).

• 3. The model is fully turbulent. The whole snow transport and the computation of the shear stresses is affected by turbulence.

• 4. The mass transport is assigned by the balance equations of multiphase flows.

Algebraic Slip Model

The algebraic slip model is based on the balance equations for interpenetrating phases, where the phases are allowed to move with different velocities (Fluent Inc., 2005). Hence, we account for the snow particles in saltation and suspension as a continuous phase. The saltation and the suspension layers are not separated numerically (as in Gauer, 1999, or Doorschot, 2002). Physically, the difference between those transport modes is due to different grain sizes. Saltating grains are too large to be lifted by the flow. Thus, snow transport is provided by the influences of the flow, averaged particle properties (e.g. size, shape, etc.) and gravity. Bagnold’s impact and erosion criteria (Bagnold, 1941) distinguish between zones of erosion and zones of deposition. The great advantages of the continuum approach are that the model is fully turbulent (so the whole snow transport is affected by turbulence effects) and that particle interaction is modelled by the algebraic slip model.

Finally, saturation effects within snow transport are considered in the model. Empirical relationships (Naaim-Bouvet and others, 2004), such as the height of the saltation layer which influences the saturation volume fractions in the finite control volumes of the computational domain (Schneiderbauer, 2006), are of vital importance for snowdrift computation.

Identification of zones of accumulation and erosion

Bagnold (1941) derived an impact and aerodynamic entrainment threshold friction velocity (u_{* t)i,e} :

(3)

where ρ _p and ρ _a are the densities of the snow phase and air, respectively, d _p is the particle diameter and A _i,e are dimensionless empirical parameters for which the subscripts i and e denote the impact and erosion process, respectively. Equation (3) was determined from measurements of a logarithmic wind profile at height z as in Equation (1), which assumes flat terrain. The friction velocity u_* is a mathematical symbol for the shear stress of a logarithmic wind profile (cf. Equation (3)):

(4)

where τ _a is the shear stress applied by a fully developed turbulent wind over a flat surface onto the snowpack (Andreotti and others, 2002). Therefore, from Equation (3) we obtain the shear stress thresholds τ _ci,e for the impact and erosion process:

(5)

These thresholds are directly proportional to the gravity force F _g acting on the particle with diameter d _p and cross-section A _p:

(6)

Criteria (5) and (6) are valid for typical sizes of saltating snow grains (Schneiderbauer, 2006, p. 55).

For a plane with slope angle γ and more complex terrain, the logarithmic wind profile (Equation (1)) is modified, and the air-induced shear stress τ _a is determined by

(7)

where μ denotes the viscosity of air and ∂u/∂n the velocity gradient in the surface normal direction. For particle-laden flows, impacting particles induce an additional contribution to the total shear acting on the surface,

(8)

where u _m is the velocity of the air–snow mixture and μ _m is its viscosity, defined as

(9)

where α denotes the volume fraction of the snow phase. In our approach, the snow phase is considered as a very viscous fluid with viscosity μ _p. An additional shear stress τ _g induced by the gravity force acting on the particles at the surface of the snowpack appears in the force balance. Hence, the total shear stress acting on the surface is

(10)

We assume that this additional shear, as illustrated in Figure 2, is a function of the density of the snowpack ρ _s, the particle volume V _p, the particle cross-section A _p and slope angle γ:

(11)

Fig. 2. Illustration of the different forces acting on a snow particle as well as the accumulation and erosion processes (from Schneiderbauer, 2006).

where n(γ) denotes the surface unit normal vector as a function of the slope angle γ. Thus, the gravity-induced shear Equation (11) can be considered as the grade resistance scaled with A _i,e similar to the shear stress threshold Equations (5) and (6). The gravity-induced modification of the shear stress can also be considered as a shift of the shear stress thresholds. Hence, for a spherical particle with diameter d _p and cross-section the identification of accumulation zones is

(12)

and the identification of erosion zones is

(13)

Here, τ _i is the modified air-induced shear stress for the accumulation process and τ _e is the modified air-induced shear stress for the erosion process.

Accumulation of drifting snow

The amount of accumulated snow per unit time ṁ _s, as shown in Figure 2, can be expressed as

(14)

where S denotes the surface of the snowpack, α the volume fraction of drifting snow and u _p the velocity of the snow particles within the continuous snow phase. The relative velocity u _ap between the snow phase and air phase is given by (Fluent Inc., 2005)

(15)

with

where t _k denotes the particle relaxation time, f _drag is the drag function and a is the snow particle’s acceleration. According to Schiller and Naumann (1933), the drag function is

(16)

where Re is the dimensionless particle Reynolds number. A more detailed description is given in the Fluent 6.2 Users’ Guide (Fluent Inc., 2005, ch. 23.3.4).

The mass-balance equation holds for the accumulation process:

(17)

where volume V and surface S in Equation (17) are indicated in Figure 2. The change of the volume fraction of the snow phase within the calculation of ṁ _s can be expressed as

(18)

The rate of change of the volume of the snowpack V̇_s can be determined from Equation (18) if we divide the mass flux ṁ _s by the density of the snowpack ρ _s:

(19)

The maximum packing limit α _pl in Equation (19) is given by the ratio ρ _s /ρ _p.

By considering the impact threshold Equation (5) and the modification of the threshold due to gravity Equation (12), the mass flux to the snowpack can be expressed as a function of τ _i and τ _c,i , i.e.

(20)

where θ(x) denotes the Heaviside step function, defined as

Rebounding of the snow grains is neglected here, because of the small particle velocities below the impact threshold.

Snow-entrainment process

The exact underlying mechanism responsible for the initiation of the snowdrift process is not completely known. Anderson and Haff (1991) estimated that the number of entrained grains per unit area (denoted by N _e) and unit time depends linearly on the excess shear stress:

(21)

where τ _c,e denotes the critical shear stress from Equation (5) and ξ is an empirical constant with dimensions of (force×time) ^{− 1}.

Using Equations (10) and (13),we can calculate the amount of entrained snow per unit time from Equation (21), i.e.

(22)

where ρ _p is the particle density and V _p the particle volume. From the calculated mass flux we can determine the change of volume of the snowpack in the same way as in Equation (19):

(23)

The change of volume fraction of the snow phase in the control volume V is accounted for in the same way as in Equation (18):

(24)

Saturation volume fraction

Pomeroy and Male (1992) estimated the concentration of snow in the saltation layer c _salt as a function of the air-induced shear stress,

(25)

Hence, no saltation will be enabled if the shear stress threshold τ _c is not exceeded. However, one has to keep in mind that the airflow can only carry a certain amount of snow particles (Gauer, 1999) which is of the order:

(26)

(kgm ^{− 3}). If the volume fraction assigned by the flow is higher than the saturation, an additional mass flux to the snowpack occurs:

where α̇ _s can be approximated as

(27)

Height of the saltation layer

Owen (1964) proposed that the height of the saltation layer h _salt scales with the square of the friction velocity u_* :

(28)

where β _salt is a dimensionless scale factor. In addition, for the surface roughness length,

(29)

Naaim-Bouvet and others (2004) determined β _salt to be equal to 1.6.

In our model, u_* is calculated as a function of the spatial coordinate x and is given by Equation (4), i.e.

(30)

The height of the saltation layer and the saturation volume fraction give an upper limit for the concentration of the snow phase in the finite control volumes (e.g. Schneiderbauer, 2006).

Turbulence modification of accumulation process

One major problem affecting the simulation of snowdrift occurrences is the discretization of the calculation domain. All flow variables are averaged values over finite control volume V:

In the PARTRACE^TM simulation (Beimgraben, 1997), which can calculate the sedimentation of dissolved dust particles, the following expression for the turbulent kinetic energy is applied:

(31)

where u^' _i represents the fluctuating parts of the velocity components u_i (Fluent Inc., 2005). For a given turbulent kinetic energy and for a uniformly distributed fluctuating velocity u^' , an additional turbulent velocity in the direction of the unit surface normal is given by

(32)

Hence, the snow velocity u _p used for the mass flux calculations is corrected as follows:

(33)

where n denotes the outward unit surface normal vector.

Numerical Model

In our algebraic slip model approach, the saltation and suspension layers are not treated separately in different calculation domains as by Gauer (1999) and Doorschot (2002). (For details of the algebraic slip model in FLUENT see Fluent Inc., 2005.) Rather, the saltation layer is modelled by volume fractions of the snow phase in the volume cells adjacent to the snow cover. The mass fluxes are obtained from the flow field of the snow phase. Bagnold’s ‘stick–slip’ criteria (Bagnold, 1941) are used to distinguish zones of deposition and aerodynamic entrainment. A saturation volume fraction was introduced to avoid non-physical effects in the flow field. The deformations of the snow cover are predefined by the mass fluxes and the unit surface normals. A dynamic mesh model will remesh the domain if yield criteria are exceeded. In addition, precipitation can be considered.

The numerical modelling of snowdrift using an algebraic slip model approach can be divided into two separate calculations for each time-step: (1) calculation of the flow field by the commercial computational fluid dynamics (CFD) software FLUENT and (2) computation of the mass fluxes between snowpack and saltation transport mode in order to use Bagnold’s erosion and deposition criteria, respectively.

The flowchart of our snowdrift model is shown in Figure 3.

Fig. 3. Flowchart of the numerical snowdrift model (from Schneiderbauer, 2006).

Computation of the flow field

The flow field is calculated using the algebraic slip multi-phasemodel provided by FLUENT. Hence, the snow particles are characterized by a continuum, and the relative (slip) velocities between air and snow are given by algebraic Equation (15).We use three different phases to simulate snowdrift and precipitation:

• 1. air (primary phase);

• 2. precipitation of snow (secondary phase, not considered within this work); and

• 3. snowdrift (secondary phase).

In addition, physical boundary conditions such as velocity profile provided by Integrated Nowcasting through Comprehensive Analysis (INCA) are applied to obtain a realistic flow-field computation.

Mapping of deposition patterns

The grid update for a node x_j on the boundary between fluid and snowpack is given by the average displacements of the adjacent faces l. These displacements h_l for a face l are computed from the mass fluxes as

(34)

where the subscripts i and e denote an impact and an erosion process, respectively. ΔV_l indicates the change in the volume of the snowpack that follows from Equations (19) and (23). If there are k faces in the neighbourhood of node xj, the update is given by

(35)

where n_m denotes the unit surface normal of face m and is defined to point out of the domain. In general, the motion of a node x_j is not very fast compared with the velocities of the flow due to very low volume fractions. Hence, we can decouple the timescales of the grid update and of the flow-field computation.

Decoupling of the timescales

For future applications of this fully developed 3-D snowdrift model (e.g. operational avalanche warning systems) saving computational time is an important aspect to obtaining simulations in real time. Furthermore, solving the non-linear Navier–Stokes equation contributes predominantly to the computational effort in our snowdrift model.

We introduce a factor S, which decouples the timescale T_flow of the flow field and the timescale of the update process _gridupdate:

(36)

Therefore, the motion of a node x_j is accelerated by a factor . If we consider _gridupdate as the physical timescale, the flow field will be kept constant for a time period of _gridupdate. Note that the decoupling will be applicable only if

is satisfied. The set N contains all nodes on the snow cover.

To find an upper limit for the scaling factor , Schneiderbauer (2006) performed a parameter study of a simple 2-D test case. It was found that for < 100 the snowdrift simulation is not affected significantly.

Local Weather Model INCA

INCA is being developed at the Zentralanstalt für Meteorologie und Geodynamik, Vienna, Austria, and is used to obtain time-varying boundary conditions. A fine resolution of 1 km×1 km is made possible by using several data sources during the calculation of the different meteorological fields (e.g. wind, pressure, etc.) as shown in Figure 4 (Haiden and others, 2007). In this work, we take advantage of the Nowcasting mode of INCA (there are two other modes, Analysis and Short-range; Haiden, 2005) to provide an (up to) 6 hours forecast of the wind speeds and wind direction for the simulation of snowdrift in alpine environments. Furthermore, we concentrate on the appropriate integration of wind fields, neglecting other meteorological data.

Fig. 4. INCA uses several different sources for the calculation of meteorological fields; the single most important data source is a network of about 150 automated surface stations called Tawes. The numerical weather prediction model Aladin is used as a first guess for the 3-D analysis of temperature, humidity and wind.

INCA is implemented by overlaying the INCA dataset onto the topographic dataset and carrying out a bidirectional interpolation onto the boundary using

(37)

with

(38)

where u _B denotes the west-wind component to be calculated on the boundary (the same formula applies for the calculation of the south-wind component v _B), u_ij are the corresponding west-wind data of the INCA dataset, x and y are the Cartesian coordinates of the boundary element in the actual scope and [. . .] indicates the truncation operator.

Results

In order to justify the effort needed to implement INCA, we compared computed snowdrift patterns at Grimming mountain, Austria, to show the difference between real time-varying wind fields and time-averaged wind fields. The time averaging of the 6 hours INCA dataset of wind speeds leads to u _B = 11.8ms ^{− 1} and v _B = 0ms ^{− 1}.

At the bottom, an erodible layer of cohesionless snow particles with a thickness of 2 m was assumed. The typical particle diameter was set to 300 μm, and the density of the particles was set to 900 kgm ^{− 3}. The parameters A _i and A _e in Equation (2) were set to 0.14 and 0.16, respectively. Hence, the critical wall shear stresses (Pa) are given by

For this simulation, a k-ϵ turbulence model was used and we assumed the turbulence parameters of turbulent intensity I and turbulent viscosity μ _t at the boundaries:

where μ is the viscosity of air. Finally, to save computational time, the scale factor S was set to 50. The required computational times (including the wind-field computation), using a Dell Precision 300 personal computer with an Intel® ^R Core^TM 2 2.4GHz processor, were (1) for time-averaged boundary conditions, 220 min and (2) for time-varying boundary conditions, 622 min.

The dimensions of the simulation area were 10.6 km in the east–west direction and 8.8 km in the north–south direction. The upper boundary was located 2 km above the highest point of the mountain. The horizontal resolution within the area of interest is 15 m, with a linear increase of 10% towards the boundaries. In Figure 5, the horizontal distribution of the grid is shown. The vertical thickness of the first gridcells is 1 m with a vertical stretching of 10%. Finally, the whole grid contains 296 450 hexahedral cells.

Fig. 5. Horizontal distribution of the gridcells used in the simulations. The arrow points northwards and the colours indicate the elevation in m. The generation of the grid takes about 5 s using a 25 m resolved digital elevation model (DEM) dataset.

The difference in snow heights caused by the use of real meteorological data is significant. Less snow is eroded using averaged wind fields because of lower local wind speeds (due to averaging, wind-speed peaks disappear) (Figs 6 and 7). At most, 5% of the initial snow depth is eroded and therefore deposited. For example, the INCA dataset contains wind speeds up to 25 ms− ¹ whereas the averaged boundary conditions are given for a wind speed of 11.8 ms− ¹. In comparison, the use of the INCA dataset leads to more snowdrift where snow is displaced from mountain ridges towards mountain chutes (Figs 8 and 9). At most, 15% of the initial snow depth is eroded and the maximum increase of the snow heights is also about 15%. Figures 6–9 highlight the difference between using an averaged wind field or a real INCA dataset. For Figures 6–9, the colour bar corresponds to the snow depth h/h _ref, where h _ref = 2 m.

Fig. 6. Snowdrift pattern at the northeast slope of Grimming mountain, after 6 hours, based on the time-averaged wind field. The arrow indicates the mean wind direction (westerly).

Fig. 7. Snowdrift pattern at the southeast slope of Grimming mountain, after 6 hours, based on the time-averaged wind field. The arrow indicates the mean wind direction (westerly).

Fig. 8. Snowdrift pattern at the northeast slope of Grimming mountain, after 6 hours, based on the time-varying wind field.

Fig. 9. Snowdrift pattern at the southeast slope of Grimming mountain, after 6 hours, based on the time-varying wind field.

As well as the higher accuracy, the time needed to obtain a solution increases about 280% due to lower convergence of the fluid dynamic solver because of varying boundary conditions.

Conclusions

Figures 6 and 8 show the differences when computing the snow heights using an averaged wind field and a varying wind field. These differences cannot be neglected if we wish to obtain realistic snowdrift patterns in alpine environments. The much-improved result achieved by using numerical weather prediction models such as INCA justify the huge increase of computational time.

The simulation of snowdrift processes by our method seems to be a good tool to support operational avalanche warning services. For future applications, a tuning of the model parameters will be necessary to obtain quantitative and qualitative results of snow heights. Furthermore, additional research has to be carried out to find appropriate boundary conditions to include precipitation events in the simulations. Until now, temperature has been neglected in the snowdrift model. However, temperature affects the accumulation and erosion criteria significantly (e.g. cohesion force between the snow grains). In the present state, the operational use is limited by computational times. The area considered is therefore limited by spatial size and resolution. In order to save computational time, the implementation of the snowdrift model has to be parallelized to run on host systems. Finally, the quality of the NWP model data influences the accuracy of the snowdrift simulations.

The presented model provides a qualitative measure of a potential avalanche area in a mountain chute on the south-east slope (Fig. 9). This area is under observation by radar measurements by the avalanche alert service of Styria, Austria, to detect avalanche releases for road closing. However, the computation using an averaged wind field does not detect any snow deposition there (cf. Fig. 7). It remains of vital importance to consider time-varying weather data for the prediction of snow depositions for operational avalanche warning.