A snow-transport model for complex terrain

Abstract

As part of the winter environment in middle- and high-latitude regions, the interactions between wind, vegetation, topography and snowfall produce snow covers of non-uniform depth and snow water-equivalent distribution. A physically based numerical snow-transport model (SnowTran-3D) is developed and used to simulate this three-dimensional snow-depth evolution over topographically variable terrain. The mass-transport model includes processes related to vegetation snow-holding capacity, topographic modification of wind speeds, snow-cover shear strength, wind-induced surface-shear stress, snow transport resulting from saltation and suspension, snow accumulation and erosion, and sublimation of the blowing and drifting snow. The model simulates the cold-season evolution of snow-depth distribution when forced with inputs of vegetation type and topography, and atmospheric foreings of air temperature, humidity, wind speed and direction, and precipitation. Model outputs include the spatial and temporal evolution of snow depth resulting from variations in precipitation, saltation and suspension transport, and sublimation. Using 4 years of snow-depth distribution observations from the foothills north of the Brooks Range in Arctic Alaska, the model is found to simulate closely the observed snow-depth distribution patterns and the interannual variability.

1. Introduction

In the Arctic, on the prairies and in high-alpine regions, the combination of wind and snowfall play an important role in determining the winter snow distribution. Snow is scoured from ridges and windward slopes, and deposited in the lee of bluffs and other topographic obstructions; this wind-driven transport creates snow depths that can vary by a factor of ten or more over short distances. The interplay between wind and topography is the primary driver of this lateral snow movement, but vegetation and the moving snow itself also affect the transport through complex feedbacks and interactions. There has been considerable effort applied towards understanding and modeling snow-transport processes, with the goal of simulating and predicting snow distributions under a given set of topographic and weather conditions. These interactions are conceptually simple but, when trying to model all of the interactions occurring under a wide range of naturally occurring meteorological conditions and topographic configurations, model formulation and implementation can be challenging.

There are numerous reasons for needing to simulate or predict the wind redistribution of snow. In alpine regions, wind loading can lead to dangerous avalanche conditions and predictive snow-transport models can indicate when to close threatened highways, railroads and ski runs. In Arctic tundra regions, the uneven distribution of snow exerts strong control over plant-community distribution (Evans and others. 1989). This plant distribution, in conjunction with the snow-depth distribution, determines where caribou, musk-oxen and other animals can find winter forage. Assessing sustainabiliiy of animal harvest for human consumption or predicting the impact of blizzards on the animals requires models capable of simulating the snow distribution. At the largest scale, seasonal snow cover, because of its low thermal conductivity, high albedo and high degree of spatial and temporal valiability, plays an important role in governing the global radiation balance. In the Northern Hemisphere, the mean monthly land area covered by snow ranges from 7 to 40% during an annual cycle, making snow the Earth’s most rapidly varying large-scale surface feature (Hall, 1988). In regimes where snow drifting creates non-uniform snow accumulations, melting in spring creates a patchy mosaic of vegetation and snow that is particularly complex to describe and model, and has a direct impact on energy exchanges between the land and atmosphere (Liston, 1995). This patchy snow distribution can play an important role in determining the timing and magnitude of snowmelt run-off. To capture these processes in climate and hydrologie models, accurate descriptions of snow distribution are necessary.

The snow-transport model described in this paper builds on many previous snow-transport studies. These studies show that snow saltation and suspension are the main mechanisms for moving snow (Mellor, 1965; Kobayashi, 1972; Radok, 1977; Male, 1980; Takeuchi, 1980; Kikuchi, 1981; Schmidt, 1982, 1986; Maeno and others, 1985; Pomeroy, 1988; Pomeroy and Gray, 1990), and the empirical and analytical relationships between wind speed and transport rates developed by these researchers are found in their snow-transport models and those developed by others (e.g. Tabler and others, 1990; Uematsu and others, 1991; Liston and others, 1993a; Bintanja, 1998a, b; Sundsbø, 1997; Purves and others, 1998). Other studies have shown that sublimation of wind-transported snow can be important, reducing the snow-cover water balance by a significant fraction under certain conditions (Goodison, 1981; Benson, 1982; Benson and Sturm, 1993; Pomeroy and others, 1993, 1997). The rate of sublimation, though complicated, has been found to be a function of wind speed, air temperature, humidity, particle size and solar radiation reaching the particle surface (Thorpe and Mason, 1966; Schmidt, 1972, 1982, 1991; Lee, 1975; Tabler, 1975a, 1987; Male, 1980; Pomeroy, 1988: Déry and Taylor, 1996). Studies have also shown that there are changes in aerodynamic roughness with increasing wind velocity, due to the drifting snow itself (Owen, 1964; Budd and others, 1966; Kobayashi, 1979; Schmidt, 1980, 1982; Kikuchi, 1981; Greeley and Iversen, 1985; Pomeroy, 1988; Déry and Taylor, 1996), and the resulting formulations are frequently used in snow-transport models.

The patterns of accumulation and erosion around variable topographic features have been described (Mellor, 1965; Tabler and Schmidt, 1986; Benson and Sturm, 1993), simulated in wind-tunnels or using small-scale models (Finney, 1939; Kind, 1976; Tabler and Jairell, 1980; Iversen, 1981), and modeled numerically. Two-dimensional, or cross-section, models have been developed (Berg and Caine, 1975; Tabler, 1975b; Uematsu and others, 1991; Liston and others, 1993a; Sundsbø, 1997) that can predict the evolution of snowdrifts and the distribution of snow-related features along the wind-flow direction; for example, both upwind and downwind of a snow fence or ridge. Many of these two-dimensional models rely heavily on the work of Pomeroy, who developed a comprehensive snow-transport modeling system (the Prairie Blowing Snow Model) to simulate blowing- and drifting-snow processes on the Canadian Prairies (Pomeroy, 1988, 1989; Pomeroy and Gray, 1990, 1995; Pomeroy and others, 1993).

The requirements of three-dimensional snow-transport models, ones capable of producing spatially distributed maps of snow depth, have only recently been addressed. In a simulation of the three-dimensional snow distribution within a 68 km² Arctic Canada watershed, Pomeroy and others (1997) divided the domain into blowing-snow source and sink regions (or elements) based on topography and vegetation characteristics. Using a simplified version of the Prairie Blowing Snow Model, which can be driven with monthly mean climatological data, a mass-balance computation was performed over each landscape element, yielding an end-of-winter snow distribution and simulation of the associated snow-transport-related fluxes. In another study, Purves and others (1998) developed a simple, rule- and cell-based model of snow transport and distribution.

