Turbulence Model

Based on the universal conservation laws of mass and momentum, the governing equations for a Newtonian fluid are commonly referred to as the Navier-Stokes equations. Reynolds, or time averaging, introduces additional terms which must be modeled to obtain a closed system of equations. Restricting our interest to steady, two-dimensional flows, and adopting the concept of a turbulent-viscosity leads to the following continuity, and x and z-momentum equations:

(1)

(2)

(3)

where u and w are the mean velocity components in the x (horizontal) and z (vertical) directions, respectively; p is pressure; ρ is the fluid density; and ν_t is the turbulent-viscosity coefficient which is a property of the flow and may vary greatly over the flow field.

To close this system of equations, ν_t must be formulated in terms of the mean flow variables. Existing turbulence models range from simple models based on Prandtl’s mixing-length hypothesis, to those that use differential transport equations for the individual turbulent stresses and fluxes. The simple models based on the Prandtl mixing-length hypothesis suffer from several drawbacks. One primary failing is that the hypothesis implies that the turbulence is in local equilibrium, consequently the model is unable to account for the transport and history effects of turbulence. In addition, since it is difficult to prescribe the mixing-length distribution in any but the simplest flows, these models are not applicable to more complex, separating flows (Reference Rodi and KollmannRodi, 1980). A promising alternative is the k-∊ turbulence model of Reference Jones and LaunderJones and Launder (1972). This is a two-equation model which uses partial differential transport equations to describe the distribution of turbulent kinetic energy, k, and the dissipation rate, ∊, of the turbulent kinetic energy.

Since ν_t is proportional to k²/∊, it is assumed that

(4)

where c_µ is a constant. For the flow of interest here, the turbulent kinetic energy transport equation is

(5)

where term (i) describes convective transport of k, term (ii) is the diffusive transport, term (iii) describes the production of k through shear, and term (iv) is the rate of viscous dissipation. The ∊ equation is

(6)

where each term has a meaning similar to those found in the k equation. Based on the work of Reference Jones and LaunderJones and Launder (1972) the empirical constants in these equations are: c_μ = 0.09, c_1∊ = 1.44, c_2∊ = 1.92, σ_k = 1.0, σ_∊ = 1.3. For atmospheric flows, as opposed to engineering or hydraulic flows, the values of C_μ = 0.03 and C_1∊ = 1.16 have been suggested (Sutton and others, Reference Sutton, Brandt and White1986). Since the mid-1970s the k-∊ model has been tested extensively and has found success for a variety of different flows, including free shear flows, wall boundary layers, and flows involving recirculation.

Boundary Conditions

To solve the governing equations, boundary conditions must be applied to all boundaries of the computational domain for the dependent variables, u, w, p, k and ∊. Representative atmospheric profiles of velocity, turbulent kinetic energy, and dissipation rate are required for the inflow boundary conditions. Wind profiles within a neutrally stratified boundary layer can be approximated by

(7)

where u_* is the shear or friction velocity, z is the height above the surface, z₀ is the roughness length, and κ = 0.4 is von Karman’s constant. The rate of dissipation above the viscous sublayer is given by

(8)

(Reference SorbjanSorbjan, 1989). Reference Frost, Harper and FichtlFrost and others (1975) provided atmospheric profiles of turbulent kinetic energy k for the same conditions, where

(9)

Boundary conditions for the solid surfaces are formulated by choosing a computational domain which is slightly removed from the solid boundary, and then applying semi-empirical wall functions to the gap between the solid boundary and computational boundary. The logarithmic profile given by Equation (7) is used as the u velocity boundary condition on all horizontal solid surfaces including the top of the fence and on the tops of snowdrifts. In the presence of saltation, Reference OwenOwen (1964) suggested that the effective roughness length of the saltating surface is proportional to u_* ²/2g, where g is the gravitational acceleration. In light of the saltation modeling assumptions imposed in the current study, as discussed below, the implementation of a non-linear relationship between u and u_* does not appear justified. A complete summary of the boundary conditions applied to this problem is found in Table 1. Combining Equations (1)-(6) and the boundary conditions yields the coupled, highly non-linear system of equations which models the steady-state, two-dimensional turbulent flow field.