In this paper, we present a physically based numerical snow-transport model (SnowTran-3D) that can be used to simulate snow-depth evolution over topographically variable terrain. The model was developed for, and tested in, an Arctic-tundra landscape but is applicable to other treeless areas such as prairies and alpine regions. It builds on previous snow-transport studies and modeling efforts but also includes several improvements. First, it is not restricted to two dimensions (i.e. two-dimensional snow cross-sections or profiles) but produces spatially distributed maps of snow-depth or snow water-equivalence. It introduces the equations required to compute the saltating and turbulent-suspended mass flux in non-equilibrium transport conditions, such as occur when flow is accelerating or decelerating near a topographic obstruction. In addition to being able to handle wind speeds which vary over the domain, the model simulates snow deposition and erosion under conditions of variable wind direction over the domain (i.e. the model will handle convergent and divergent wind fields). As part of the simulation procedure, the model keeps track of each component of the transport (saltation, suspension and sublimation) as a function of time and can be used to assess the relative importance of these processes in the development of the final snow distribution. The model, which requires time-independent inputs of vegetation type and topography, and time-dependent inputs of air temperature, humidity, precipitation and wind speed, and direction, has been validated against 4 years of detailed snow-depth maps from a 2 km by 3 km domain in the foothills north of the Brooks Range in Arctic Alaska.

2. Model Description

The snow-transport model is fully three-dimensional, in that it is implemented in two horizontal dimensions (x and y), and evolves the snow and snow water-equivalent depth (the z dimension) over a topographically variable domain. To simulate the x and y snow-transport contributions, the model formulation makes use of an additional coordinate system, x ^* (m), defined by the wind-flow direction (increasing downwind). The snow fluxes are computed along x ^* and then separated into their x and y components. The model considers only non-equilibrium transport from the perspective of spatially accelerating and decelerating flow; non-equilibrium transport resulting from temporal wind-speed accelerations and decelerations are not accounted for. The model is applicable to horizontal domains ranging from several tens of meters to a few hundred kilometers on a side. The topography within the domain can vary from flat, to gently rolling, to highly varying, such as regions where flow separation might occur over sharp ridges, gullies or valleys.

Figure 1 illustrates the key input parameters (solar radiation, precipitation, wind speed and direction, air temperature, humidity, topography and vegetation snow-holding capacity), the key processes (saltation, turbulent suspension and sublimation) and the key outputs (spatial distribution of snow erosion and deposition) from the model. The six primary components of the snow-transport model are: (1) the computation of the wind-flow forcing field; (2) the wind shear stress on the surface; (3) the transport of snow by saltation; (4) the transport of snow by turbulent-suspension; (5) the sublimation of saltating and suspended snow, and (6) the accumulation and erosion of snow at the snow surface, a lower boundary that is allowed to move with time.

Fig. 1. Key features of the snow-transport model applied to topographically variable terrain.

a. Deposition and erosion

The foundation of this snow-transport model is a mass-balance equation which describes the temporal variation of snow depth at a point. Deposition and erosion, which lead to changes in snow depth at this point are the result of (1) changes in horizontal mass-transport rates of saltation, Q _s (kg m^-1 s^-1), (2) changes in horizontal mass-transport rates of turbulent-suspended snow, Q _t (kg m^-1 s^-1), (3) sublimation of transported snow particles, Q _v (kg m^-2 s^-1) and (4) the water-equivalent precipitation rate, P (m s^-1). Combined, the time rate of change of snow depth, ζ (m), is

(1)

where t (s) is time; x (m) and y (m) are the horizontal coordinates in the west–east and south–north directions, respectively; and ρ _s and ρ _w (kg m^-3) are the snow and water density, respectively. Figure 2 provides a schematic of this mass-balance accounting. Equation (1) is solved for each individual grid cell within a domain and is coupled to the neighboring cells through the spatial derivatives (d/dx, d/dy).

Fig. 2. Schematic of the snow-transport model mass-balance computation.

b. Surface shear stress

The shear stress that wind exerts on the surface largely determines whether the wind speed is sufficient to transport snow. This shear stress is formulated in terms of the wind-shear (or friction) velocity, u _* (m s^-1), which is given by

(2)

where u _r (m s^-1) is the wind speed at reference height z _r (m), z ₀ (m) is the surface roughness length, and κ is von Kàrmán’s constant. If the wind-shear velocity exceeds the threshold shear velocity, u _*t (m s^-l), of the snow cover and, if snow is available, saltation will begin.

The availability of snow for transport is determined by defining a snow-holding capacity, C _v (m), for each vegetation type or land-surface classification within the domain. This capacity is a function of vegetation height and density, and/or microtopographic relief. Precipitation is converted to snowfall using the snow density and accumulates until the snow depth exceeds the vegetation snow-holding capacity. Once exceeded, any additional snow is available for wind transport.

To account for regions within the domain where the snow depth does not completely cover the vegetation, z ₀ is determined by linearly weighting the vegetation roughness length, z _{0_}, and the roughness length defined for snow, z _{0_snow}, according to the depth-fraction of vegetation covered by snow, F _s

(3)

where, in this case, the snow depth, ζ, is less than or equal to the vegetation snow-holding capacity. C _v. Thus, the roughness length, z ₀, is given by

(4)

This allows the computation of a spatially continuous wind field, including regions in the domain where the vegetation snow-holding rapacity has not been reached.

Given sufficient winds and snow depths above the vegetation snow-holding capacity, snow will begin to saltale. When it docs, the wind-shear velocity, u _* (m s^-1) and the roughness length, z ₀, are coupled (Owen, 1964; Tabler, 1980) and described by Equation (2), and

(5)

where g (m s^-2) is the gravitational acceleration. To determine whether snow is saltating, we first assume that saltation is present and then solve Equations (2) and (5) for u _* and z ₀ using Newton–Raphson iteration. If the resulting u _* exceeds a threshold shear velocity, u _*t, then saltation is present. If not, z ₀ = z _{0_snow} is defined for the stationary snow cover and u _* is computed using Equation (2) alone. In contrast to Pomeroy (1988), who was interested in agricultural fields containing large fractions of wheat stubble, we have assumed, largely due to lack of data, that stalks of vegetation protruding above the snow-holding-capacity height have a negligible influence on the roughness-length computation.

c. Saltation

Under equilibrium conditions, Pomeroy and Gray (1990) found that the saltation-transport rate, Q _{s_max} (kg m^-1 s^-1), can be described by

(6)

where ρ _a (kg m ³) is the air density. In general, u _* will vary in space and time, so Q _{s_max} will also vary. In addition, Q _{s_max} will be achieved only if there is sufficient fetch distance to entrain the full amount of snow that the wind speed, via u _*, is capable of carrying.

For the non-equilibrium (and equilibrium) conditions of interest in the current study, the saltation-transport rate, Q _s (kg m^-1 s^-1), is assumed to follow

(7)

and

(8)

where x ^* (m) is the horizontal coordinate in a reference frame defined by the wind-flow direction (increasing downwind); Δx ^* (m) is horizontal grid resolution; μ is a non-dimensional sealing constant, and f (m) is the equilibrium-fetch distance, assumed to be somewhere between the 300 m suggested by Takeuchi (1980) and the 500 m adopted by Pomeroy and others (1993).

For the case of decreasing wind-shear velocity with distance downwind, Equation (8) sets the saltation flux to be the minimum of either (a) that which can be supported by the wind-shear velocity, or (b) the saltation flux available from the upwind grid cell. For the case of constant or increasing wind-shear velocity with distance, Equation (7), discretized as

(9)

considers the upwind saltation flux (at x ^* – Δx ^*) and adds a bit more to it, to gel the flux at x ^*. The amount added decreases as Q _{s_max} (x ^*) is approached, until the fetch distance, f, has been reached or exceeded, at which point the saltalion flux becomes constant at Q _s(x ^*) ≈ Q _{s_max}; thus, Equation (6) is recovered as equilibrium is reached. To solve Equations (8) and (9) (and Equation (1), as discussed later), the saltation flux at upwind boundaries of the domain must be provided. At such a boundary, we assume that sufficient fetch and snow supply is available to reach equilibrium, and the inflow boundary condition is given by Equation (6). Because non-equilibrium transport resulting from temporal wind-speed accelerations and decelerations are not considered, initial conditions in the form of transport fluxes are not required and the transport fluxes at one time-step do not influence the transport at the next time-step.

As a further illustration of the non-equilibrium saltation-transport model behavior, consider an example where the saltation flux at the upwind boundary is zero, Q _s(x ^* = 0) = 0; such a condition might occur at the downwind edge of a drift trap. With this boundary condition and a constant wind-shear velocity (Q _s(x ^*) = Q _{s_max} = constant), Equation (7) has the solution

(10)

which is shown in Figure 3, where μ = 3.0 has been chosen such that Q _s(x ^* = f = 500 m) equals 95% of Q _{s_max}.

d. Suspension

Saltation must be present in order to have turbulent-suspended snow particles, and the saltation-transport rate provides the lower boundary condition for determining whether, and the degree to which, suspended transport occurs. The transport rate of turbulent-suspended snow, Q _t (kg m^-1 s^-1), is given by

(11)

where z is the height coordinate, φ _t (kg m^-3) is the mass concentration of the particulate cloud, u (m s^-1) is the wind velocity given by Equations (2) and (5), and the integration limits are the top of the saltalion layer, h _* (m), and the top of the turbulent-suspension layer, z _t (m), where the concentration is zero.

Fig. 3. Heterogeneous progression of saltation transport given by Equations (7) (and (9) ) for conditions of μ — 3.0, constant surface-shear velocity and an upwind boundary condition flux (at x^* = 0) equal to zero. For this case, 95% of equilibrium transport occurs at a fetch distance of 0.5 km.

The concentration of blowing snow within the turbulent-suspension layer, φ _t, is defined according to Kind (1992) as

(12)

where φ _r is the mass concentration at some reference level, z _tr (m) near the surface: s (m s^-1) is the particle-settling velocity, assumed to be 0.3 m s^-1 (Schmidt, 1982); and φ _* is a concentration scaling parameter. This formulation assumes that the concentration is in local equilibrium with the reference-level concentration and does not take into account any advective or non-steady affects. The ratio of the scaling parameter to the reference mass-concentration is approximated by

(13)

where β is a coefficient assumed to be 0.5 (Kind. 1992).

To solve Equation (12), the concentration reference level, z _tr (m), is defined to be the top of the saltation layer, at height h _* (m). This approach ensures continuity at the saltation/turbulent-suspension interface, and provides a method to define the lower boundary condition of the turbulent-flow regime based on our understanding of the conditions within the saltation layer. The saltation-layer height is estimated following Greeley and Iversen (1985),

(14)

Pomeroy and Gray (1990) suggested that the horizontal particle velocity within the saltation layer, u _p (m s^-1), is constant with height and approximated by

(15)

Under the further assumption that the saltation-layer mass flux is also independent of height, the saltation-layer reference-level mass concentration, φ _r (kg m^-3), is given by

(16)

where Q _s is given by Equations (8) and (9). Figure 4 illustrates the vertical profile of mass concentration for both saltation and suspended transport, for the case of u ^*t = 0.25 m s^-1, Q _s = Q _{s_max} from Equation (6), and u _* ranging from 0.3 m s^-1 to 1.0 m s^-1 in steps of 0.1. While the height-independent saltation-layer mass-concentration is an oversimplification (Kobayashi, 1972; Takeuchi, 1980; McKenna-Neuman and Nickling, 1994), this approximation is used until a more realistic description is available.

Fig. 4. Vertical profiles of mass concentration for saltation and turbulent-suspended transport, for the case of u_*t = 0.25 m s^-1, Q _s = Q _{s_max} from Equation (6), and a variable wind-shear velocity. The increase in saltation-layer height with increasing wind-shear velocity is evident in the profiles.

Figure 5 shows the variation of saltation and turbulent-suspension transport with wind-shear velocity, for the case of u _*t = 0.25 m s^-1 and Q _s = Q _{s_max} from Equation (6). In this example, suspended transport is zero between wind-shear velocities of 0.25 and 0.31, and saltation is the only transport mechanism. At a wind-shear velocity of 0.44, saltation and suspension transport become equal, and above that value suspension rapidly begins to dominate transport.

e. Sublimation