Table 1. Boundary conditions used in the computational model, where n if the direction perpendicular to the boundary

Saltation Model

The presence of saltation is dependent upon the relative magnitudes of the shear strength of the snow surface and the shear stress at the surface produced by the wind above. Experimental studies have found threshold shear-velocity values ranging from < 0.1ms⁻¹ for light, dry snow, to >0.4ms⁻¹ for old, hardened snow. Under moderate wind speeds, measurements above a horizontal surface show a very rapid decrease in mass flux with height through the saltation layer. Reference Kind, Gray and MaleKind (1981) reported that, for a wind speed of 10ms⁻¹, approximately 90% of the total saltating mass flux was contained within the first 3 cm above the surface. At higher wind speeds the suspension of particles in the turbulent-flow regime plays a significant role in the mass transport (Reference KindKind, 1992), and is strongly dependent upon the threshold shear velocity (Reference Pomeroy and GrayPomeroy and Gray, 1990).

Saltation models have been developed which describe the steady-state transport of saltating snow over uniform horizontal surfaces. For example, Reference Pomeroy and GrayPomeroy and Gray (1990) suggested

(10)

where Q_s is the total mass-transport rate per unit lateral dimension, u_*t is the threshold shear or friction velocity, and 0.68 is an empirical constant.

The saltation model implies that decreasing shear velocity leads to decreasing transport, and that when the shear velocity drops below the threshold velocity, mass transport will cease. Since Equation (10) was developed assuming a fetch long enough to provide a steady-state supply of particles, it is not directly applicable to the case of two-dimensional flow where the surface shear velocity can vary significantly over short distances. For example, consider the increase in shear velocity as the wind travels over a small bump: while the saltation model implies an increase in mass transport, that mass is not necessarily available. To apply the saltation model to the two-dimensional flow field of the present study, only the switching feature of the model is used, thus it only determines the presence or lack of transport. It is assumed that saltation ceases and snow accumulates at the location where the shear velocity drops below the threshold. In this simplified form, transport rates and the time required for drift traps to reach equilibrium cannot be computed.

The turbulent flow model is solved by discretizing the flow domain into an array of rectangular grid boxes. The surface shear velocity, u_*, is computed through the turbulent flow field calculation, and is then used to run the saltation model upon providing a threshold shear velocity, u_*t. Snow-accumulation patterns are developed by assuming that snow will be deposited within the surface grid box where the shear velocity first drops below the threshold. When enough snow accumulates to fill that grid box, it is blocked out and a new lower-boundary configuration results. The flow field is then recomputed and the process is repeated until an equilibrium drift has been formed which contains no regions of shear velocity below the threshold.

Flow Over a Solid Fence

Two experimental studies of atmospheric flow over a long, thin solid fence are used to verify the turbulent-flow model output. The first study is a full-scale field experiment (Reference JacobsJacobs, 1984), while the second study was conducted in an atmospheric wind tunnel (Reference PereraPerera, 1981). In Jacobs’ experiment, wind-flow profiles were measured at nine locations aligned perpendicular to a thin solid fence of 2 m height, 64 m length and 0.02 m thickness. The surface cover consisted of mowed heather with roughness length of 0.035 m. Wind speeds were measured using cup anemometers.

Duplicating this experiment with the computational model shows that all major features of the flow have been simulated. Jacobs found a recirculation eddy at the front of the fence, starting near the surface at about −0.5 h, where h is the fence height, and reattaching on the fence at a height between 0.5 h and 1.0 h. The model places the eddy at approximately −0.8 h for the surface and 0.6 h for the fence. Behind the fence Jacobs found a recirculation eddy which started at the top of the fence and extended a distance of between 5 h and 10 h to the ground surface. The model produced a reattachment point a distance 5.8 h from the fence. Figure 1 shows a comparison of the velocity fields measured by Jacobs and computed by the model. Displayed are the velocity profiles measured at locations x/h = −3.0, −1.0, 1.0, 3.0, 5.0 and 15.0, where the velocity has been scaled by the incoming velocity at fence height h. The model produces a realistic depiction of the measured velocity field. The largest discrepancies occur in the lee recirculation eddy where, although identified visually in Jacobs’ study, the cup anemometers have not resolved the flow direction. Further inaccuracies in field measurement techniques and variations in such features as surface roughness and wind-approach direction do not warrant a more detailed analysis of the velocity differences.