During transport by saltation and turbulent suspension, snow particles lose mass to the atmosphere through sublimation. In addition, the snow-cover surface can exchange moisture with the atmosphere by sublimation under non-transport conditions, but at a much lower rate than during transport. In this study, we only consider the snow mass lost during snow transport by the wind. The sublimation rate of transported snow, per unit area of snow cover, Q _v (kg m^-2 s^-1), is given by

(17)

where Ψ (s^-1) is the sublimation-loss rate coefficient, φ (kg m^-3) is the vertical mass-concentration distribution, and the integration limits are from the snow surface through the saltating and turbulent-transport regimes. If the saltation-layer contribution (z = 0 to z = h _*), where we have assumed constant flow properties with height, is separated from the turbulent-suspension layer contribution, Equation (17) becomes

(18)

where here the subscripts s and t refer to the saltation and turbulent-suspension layers, respectively, and the saltation-layer mass concentration, φ _s (kg m^-3), is given by Equation (16).

Fig. 5. Variation of saltation and turbulent-suspension transport with wind-shear velocity for the case of u_*t = 0.25 m, and Q _s = Q _{s_max} from Equation (6).

Our formulation for the sublimation rate of wind-transported snow, Q _v, follows that of Schmidt (1972, 1991), Pomeroy and others (1993) and Pomeroy and Gray (1995). The sublimation-loss-rate coefficient, Ψ, describing the rate of particle mass-loss as a function of height within the blowing and drifting-snow profile, is a function of (1) temperature-dependent humidity gradients between the snow particle and the atmosphere, (2) conductive and advective energy-and moisture-transfer mechanisms, (3) particle size, and (4) solar radiation intercepted by the particle. In addition, we assume that (a) the mean particle size decays exponentially with height, (b) the relative humidity follows a logarithmic distribution, decreasing with height, (c) the air temperature in the snow-transport layer is well mixed and constant with height, (d) the variables defined within the saltation layer are constant with height, and those in the turbulent-suspension layer vary with height, and (e) the solar radiation absorbed by snow particles is a function of the solar-elevation angle (latitude, time-of-day, day-of-year) and fractional cloud cover, σ ₀. The details of this formulation are given in the Appendix.

f. Wind field

To drive the snow-transport model, a reference-level wind-flow field, u _r, over the domain is required at each model time-step. This wind field can be generated by several methods, including (in decreasing order of complexity or difficulty), (1) applying a physically based, full atmospheric model (e.g. the Regional Atmospheric Modeling System (RAMS); Pielke and others, 1992), or a boundary-layer circulation model (e.g. Liston and others, 1993b; Liston, 1995), both of which satisfy all relevant momentum and continuity equations, (2) applying an atmospheric model in which only mass continuity is satisfied (e.g. Sherman, 1978; Ross and others, 1988), (3) interpolation using extensive (windward-and lee-slope) wind speed and direction observations, or (4) interpolation using wind-speed and direction observations in conjunction with empirical wind-topography relationships (e.g. Ryan, 1977).

To minimize computational expense, the current study uses method (4). Observed wind speeds in the domain are interpolated to the model grid and then the wind-speed field, u _r, is modified to account for topographic variations by multiplying by an empirically based weighting factor,W,

(19)

where Ω_s and Ω_c are the topographic slope and curvature, respectively, in the direction of the wind, and γ _s and γ _c are positive constants which weight the relative influence of Ω_s and Ω_c on modifying the wind speed. The slope and curvature are computed such that lee and concave slopes produce Ω_s and Ω_c less than zero, and that windward and convex slopes produce Ω_s and Ω_c greater than zero. Thus, lee and concave slopes produce reduced wind speeds, and windward and convex slopes produce increased wind speeds. Values of γ_s = 11 and γ_c = 360 have been found to produce wind-speed variations, in response to topography, which are consistent with observations of wind-microtopographic relationships (Kreutz and others, 1955; Nägeli, 1971; Yoshino, 1975), and have been used in this study. In this application, because of the relatively gentle nature of the local topography, the wind direction is assumed to be unaffected by topography; a point verified by measurements of surface-drift features in our Arctic lest domain (König, 1997). The topographically modified wind field is used to compute the u _* distribution by applying Equations (2) and (5), and Q _{s_max} is computed for each gridcell using Equation (6).

The horizontal coordinate x ^*, used in Equations (2) through (19), is aligned with the wind direction. To simulate the snow-depth evolution on a rectangular grid, the snow fluxes are computed along x ^* and then separated into their x and y components. The wind-field, u _r (Fig. 6a), is broken into its vector components, u (east–west) and v (north south) (Fig. 6b and c). These are used to identify sub-domains in which the u and v components come from the same direction: cast or west sub-domains for u, and north or south sub-domains for v (Fig. 6b and c). The magnitudes of the wind-velocity components are used to partition Q _{s_max} into an x and y component such that the vector sum of the u and v contributions equal Q _{s_max}

(20)

and

(21)

Equations (8) and (9) are then applied in the wind-flow direction for each sub-domain, where Q _{s_max} Equations (8) and (9) is replaced by Q _{s_max} (x) or Q _{s_max} (y) (Equations (20) or (21)), for the x and y directions, respectively. These components are then used in the solution of Equation (1).

3. Computational Results

The snow-transport model was applied to a 6 km² area that includes Imnavait Creek, in Arctic Alaska (Fig. 7). Imnavait Creek is a 2.2 km² basin located between the headwaters of the Kuparuk and Toolik Rivers at 68^° 37^′ N, 149^° 17^′ W, and an elevation (if approximately 900 m. The topography, physiography and vegetation of the watershed are representative of the area north of the Brooks Range and south of the Arctic Coastal Plain (Wahrhaftig, 1965). The vegetation covering the site is composed of low-growing sedges and grasses roughly 15 cm in height, with occasional groupings of taller willows, approximately 40 cm high, located in hillside water tracts and valley bottoms (Fig. 7). Tussock-tundra covers much of the area, with swampy features in the valley bottoms, and dry rocky outcrops on the exposed ridges. The topography is characterized by gently rolling ridges and valleys having wavelengths of 1–2 km and amplitudes of 25–75 m (Fig. 7). These ridges are aligned nearly north–south. Also, within the domain are several more-pronounced topographie features having much shorter and steeper slopes (up to 30^° slopes over distances of a few lens of meters), such as the “S2-Drift” (Fig. 7).

Fig. 6. (a) Example two-dimensional (x and y) gridded convergent/divergent wind-flow field. (b) u component of example flow field, with entire domain having flow from the east. (c) v component of example flow field, with two sub-domains, one having flow from the south (top sub-domain) and one having flow from the north (bottom sub-domain).

Muriel simulations were run from 1 September through 30 April, for the winters 1986–87, 1987–88, 1988–89 and 1989–90. In the discussions which follow, these periods are referred to as 1987, 1988, 1989 and 1990, respectively. The simulations were run with a 20 m by 20 m grid spacing and a 1 day time-step. Table 1 summarizes the user-defined constants used in the simulations.

Fig. 7. Vegetation and topography for the 2 by 3 km Imnavait Creek simulation domain. Topographic contour interval is 5 m. Also shown is the location of the Imnavait watershed boundary. (Data courtesy of N. Auerbach and S. Walker, Institute of Arctic and Alpine Research, Boulder, Colorado.)

a. Atmospheric forcing

Climate data from two meteorological towers in Imnavait Creek were obtained from the Water Research Center, University of Alaska, Fairbanks, and used to provide atmospheric forcing for the model simulations. A continuous, daily meteorological record of air temperature, humidity, and wind speed and direction, was produced by merging data from one tower to another if one lower failed, and by linearly interpolating between data values on the few occurrences when both meteorological towers failed (Fig. 8a–d; Table 2).

Table. 1. User-defined constants used in model simulations

Reliable (or even unreliable) winter precipitation measurements for the Arctic in general, and the test domain in particular, do not exist (Benson, 1982). To generate a daily precipitation dataset for the model simulations, it was necessary to adopt the following procedure: snow precipitation was assumed to fall when the air temperature was below-freezing and the relative humidity was greater than 80%. For each day meeting these criteria, the number of humidity percentage points above 80% was divided by the total number of percentage points above 80% for the entire winter. This produced a snow-precipitation factor equal to the fraction of total winter snowfall that occurred on each day. The daily factor, when multiplied by an estimate of the total end-of-winter precipitation, produced the precipitation amount for the day. Using this precipitation index factor reduced the problem of defining daily winter precipitation to one of only-needing to know the total winter precipitation. Unfortunately, that too was unknown, but we did know, from our field observations, the end-of-winter snow water-equivalent depth. So, the model was run in an iterative “data assimilation” mode, where the total winter precipitation was adjusted until the model-produced end-of-winter snow water-equivalent depths closely matched those measured in Imnavait Creek (Table 3) (Kane and others, 1991; Hinzman and others, 1996; Benson, unpublished data). In the simulations, precipitation was assumed to be spatially uniform over the domain. The total winter precipitation values necessary to produce the spatially averaged snow water-equivalent depths given in Table 3 are included in Table 2. In addition, the daily snow-precipitation forcing used in the model simulations is included in Figure 8a–d.

We found the model to be relatively insensitive to the value chosen for the threshold relative humidity (80%), which effectively determined when it snowed. In the Arctic, low winter temperatures mean that the shear strength of new snow is relatively low long after it snows. Thus, the timing between precipitation and wind events is not as critical as in warmer climates, where snow that accumulates a few days (or even hours) before the wind blows, may have strengthened to the point that it is unavailable for transport at naturally occurring wind speeds. In applications where the snow shear strength increases significantly following precipitation events, the method we use to determine precipitation timing and quantity would be inappropriate.

b. Observed snow distributions

Observed end-of-winter snow distribution maps for 1987, 1988, 1989 and 1990 for the domain shown in Figure 7 are presented in Figure 9. These maps, part of a sequence that began in 1985, have been prepared from a combination of ground-based depth and density measurements, and aerial photographs (Liston, 1986; Hinzman and others, 1996; König 1997; Benson, unpublished data). As an example of how these maps were developed, in 1985 three sets of oblique aerial photographs were taken: the first set was taken prior to any melting, when the snow depth was at a maximum; the second set was taken 6 days into the melt when 90% of the site was snow-covered; and the third set was taken 13 days into the melt when 20% of the site was snow-covered. For each of the three sets of aerial photographs, maps were drawn outlining the snow–vegetation boundaries. These three maps were merged into a common snow-pattern map, to which the snow water-equivalent depth observations were added. Analysis of the ground-based observations and the snow-pattern map showed that the map could be used to extrapolate snow-depth data into regions without direct measurements.

The accuracy of the observed snow-distribution maps (Fig. 9) can be variable, depending on the number of ground-based measurements and the number and quality of aerial photographs. In general, boundaries on the maps are somewhat simplified and idealized, though actual snow water-equivalent depth values within mapped areas are assumed fairly accurate. Small-scale features, like hollows and depressions that might hold deeper snow, tended to fall below the spacing of depth and snow water-equivalent measurements and are not depicted in the maps. For simplicity in discussion and comparing model results to measured values, we refer to these maps as the “measured or observed snow distribution” but, in reality, the true distribution of snow is not known. In some cases, we suspect that the measured snow-distribution maps contain errors in boundary-locations and depths that are greater than the errors within the modeled snow distributions.

Fig. 8. (a) September 1986 through April 1987 daily average atmospheric forcing data of air temperature, relative humidity, wind speed and direction, and snow water-equivalent precipitation used in the model simulations. (b) Same as Figure 8a, except for the period September 1987 through April 1988. (c) Sam as Figure 8a, except for the period September 1988 through April 1989. (d) Some as Figure 8a, except for the period September 1989 through April 1990.

Table. 2. Summary of atmospheric forcing used in model simulations. All values are compiled from daily averaged meteorological observations

Table. 3. Observed and model-simulated end-of-winter, snow water-equivalent depths (cm)

c. Modeled snow distributions

The simulated end-of-winter snow water-equivalent distribution for 1987, 1988, 1989 and 1990 are shown in Figure 10. Over the gently rolling topography that comprises most of the domain, the model simulates deeper snow on north- and east-facing slopes, with noticeably thinner snow on the south- and west-facing slopes. Since the wind is generally from the south or southwest (Table 2; Fig. 8a–d), these correspond to lee and windward slopes, respectively. The measured snow-distribution maps (Fig. 9) confirm that this is the pattern that develops. Even the trends of snow water-equivalent contours resulting from model simulations are similar to those of the measured distributions. Deep drifts (say over 30 cm) are also correctly simulated by the model. The model indicates that most of these lie in the lee of the steep topography located in the northern quarter of the domain, a fact confirmed by the field measurements. For example, the model predicts an excess of 30 cm snow water equivalent at a location where we bave observed drift accumulation since 1985. In this location (labeled “S2-Drift” in Figure 7), a large wedge-shaped accumulation, ranging up to 2 m thick, forms every winter. Overall, the model indicates that these drift accumulations of over 30 cm snow water equivalent covered 1%, 1%, 2% and 1% of the domain in 1987, 1988, 1989 and 1990, respectively. These areal coverages compare well with the 2%, 1%, 2% and 1 % indicated by the observed distributions.

The model simulations can also be used to assess the various factors producing the end-of-winter snow distribution. The individual domain-averaged, winter-total contributions of precipitation, saltation, turbulent-suspension and sublimation are shown in Table 4. The sign of the values in Table 4 indicates either accumulation (+) or loss (–) of water from the domain. Precipitation minus sublimation, minus the transport out of the domain by saltation and suspension, equals the depth of snow on the ground at the end of winter. As an example, in 1987,71 % of the winter precipitation remained on the ground, 22% was returned to the atmosphere through sublimation, and 5% and 2% was transported out of the domain by saltation and suspension, respectively. The spatial distributions of each of these water-balance components are displayed in Figures 11, 12 and 13.

Inconsistencies between modeled and measured snow distributions appear to arise, at least in part, from inaccuracies in the wind field computed by Equation (19). One type of inconsistency is where the modeled snow depth is correct but its location is shifted laterally with respect to the measured distribution. The precise character of increasing and decreasing wind speed as it flows over a subtle ridge is difficult to capture with the simple wind model. In contrast, the influence of the flow over a sharp ridge, such as the “S2”drift trap in the north of the domain, is easy to simulate and the model consistently locates the resulting deep drifts correctly. Equation (19) is also inadequate to simulate the wind field required to produce the observed deep accumulation located at approximately 2.3 km north and 1.7 km east, occurring in all four winters (cf. Figs 9 and 10). This accumulation forms because the wind coming from the south or southwest “feels” the relatively large topographic obstruction ahead and, as the streamlines begin to separate from the snow surface, the shear velocity is reduced and snow is deposited much like the accumulai ion upwind of a snow fence (Mellor, 1965; Tabler, 1980; Liston and others, 1993a, b). To simulate such flow features requires a more sophisticaied and computationally expensive wind model.

4. Analysis of Model Simulations

Based on data collected for the Imnavait Creek watershed between 1985 and 1993 (Kane and others, 1991; Hinzman and others, 1996), the period over which we ran our simulations includes two “average” accumulation winters (1987 and 1990), one “light” accumulation winter (1988) and one “heavy” accumulation winter (1989) (Table 4). For each of these winters we have calculated, from model simulations, the area-averaged snow water equivalent remaining on the ground at the end of winter, and how much snow (in snow water-equivalent units) has been lost from the domain by sublimation, saltation and turbulent-suspension (Table 4). We have also mapped how these gains and losses were distributed over the domain (Figs 11, 12 and 13). The results indicate that, averaged over the domain, 9–22% of the total-winter precipitation was returned to the atmosphere by sublimation. This range is comparable to values found by Pomeroy and Gray (1995) for the Canadian Prairies and is large enough to impact the end-of-winter distributions significantly. In 1987, for example, there was 27% greater precipitation than in 1990, but at the end of winter, both the area-averaged snow water equivalent (Table 4) and the snow-distribution patterns (Figs 9 and 10) were similar. The primary cause leading to these similar end-of-winter conditions was sublimation losses, which were 2.9 times greater in 1987 than in 1990. The greater amount of sublimation in 1987 was the result of more days with average wind speeds in excess of 10 m s^-1 (Table 2). As a further example, the precipitation in 1988 and 1990 were nearly equal but, with sublimation losses running nearly twice as high in 1988, the 1988 end-of-winter snow water-equivalent values were notably lower. The 1988 sublimation losses were large enough to change the winter snow-accumulation classification from “average” to “light” (Table 4). In this case, the number of days having average wind speeds in excess of 10 m s^-1 was not greatly different between the two winters, but in 1988 the daily average wind speed stayed almost continuously above 5 m s^-1 through January and part of February (Fig. 8b), while in 1990, higher wind speeds were limited to isolated events (Fig. 8d). During the same period, the 1988 temperature was higher and the relative humidity lower, both conditions that favor sublimation.

Fig. 9. Observed snow water-equivalent distributions (cm) for the winters 1987, 1988, 1989 and 1990 in the Imnavait Creek area (Liston, 1986; Hinzman and others, 1996; Benson, unpublished data ).

Fig. 10. Model-simulated end-of-winter snow water-equivalent distributions (cm) for the winters 1987, 1988, 1989 and 1990. Also shown is a 1 km box in the 1987 panel outlining the domain coincident with Figure 14.

Model simulations indicate that maximum sublimation losses were located on south- and west-facing slopes (Fig. 11), in the same locations where saltation and suspension-transport losses were highest. Since the sublimation rate is tied to both the saltating and suspended mass flux through Equation (18) (also see the Appendix), this is expected; these windward slopes are where transport processes are most vigorous (Figs 12 and 13). Data collected during the winter of 1995–96 from a 1 km by 1 km section within the model domain support the simulations. For this area, the snow-depth distribution was mapped on 6 November 1995, and again on 26 March 1996. The difference in depth between the two observation times is shown in Figure 14; the domain covered by Figure 14 corresponds to the 1 km box drawn in the 1987 panel of Figure 10. The observed pattern in Figure 14 is strikingly similar to those found in Figure 10, and indicates that on the windward slope little change in depth occurred during the 5 month period, while on lee and valley-bottom slopes, increases in snow depth as large as 55 cm occurred. During the same period, snowfall was about average with numerous snowfall events. The minimal windward-slope depth changes and small downwind (to the north-northeast) snow-depth enhancement, in conjunction with the model results, suggests that the majority of winter precipitation falling on windward slopes was lost by sublimation.

Table. 4. Domain-averaged, winter-total, snow water-equivalent quantities (cm) from model simulations

The saltation-transport patterns (Fig. 12) highlight the importance of topography in controlling this process. During the years 1987, 1988 and 1989, snow was generally removed from south- and west-facing slopes and deposited on lee slopes, collecting in large amounts in the “S2” and other drift-accumulation areas. In contrast, in 1990, when the wind tended to blow more from the southeast (Table 2) and roughly parallel with the ridge and valley lines, there was virtually no change in depth due to saltation over most of the gently rolling topography of the domain (Fig. 12). At the same time, the “S2” and other drift areas continued to accumulate snow. These drift-accumulation areas are much smaller than their upwind supply zones, and are thus able to acquire significant accumulations while the supply-area depths change only slightly.

The patterns of snow-depth change resulting from suspended transport (Fig. 13) are similar to those of saltation, with appreciable amounts transported in 1987, 1989 and 1990. In many regions of the domain, for these 3 years, the suspended-transport amounts exceed those of saltation (cf. Figs 12 and 13). As indicated by Figure 5, this behavior is expected at higher wind speeds and is consistent with findings by Pomeroy (1989), who concluded that, as wind speed increases, suspension becomes the dominant transport mechanism. In 1988, when the daily average wind speed never exceeded 10 ms^-1 (Table 1), transport due to suspension was minimal (Fig. 13).

5. Discussion

The snow-transport model has been used to map the redistribution of snow (or snow water equivalent) by wind over a realistic landscape consisting of variable vegetation, ridges and valleys. The resulting end-of-winter snow distribution maps (Fig. 10) compare well with ones based on measurements (Fig. 9). The model could produce accurate maps, yet still have errors in its sublimation, suspension and saltation routines. However, these routines are based on well-researched process studies and, where we have been able to compare individual simulated transport components with observations (e.g. Fig. 14), the agreement is good. Where differences between the measured and simulated maps exist, they appear, in part, to be largely due to inaccuracies in the simulated wind field. High-quality atmospheric forcing (precipitation, wind speed and direction, air temperature and relative humidity) are required to simulate the observed snow distributions but are not always available.

Each component of the model is necessary to produce the correct end-of-winter snow distribution but in the Arctic the sublimation routine is crucial. In some regions of Imnavait Creek, sublimation removed as much as half of the snow on windward slopes (Figs 10 and 11). This suggests that a substantial amount of the large-scale distribution patterns shown in Figures 9 and 10, and possibly throughout the Arctic, are the result of sublimation. If so, the finding has practical implications regarding the origin of perpetually snow-free zones in the Arctic. For example, the snow-free valleys of the Brooks Range in Alaska (e.g. the Itkillik, the Anaktuvak and the Atigun), where katabatic down-valley winds are frequent, may be the result of blowing-snow sublimation, not a lack of winter precipitation nor the transport of snow out of the valley. In the coastal regions of Arctic Alaska, where wind speeds are generally much higher, sublimation losses are expected to be greater than those found in Imnavait. As an example, increasing the wind speed from 4.0 m s^-1 (the measured Imnavait 4 year winteraverage) to 5.1 m s^-1 (the Barrow winter average; Hare and Hay, 1974), increases the 4 year average sublimation loss from 17% to 27% of the winter precipitation.

Using the model to partition the end-of-winter snow distribution into its individual components also provides insight into deep-drift accumulation processes. Drifts form even in years when saltation or suspension transport are relatively small. This is because drifts occupy a small area compared to their source area. In addition, once the snow is trapped in the drift, it suffers virtually no loss from sublimation (Fig. 11), ensuring the snow depth remains high.

Topography dominates the wind redistribution of snow in a way that makes the general snow-distribution patterns more sensitive to changes in wind direction than wind speed. The distribution patterns (Figs 9 and 10) are similar as long as the average winter wind comes from a generally fixed direction but any significant shift in wind direction relative to topography alters the distribution pattern. As part of model-sensitivity simulations we have imposed wind-direction perturbations ranging from ±15^° to ±180^° away from the observed wind directions. In these simulations, the patterns of snow depth, saltation, suspension and sublimation can all be modified significantly, depending on the relationship between the new wind direction and topography; the most dramatic snow-pattern changes occur when the new wind direction changes a windward slope to a lee slope, or a lee slope to a windward slope. Application of the model suggests that a realistic computation of the wind-speed and direction field is a dominating factor leading to a successful snow-distribution simulation, and that the simulated snow distribution closely resembles the integrated winter wind field occurring above some threshold wind speed.

Fig. 11. Total winter sublimation, in snow water-equivalent depth (cm) for the simulations given in Figure 10.

a. Model deficiencies

The model performs reasonably well in our test cases but it is still a simple approximation to a complex natural system. Model deficiencies can be divided into two areas: those that arise because of limitations and constraints on input data and those that arise due to erroneous assumptions, oversimplifications or errors in the model snow-transport physics. The former limitations are thought to be more severe than the latter.

As an example, if our topographic data do not resolve a hill or bump, the model cannot accumulate snow in the lee of that bump. Likewise, if our meteorological input fails to resolve brief but intense periods of snowfall or high winds, then the actual drift accumulation or erosion will not be simulated. In addition, the simple model used to distribute wind over the domain may be inadequate for many topographic profiles. A more detailed and sophisticated wind model, one that could handle flow over a large range of topographic configurations, including three-dimensional flow over sharp ridges, peaks and valleys, or flow around solid and porous obstructions such as trees, snow fences and rock outcrops, would improve the results. However, the trade-off is increased computational complexity, and spatially and temporally detailed input data would still be required.

Fig. 12. Total winter snow-depth changes resulting from saltation, in snow water-equivalent depth (cm) for the simulations given in Figure 10.

The choice of model grid size will influence the resulting simulation. In general, the physics and assumptions adopted in the model will not be violated if the grid size is reduced or increased, say within a range of 1 m to a few hunched meters, but the different resolutions will have certain consequences. Increasing the grid size in the current study leads to reduced definition of the “S2-type” drifts. If the resolution exceeds the general scale of that feature, it will be lost by the simulation entirely, or at least smoothed to the extent that it is not recognizable in its original form. For the domain considered here, a computational grid of 100 m would be quite appropriate and produce outputs similar to those presented here, with the expected loss of the deeper drifts (e.g. we are running the model successfully for a 230 km by 85 km domain, using a 100 m grid, that encompasses the Kuparuk River basin in Arctic Alaska). When a high-resolution, 1 m grid is used over an area like the “S2” drift, a sharp decrease in shear velocity will occur just over the topographic lip. Large snow deposits will form in the lee of the lip, and in time this snow may be higher than the elevation of the gridcell immediately upwind. In Nature, the wind would notice this snow bump and transport it farther downwind, but in the model this must be done by including a parameterization that redistributes the snow downwind. An alternative approach is to adjust the wind field as the “snow topography” evolves (Liston, 1991; Liston and others, 1993a; Sundsbø, 1997). This would replicate the natural system more accurately and would be a significant improvement in places where large changes in surface topography occur due to snow deposition. As snow accumulates, the drift trap should evolve towards an equilibrium profile. Observational studies have found an equilibrium snow-embankment angle of approximately 17% (Finney, 1939; Tabler, 1975b), below which snow no longer accumulates, suggesting the wind-shear velocity is uniform over the embankment. Such a feature should be reproduced as part of a wind-modeling system that adjusts to the “snow topography”.

Fig. 13. Total winter snow-depth changes resulting from suspension, in snow-water-equivalent depth (cm) for the simulations given in Figure 10.

Fig. 14. The difference between snow depth measured on 6 November 1995 and 26 March 1996 for a domain corresponding to the 1 km box drawn in the 1987 panel (Fig. 10).

Problems with our physical formulations may also adversely affect the simulations. These include (in decreasing order of severity):

(1) Actual snow density and threshold shear velocity vary both temporally and spatially but, as a first approximation, the model assumes both to be constant in time and space (Table 1). In order to test the effect of this simplification, a threshold shear-velocity parameterization was implemented. This parameterization kept track of the amount and location of relatively new snow and redistributed or old snow, assigning different threshold shear velocities to each. In this scheme, the new-snow depth increased with each precipitation event, then decreased following the event according to a predefined densification rate. As the snow densified, its threshold shear velocity increased. The parameterization did not lead to any measurable improvement in the accuracy of the modeled snow cover and the increase in computational complexity did not seem justified. This result suggests that in the Arctic, it is acceptable to use a spatially and temporally constant threshold shear velocity, because the threshold shear velocity changes very slowly in the relatively low air temperatures. A more complex formulation would be required for model implementations in warmer climates.

(2) Natural fluctuations in wind speed and direction produce transport fluxes that may not be in equilibrium with the average wind speed. Our model does not account for such non-steady conditions. Related to this is the assumption that the concentration profile of the suspended snow load is in local equilibrium with the saltation-layer concentration. This is also inappropriate under certain conditions. Under high wind speeds, for example, when snow is transported over a sharp ridge, the turbulent suspended plume can be transported far beyond the ridge crest, while the saltation-layer concentration falls off rapidly with the decrease in wind-shear velocity. In the model, the reduced near-surface wind speed just over the ridge crest leads to unrealistically low particle concentrations in the plume above, and consequently unrealistically high snow depositions.

(3) As indicated in the Appendix, feedbacks between the sublimation rate and the humidity field are not accounted for. We intuitively expect that sublimation should lead to increased atmospheric moisture and this in turn will decrease the sublimation rate. However, using a model that accounts for these feedbacks, and additional interactions with turbulent moisture and momentum fluxes, R. Bintanja (personal communication, 1998) has found that the absolute-humidity increase, resulting from snow-particle sublimation, is compensated for by increased upward moisture transport out of the surface layer. Clearly, such interactions need to be addressed in a model attempting to simulate realistically the details of blowing-snow sublimation. In addition, these results suggest that the quantities produced by our simple sublimation submodel may be in error but magnitude and direction of the potential error is unknown.

(4) The model assumes a vertical distribution of mean particle size. In reality, the range of particle sizes within the snow-transport profile can vary by an order of magnitude (Budd, 1966; Schmidt, 1982). This spectrum of sizes is expected to influence many of the processes associated with turbulent kinetic energy and moisture transfers occurring during snow transport (Déry and Taylor, 1996) but we do not know how it might affect the modeled results.

6. Conclusions

A physically based numerical snow-transport model (Snow-Tran-3D) has been developed and used to simulate the redistribution of snow by wind in variable topography. The model describes interactions between the atmosphere, land and snow cover in treeless areas such as tundra, prairie and alpine landscapes. The model includes an accounting of (1) the capturing of snow by vegetation, (2) non-equilibrium transport by saltation and turbulent suspension, (3) the influence of saltation on surface-shear velocity, (4) sublimation of the blowing and drifting snow and (5) erosion and deposition.

Comparison of model simulations with 4 years of Arctic Alaska field observations shows the model is able to describe the general snow-distribution features. In addition, the model can be used to quantify moisture losses due to sublimation during snow-transport events. In the foothills north of the Brooks Range in Arctic Alaska, these losses were found to range between 9% and 22% of the precipitation input. In the coastal regions of Arctic Alaska, where wind speeds are generally much higher, the model suggests that sublimation will be greater.

In the current study, the threshold shear velocity has been assumed to be a pre-defined constant which does not evolve in time or space. While we know that this is a crude approximation to the natural system, attempts to introduce more complex and physically realistic parameterizations of this important variable did not lead to improved simulations, at least within the spatial and temporal resolution and quality of our validation datasets. One reason why we are able to use this simple approximation successfully is that, during the winter months, temperatures in Arctic Alaska are low and the snow maintains a nearly constant surface-shear strength. To apply this snow-transport model in a region where near-freezing temperatures are a regular occurrence, or to follow the detailed evolution of individual storm events, a more realistic formulation of the threshold shear velocity, which includes an accounting of temperature, wind, precipitation and snow-transport histories, would be required.

The ability of the snow-transport model to reproduce the observed snow distributions and moisture fluxes is strongly dependent upon the simulated wind-flow field. While the current wind-field representation was found to be adequate, the transport model would benefit from a more general wind model which would allow straightforward application to other, generally more complex, landscapes.

The success of this snow-transport model suggests that it is a viable complement to labor-intensive field snow-mapping methods. Using standard meteorological observations, and topography and vegetation datasets, in conjunction with limited end-of-winter mean snow water-equivalent observations, the snow-cover distribution can be simulated. The physically based nature of the model allows for its application in a wide variety of climates and allows the ability to address issues of changing climate. With only a few modifications to account for processes occurring at higher temperatures, the model should be applicable in a wide range of landscapes and regions throughout the world.