Fig. 1. Comparison of velocity profiles measured by Reference JacobsJacobs (1984) and computed by the model.

Perera’s experiment was also a two-dimensional, thin fence study, but here the roughness length was significantly smaller at 0.00036 m. Wind speeds were measured in a wind tunnel using pulsed-wire anemometers. These anemometers, unlike cup and hot-wire anemometers, have the ability to distinguish between positive and negative directions of wind flow. The computational reproduction of this study is visually similar to the Jacobs flow field: all major features of the flow have been reproduced. The smaller roughness length has extended the windward recirculation eddy to −1.1 h at the surface, and reduced the lee eddy to 5.5 h. Figure 2 compares the experimental and model velocity profiles at locations x/h = −1.25, 0.625, 1.25, 2.5, 5.0 and 15.0 for Perera’s study.

Fig. 2. Comparison of velocity profiles measured by Reference PereraPerera (1981) and computed by the model.

Snow Accumulation Patterns

The two-dimensional flow over a solid vertical fence is used as a test case for the coupled flow-saltation model. The fence is 2 m high and 0.5 m thick, and the incoming velocity profile is defined by Equation (7) with a 10 m velocity of 10ms⁻¹ and a roughness length of 0.001 m. The computed flow field for this boundary configuration is given in Figure 3.

Fig. 3. Computed flow field for test-case problem of single-phase, turbulent flow over a solid fence.

The equilibrium drift configuration computed for a threshold shear velocity u_*t = 0.2 ms⁻¹ is given in Figure 4. The evolution of the profile is illustrated in Figure 5, where the numbers indicate the sequence in which grid boxes were blocked out as the drift evolved. Figure 5 shows the drift first forming windward of the fence, in the region of the recirculation eddy, where a shear velocity below the threshold is first encountered. This drift builds at a position removed from the fence and then fills in between the fence and the drift. After the windward drift reaches equilibrium, snow begins to accumulate in the lee of the fence. The shoulder or scarp formation found in the lee drift at intermediate stages, for example the one present after the 41st blocked grid box, is a commonly observed feature of snowdrift formation (Reference MellorMellor, 1965; Reference TablerTabler, 1975). As the lee drift extends downwind of the fence, the flow field is modified and may produce additional accumulation at locations windward of the main drift scarp. This is illustrated, for example, by the blocking of the 51st grid box. This result also occurs with the drift windward of the fence. Increasing (or decreasing) the threshold shear velocity will increase (or decrease) the size of the equilibrium drift.

Fig. 4. Equilibrium snowdrift and flow-field configuration computed for a threshold shear velocity of u_*t = 0.2 ms⁻¹ for the test-case problem of single-phase, turbulent flow over a solid fence.

Fig. 5. Progression of snowdrift development for the case u_*t = 0.2ms⁻¹. The numbers indicate the sequence in which the model blocked grid boxes as the drift evolved towards the equilibrium configuration.

The lee equilibrium drift surface has a slope of 10–13%. For the case of embankment slopes descending from a level plain, Reference FinneyFinney (1939) found that no snow accumulates on embankment slopes of less than 17%. Reference Berg and CaineBerg and Caine (1975) and Reference TablerTabler (1975) found only small accumulations on slopes of 17% for this same 0% approach-slope configuration. The relatively steep approach slope described by the drift windward of the fence in the computational experiment produces a longer drift than those produced in the above examples which have approach slopes of zero. This is consistent with Tabler’s study where, for an approach slope of approximately 10%, snow accumulated on a lee slope of approximately 10%. Planned model improvements which include accounting for the transport of snow due to turbulent suspension and/or precipitation are expected to aid in filling down-wind regions and help create snow-accumulation sequences which more closely fit natural accumulation patterns.