Wind is widely recognized as one of the dominant controls of snow accumulation and distribution in exposed alpine regions. Complex and highly variable wind fields in rugged terrain lead to similarly complex snow distribution fields with areas of no snow adjacent to areas of deep accumulation. Unfortunately, these complexities have limited inclusion of wind redistribution effects in spatial snow distribution models. In this study the difficulties associated with physically exhaustive wind field modeling are avoided and terrain-based parameters are developed to characterize wind effects. One parameter, , was based on maximum upwind slopes relative to seasonally averaged winds to characterize the wind scalar at each pixel location in an alpine basin. A second parameter, , measured upwind breaks in slope from a given location and was combined with an upwind application of to create a drift delineator parameter, D0, which was used to delineate sites of intense redeposition on lee slopes. Based on 504 snow depth samples from a May 1999 survey of the upper Green Lakes Valley, Colorado, the correlation of the developed parameters to the observed snow distribution and the effect of their inclusion in a spatial snow distribution model were quantified. The parameter was found to be a significant predictor, accounting for more of the variance in the observed snow depth than could be explained by elevation, solar radiation, or slope. Samples located in D0-delineated drift zones were shown to have significantly greater depths than samples located in nondrift zones. A regression tree model of snow distribution based on a predictor variable set of , D0, elevation, solar radiation, and slope explained 8%–23% more variance in the observed snow distribution, and performed noticeably better in unsampled areas of the basin, compared to a regression tree model based on only the latter three predictors.
In mountainous basins, snow distribution exhibits tremendous spatial heterogeneity (Elder et al. 1991; Doesken and Judson 1996; Luce et al. 1998; Balk and Elder 2000). Intrabasin snowmelt differences, which strongly influence discharge (Marks et al. 2001), basin ecology (Barron et al. 1993), water chemistry (Woolford et al. 1996), and hillslope erosion (Tarboton et al. 1991), can be largely dependent upon snow accumulation differences (Seyfried and Wilcox 1995; Luce et al. 1998). In order to accurately model intrabasin runoff it is essential that, whenever present, disparate distribution patterns be accounted for (Seyfried and Wilcox 1995). The redistribution of snow by wind has often been cited as one of the strongest influences on snow accumulation differences in alpine basins (e.g., Elder et al. 1991; Blöschl and Kirnbauer 1992; Luce et al. 1998).
The effects of wind on snow have been widely researched (e.g., Föhn 1980; Schmidt 1982; Tabler 1994), and the significance of wind in determining the spatial distribution of snow in open, windswept areas has been firmly established (e.g., Elder et al. 1991; Kane et al. 1991; Pomeroy 1991; Tabler 1994; Luce et al. 1998; Liston and Sturm 1998). While physically based, continuous models of snow redistribution have been developed for flat to gently rolling terrain (e.g., Pomeroy et al. 1993; Liston and Sturm 1998), extension of these findings to alpine terrain has been limited by the complexity of alpine wind fields. Moreover, these latter models are computationally intensive, hindering their functionality as an operational tool.
It is known that airflow responds to changes in topography following the principle of continuity (Berg 1977). Vertical flow constrictions, such as movement over a topographic obstacle, cause wind speeds to increase, thereby increasing scour rates and augmenting snow transport. Conversely, expansion of flow occurs downwind of topographic obstacles, leading to reduced wind speeds, decreased snow transport fluxes, and potential snow deposition (Berg 1977; Kind 1981; Barry 1992). If wind speeds at the topographic obstacle are strong enough and an ample slope break exists downwind of the obstacle, flow separation can occur, whereby the airstream is no longer in contact with the ground, forming a lee eddy. Often large snowdrifts, which can be an important component of alpine water storage (Barry 1992), will form in these eddies. It is also possible for a flow-separation eddy to develop at the base of a windward-facing topographic obstacle (Barry 1992; Whiteman 2000).
Recognizing the importance of wind in determining snow accumulation patterns, the link between topography and wind, and the comparative ease of terrain modeling, other researchers have also sought to parameterize the effects of wind redistribution through terrain analysis. Several researchers have taken existing terrain parameters and applied them to snow modeling. Curvature, the second derivative of a surface fitted to elevation, is a measure of the concavity or convexity of the terrain at a particular point or pixel and is commonly available as a function in Geographic Information System (GIS) software packages. Observationally, Golding (1974) and Woo et al. (1983) both found higher than average snow accumulations in concave regions such as valley bottoms and gullies, and lower than average accumulations in convex areas such as ridgetops. From a modeling perspective, Blöschl et al. (1991) achieved good results interpolating snow water equivalents as a function of slope and curvature, while Liston and Sturm (1998) applied topographic slopes and curvatures to adjust wind speed, thus controlling erosion and deposition rates within the most basic application of their snow-transport model, SnowTran-3D. Hartman et al. (1999) used the TOPMODEL-derived topographic similarity index or wetness index to redistribute snow from steep slopes to areas of topographic convergence and low gradients. The curvature parameter has a 360° sphere of influence, giving equal weight to both downwind and upwind features, while the Hartman et al. application accounts for upslope terrain regardless of its aspect relative to the wind. Winstral (1999) has shown, as Hartman et al. (1999) suggested, that terrain parameterizations capable of focusing on upwind conditions are preferred to global parameterizations (e.g., curvature and the TOPMODEL topographic similarity index).
Several researchers have developed terrain parameters explicitly designed for modeling snow redistribution. Cline (1992) analyzed both upwind and downwind terrain features to model the redistribution patterns of an early season snowfall event on Niwot Ridge, Colorado. Using parameters that described the topographic relationship between each cell and its immediate neighbors, he found that a terrain-based modeling approach was appropriate for modeling wind redistribution effects. Lapen and Martz (1993) developed two parameters to characterize shelter and exposure that were based upon analysis of cells located along cardinal direction vectors emanating from each grid cell. Similarly, they concluded that their simple empirical measures demonstrated the potential importance of topographic variables as surrogates for modeling wind-induced snow redistribution. Purves et al. (1998) developed a theoretical sheltering index, based upon slope and aspect relative to prevailing winds, to spatially distribute wind speeds for input to a blowing snow model. Though quantitative analysis was precluded due to limited snowfall during their study period, qualitative analysis based upon visual assessments of expected large-scale features showed good agreement.
The study described herein aimed to develop new means of quantifying terrain in a simple, yet spatially robust, manner in order to improve upon earlier snow redistribution modeling efforts. The robustness of the terrain parameters refers to an ability to search and analyze at varying spatial extents. The user can control the spatial orientation and lateral extent of terrain to be included in the analysis. We derived parameters that, at each grid cell, determined the degree of topographic shelter/exposure to prevailing winds, quantified upwind breaks in slope to locate potential flow separation zones, and assessed the exposure potential of the upwind slope breaks. While the parameters can also be used to analyze downwind topography to potentially distinguish flow separation zones on the upwind side of obstacles, preliminary work indicated that this dataset was not appropriate for such an analysis. In this study, we focused on the upwind terrain, seeking to characterize the wind scalar at each grid cell while also delineating sheltered sites prone to enhanced deposition in lee-eddy zones.
Based upon data collected from an intensive snow survey conducted in the upper Green Lakes Valley (GLV), Colorado, we assessed the significance of the relationships between the topographically derived redistribution parameters and the observed snow distribution. We also tested the significance of some of the commonly applied predictors of snow distribution—elevation, radiation input, and slope—for comparison with the redistribution parameter results. Elevation is commonly cited as a strong positive influence on snow accumulation rates (McKay and Gray 1981), solar radiation is an important positive component of the snow cover energy balance whose spatial heterogeneity is accentuated by rugged terrain and low wintertime solar incidence angles (Elder et al. 1991; Marks et al. 1992), and slope acts as a control on downslope redistribution due to sloughing and avalanching (Blöschl et al. 1991; Elder et al. 1991). As a final step, a regression tree model of snow distribution that included the redistribution parameters as well as elevation, potential radiation input, and slope was compared to a like model based on only the latter three predictors in order to assess the effect of the redistribution parameters upon modeled snow distribution.
Earlier versions of the presented redistribution parameters were applied in regression models of late season snow distribution in Tokopah Basin (Winstral et al. 1998) and Blackcap Basin, California (Winstral 1999), and for a prior year in the Green Lakes Valley (Winstral 1999). In each application, the tested parameters were significantly related to the observed snow distribution and improved spatial models. In the California basins, however, the prevailing storm flow of southwest to northeast formed a strong correlation between wind and solar radiation loadings that precluded a definitive independent test of the parameters. In the Green Lakes Valley, prevailing winds have a strong westerly component presenting a best-case scenario for separation of wind and solar effects in a regression-based approach. While the results of the earlier Green Lakes Valley study were encouraging, the snow data available at the time exhibited a high degree of spatial aggregation that compromised some of the findings. In this current work, data collection was confined to the upper section of the basin, with samples obtained from all accessible terrain. The orthogonal relationship between wind and solar radiation loadings in this windswept basin, in combination with the spatially extensive dataset, presented a solid testing ground for validation of this terrain-based approach to modeling wind-induced snow redistribution.
2. Field methods
a. Study site
The Green Lake 4 study area is part of the Niwot Ridge (NWT) Long-Term Ecological Research (LTER) project site located approximately 35 km northwest of Boulder, Colorado (Fig. 1). The Green Lake 4 basin flanks the eastern side of the Continental Divide, functioning as a headwater catchment for North Boulder Creek, a primary water source for the city of Boulder. The research area is 2.25 km2 in size, has an elevation range of 3560–4090 m, and is entirely above treeline. The terrain is steep and rugged, with an average slope of 28°. Seasonal snow accumulation typically begins in October. Mean monthly minimum temperatures are below 0°C from October through May, with mean monthly maxima below freezing from November to April (Berg 1986). The continental location and high altitude combine to generate low-density snowfall events and atmospheric conditions that retard thermal densification of the surface snow, creating an environment where snow transport is rarely limited by snow cover surface conditions (Berg 1986). The basin's characteristic high winds, steep topography, and continental snow climate produce an extremely heterogeneous snow cover governed largely by wind redistribution.
b. Snow data
The snow survey data were collected, prepared, and made available for this study by the National Science Foundation (NSF) supported NWT LTER project and the University of Colorado Mountain Research Station (MRS). The snow survey was conducted as part of an annual assessment of snow distribution in the study basin and was not designed specifically for the study described herein. Depth samples were taken at 504 sites in mid- to late May 1999 (Fig. 1). The locations of the sample sites were recorded and differentially corrected using a global positioning system (GPS). Whenever possible, measurements were taken at approximately 50-m intervals. Steep terrain and avalanche danger, however, limited data collection to safe and accessible areas. The 1 May estimate of snowpack conditions in the Boulder Creek drainage was 112% of average (NRCS 1999).
The modeled dependent variable in this study was snow depth. Snow-water-equivalence (SWE), the state variable of primary interest to researchers and forecasters, is obtained from the product of snow depth and snow density. Snow density is substantially more conservative than snow depth, and the variability of snow density continues to decrease as the snowpack ripens into the late spring (Elder et al. 1991). We have used snow depth as a surrogate for SWE and have assumed that distributed SWE would follow the pattern of distributed depths in this mid- to late May modeling time frame.
c. Wind data
Wind data were likewise made available by the NSF-supported NWTLTER project and the MRS. Wind measurements were made at the D-1 meteorological station located on Niwot Ridge (3744 m), an exposed east–west-trending rise in the northeast corner of the Green Lake 4 basin (Fig. 1). The wind instrumentation is located approximately 9 m above the ground surface in an exposed location. Hence, there is little snow accumulation and little variability in the instrument height through the year. The beginning of the snow accumulation season in October and relatively low May wind speeds, combined with the presence of moist spring snow, typically limit significant blowing snow events to a 6.5-month period from mid-October to early May (Berg 1986). Thus, the wind data gathered for this study were restricted to the 15 October–7 May time period.
The wind direction data consisted of the daily resultant means of instantaneous recordings taken at 30-s intervals for the accumulation season. Unfortunately, wind direction data from the D-1 site were unavailable for approximately 1 month (14 February 1999–14 March 1999) during the accumulation season prior to the survey dates. Based on historic data, wind patterns during the mid-February to mid-March period are consistent with other wintertime observations. Figure 2 illustrates the characteristic wintertime dominance of western flow at the D-1 site during the study period. The mean resultant wind direction for the study period was 265°.
A method of monitoring blowing snow events was not available during the study period. A comparison with wind speeds registered during blowing snow events recorded with an elevated snow particle counter as part of Berg's (1986) study provided insight into the degree of wind-induced snow redistribution during the 1998/99 season. At a meteorological station located on Niwot Ridge, 3 km east of the D-1 site, 53% of the hours monitored over the course of two winters (1973/74 and 1974/75) recorded blowing snow events (Berg 1986). Using wind data taken from a similar instrument height (8.75 m), Berg found that blowing snow occurred during 22% of the hourly averaged wind speeds in the 7–8 m s–1 category, during 49% of 13–15 m s–1 winds, and during 75% of winds of 29–30 m s–1. Over the course of the 1998/99 blowing snow season, 71% of the days had a mean wind speed greater than 7.5 m s–1, 29% had a mean wind speed greater than 14 m s–1, and on only two days did a recorded wind speed fail to exceed 7.5 m s–1. These data attest to the fact that wind-induced snow redistribution occurs on an ongoing basis throughout the winter in this environment.
3. Modeling methods
Due to different areal extents of the available digital elevation models (DEMs), two different modeling scales were used in this study: 10- and 30-m cells. The western border of the 10-m DEM, at its nearest, is 150–200 m west of the Continental Divide, which defines the western boundary of the Green Lake 4 basin (Fig. 1). As the upwind region of analysis of the redistribution parameters grew, topographic data extending beyond the bounds of the 10-m DEM were required, necessitating the additional coverage provided by the 30-m DEM. The 10-m DEM was created by Analytical Surveys, Incorporated (ASI; Colorado Springs, Colorado), using a combination of black-and-white and color infrared aerial photographs to collect breaklines and mass points for capture of any feature larger than 10 m (Williams et al. 1999). The 30-m coverage consisted of a combination of the U.S. Geological Survey's (USGS) Monarch Lake and Ward Mountain DEMs.
An attempt was made to merge data from the USGS 30-m DEM with the ASI-developed 10-m DEM in order to supplement the high-resolution data with smaller-scale data beyond its boundaries. However, distinct and irregular differences in elevations along the boundaries of the independently derived DEMs precluded a fully satisfactory merging. Hence, the two raster products were used independently of one another, and the scale chosen for derivation of each redistribution parameter, along with the reasoning for each choice, is given with the parameter descriptions.
Elevation, slope, and net potential radiation, calculated for use in the spatial snow models and for comparisons with the redistribution parameters, were derived from the 10-m DEM. Elevation and slope were not affected by the proximal boundary of the 10-m DEM, while the radiation calculations were negligibly affected by terrain-reflected radiation emanating from outside the basin. Net potential radiation over the snowpack was characterized by a solar radiation index based on clear-sky conditions adjusted for topography derived using the toporad algorithm (Dozier 1980; Dubayah et al. 1990), as presented in Image Processing Workbench (IPW) software (Frew and Dozier 1986; Marks et al. 1999). In order to gauge relative radiation inputs to the snowpack, solar radiation was modeled on the 15th day of each month with snow cover prior to the survey date (November–May) and summed to formulate a seasonal radiation index.
a. Redistribution parameters
The terrain-based redistribution parameters presented here quantify different aspects of the upwind terrain. Hence, the upwind terrain needed to be defined. In this study, the upwind terrain consisted of a pie-shaped area centered on the seasonal prevailing wind direction (Fig. 3). Based on the wind data gathered at the D-1 meteorological site, an upwind window of analysis centered on a prevailing wind direction of 265° was chosen. Extending upwind of each pixel, the window was laterally bound by two azimuths separated by an arbitrarily chosen width of 60° (A1 = 235°, A2 = 295°). The lateral width of the window allows for inclusion of terrain that is not directly upwind of the cell based on the D-1 observations. While this approach is somewhat robust for wind directions that vary from the seasonal mean at D-1, it still contains a natural bias to that mean. Calculated shelter/exposure indices at sites where the prevailing wind direction was systematically different (e.g., due to topographic deflection) may contain some error. Specification of a search distance (dmax) controlled the extent of upwind terrain analyzed and was varied to assess the effects of near and far terrain upon snow accumulation.
The objective of the maximum upwind slope (Sx) parameter was to quantify the extent of shelter or exposure provided by the terrain upwind of each pixel. To determine the value of the Sx parameter, we examined all cells along a vector emanating from the cell of interest, determined the cell that had the greatest upward slope relative to the cell of interest, and returned the slope between this shelter-defining cell and the cell of interest:
where A is the azimuth of the search direction, (xi, yi) are the coordinates of the cell of interest, and (xυ, yυ) are the set of all cell coordinates located along the search vector defined by (xi, yi), A, and dmax. The selection of a maximum shelter-producing pixel based on slope is analogous to the determination of solar shading within the horizon function used in radiation modeling (Dozier et al. 1981). Negative Sx values indicated exposure relative to the shelter-defining pixel (i.e., the cell of interest was higher than the shelter-defining pixel). Examples of Sx determinations are provided in Fig. 4. Theoretically, increasingly negative Sx values correspond to greater constrictions on the approaching wind flow, leading to increasing wind speeds, while increasingly positive Sx values correspond to larger and more proximal landscape obstacles, producing a greater degree of shelter and lower wind speeds. Values of Sx were determined along vectors at 5° increments within the upwind window and averaged to formulate the mean maximum upwind slope parameter, Sx, defined by
where nυ is the number of search vectors in the window defined by A1, A2, and dmax. Search distances of 50 m (Sx50), 100 m (Sx100), 300 m (Sx300), 500 m (Sx500), 1000 m (Sx1000), and 2000 m (Sx2000) were used to test the influence of near and far terrain.
It was felt that Sx calculated for within-basin sites would have been only slightly affected by the 10-m DEM's restricted view of the western topography, given the considerable western obstruction presented by the Continental Divide. Hence, Sx was calculated from the 10-m DEM.
The Sb parameter was designed to delineate topographic features capable of generating flow separation. To obtain Sb, independent representations of near and far topographic features were necessary. The difference between relative slope measures of the local and outlying topography supplied a measure of the break in upwind slope. To determine Sb, two applications of the Sx algorithm were used to determine a local Sx (Sx1) and an outlying Sx (Sxo) A user-defined separation distance (sepdist) determined the cutoff between the two regional Sx calculations, with Sx was determined along an upwind search vector with a maximum search distance (dmax) equal to the separation distance (sepdist) to determine Sx1, defined as
The distant or outlaying terrain described by Sx0 was defined as the SxA,1000 of the cell located the separation distance upwind of the cell of interest along each respective search vector:
where [(xo − xi)2 + (yo − yi)2]0.5 ≅ sepdist and (xo, yo) ∈ (xυ, yυ). Subtraction of Sxo from Sx1 yielded Sb. The average of each vector calculation of Sb within the upwind window formed a mean measure of upwind slope break, Sb, defined by
In order to determine Sxo for many pixels located near the Continental Divide it was necessary to analyze terrain beyond the western extent of the 10-m DEM. Therefore, Sxo was determined from the 30-m DEM, while calculation of Sx1 remained at the 10-m scale to retain high-resolution depictions of proximal topography.
In addition, Sxo was used to measure the relative upwind exposure of each potential slope break. With Sb providing a means of locating cells downwind of considerable slope breaks, Sxo [calculated from Eq. (2) with Sxo substituted for Sx] was independently evaluated to determine whether the described slope break was subject to the high winds and assumed correspondent high snow-transport fluxes needed to produce substantial downwind drifts. Based on threshold values of Sb and Sxo, a drift delineator variable (D0) was established whereby each sample was classified as either being within or outside of a modeled lee-slope drift zone. Applied in this manner, D0 was a first approximation that broadly defined drift areas; intradrift accumulation differences were not explicitly represented.
b. Data analysis
We sought to quantify how these redistribution parameters (Sx and D0) related to measured snow depths in two ways. First, we sought to answer the question, “Are Sx and D0 significant predictors of snow depth?” while also examining the best length scales to use for calculating Sx and D0. Linear regression and analysis of variance (ANOVA) were used to perform the aforementioned analyses. Second, we asked the question, “What effects do the addition of Sx and D0 have on a multivariate model of snow distribution?” In response to this latter question, regression tree modeling was used to formulate and compare two tree models of snow distribution: one that incorporated the new parameters with solar radiation, elevation, and local slope data as predictors, and one based solely on solar radiation, elevation, and slope.
The first analysis step was to assess the performance of the redistribution parameters as snow depth predictors. The inspection of the mean maximum upwind slope (Sx) parameter included analysis at varying spatial scales. Standard linear analysis techniques were used to determine the strength of the relationships between the parameters and the sampled snow depths. Though linear models are oftentimes inappropriate for modeling snow distribution (Elder et al. 1995), they do provide a means of comparing the relative strengths of predictor–response relationships. Comparisons were also drawn between the variance in depth explained by Sx and that accounted for by elevation, net potential radiation, and slope to test for the relative importance of each of these variables in explaining the observed snow distribution. Raster grids were then produced to visually compare the distribution of Sx values with the distribution of sampled depths. The drift delineator (D0) was applied as a binary factor or dummy variable alongside the Sx wind scalar parameter within a general linear model. F tests on the drift delineator variable would indicate whether observations with similar upwind shelter contained significantly more snow if they were in a modeled drift zone.
The second analysis step was to assess the effect that the addition of the redistribution parameters had on a model of snow distribution. Based on the aforementioned analyses, redistribution parameters were chosen for inclusion as predictor variables along with the commonly used independent variables of snow distribution modeling—elevation, radiation input, and slope—to form a “redistribution” regression tree model (Breiman et al. 1984) of the May 1999 snow distribution. A “nonredistribution” regression tree model, based on just the latter three predictors was also fit to the sample dataset. Quantitative and qualitative comparisons between the regression tree models trained on the two sets of explanatory variables were then conducted. Regression trees have been successfully used to model snow distribution (Elder et al. 1995, 1998; Balk and Elder 2000) and to determine subpixel resolution snow cover (Rosenthal and Dozier 1996).
The fitting of a regression tree model is a two-step process of growing an intentionally overfit tree with subsequent pruning to a statistically defensible size. In this application, the dependent variable was sampled snow depth, and the set of independent variables consisted of the coregistered predictor values at each sample point. The model is grown using a method of binary recursive partitioning whereby a dataset is successively split, based on values of the predictor variables, into increasingly homogeneous subsets of the dependent variable (Clark and Pregibon 1992). These subsets are termed nodes. At each split the mean of the sampled depths within the subset was the value assigned to that node. As the number of subsets (i.e., terminal nodes or tree size) increases, model R2 values correspondingly increase. Plots of R2 versus tree size provided the means to assess the comparative accuracies of the redistribution and nonredistribution models on a node-by-node basis.
Cross validation of tree models is one recommended method of determining the size of tree with the highest defensible predictive accuracy (Clark and Pregibon 1992; Venables and Ripley 1994). In this research, each snow depth sample, together with its coregistered explanatory variables, was randomly split into 10 mutually exclusive subsets for the cross-validation procedure. Nine of the subsets were combined to grow a tree and then the 10th subset was run through the tree with the model deviance recorded at each step through the tree. Cross-validation results for each of the 10 possible ways of holding a set out were accumulated and averaged, producing a tenfold cross validation (Clark and Pregibon 1992). A common feature of the cross-validation plot (model deviance versus tree size) is a lowering of model deviance (i.e., rising R2) with increasing tree size until a flattening out point is reached, followed by a subsequent increase in model deviance (i.e., decreasing R2) resulting from differences between the full dataset and the subsets used in the cross-validation procedure. The flat minimum of the model deviance plot is a suggested range for the optimal amount of terminal nodes in the final pruned tree (Clark and Pregibon 1992). Since the results of a single tenfold cross-validation run have a random element dependent on the initial splitting of the training data, successive cross-validation runs can produce inconsistent results (Venables and Ripley 1994). For a more robust estimation, the results of 100 tenfold cross-validation runs were averaged to produce the results reported here. Analysis of the cross-validation plots for the two models assessed the respective stability of each model.
The regression trees also provided the means to spatially distribute depth values throughout the basin. The spatial representations illustrated the effect that the redistribution parameters had on extrapolation of the information contained in the sampled data to the population of grid cells. Grids representing the spatial distribution of each of the predictor variables were passed through the respective tree models with the depth assigned to each cell determined by the value within the terminal node it was routed to. Regression trees essentially split the data into discrete classes of snow depth, the number of classes dependent on the amount of terminal nodes in the model, and the values contained within them. As applied in this study, the inherent structure of the regression trees served to substantially increase model scale relative to both the grid and process scales.
4. Results and discussion
Sampled snow depths ranged from 0.0 to 11.50 m, with a mean of 2.27 m and a median value of 1.65 m. The distribution of sampled snow depths was right skewed (skewness G = 1.35), which led to violations of the assumption of residual normality in linear analyses assessing the strength of the relationship between predictors and sampled snow depth. A square root transformation of snow depth produced normalized residuals within the general linear models. In the following text, snow depth should always read “square root of snow depth” when reference is made to the linear model fits.
While samples often exhibited a high variance from point to point, general trends of high and low accumulations were apparent at the basin scale (see Fig. 7). As a measure of spatial autocorrelation, the Moran's I statistic ranges from zero to one, with zero indicating no spatial autocorrelation and one signifying maximum spatial autocorrelation. Moran's I calculated for the 504 depth samples was 0.09 (p = 0), which though relatively low was indicative of significant spatial autocorrelation. The nature of alpine snow distribution—highly variable over the plot scale yet exhibiting some order at the basin scale related to drift, scour, elevation, vegetative cover, and radiation inputs—makes it exceedingly difficult, if not impossible, to collect a representative sample set that does not include some degree of spatial autocorrelation. In order to capture and model both the small-scale and larger, basin-scale variability of snow distribution, datasets such as the one applied in this study are desirable. Removal of the spatial autocorrelation from this intensive dataset would entail substantially decreasing the sample size, with a subsequent loss of information on small-scale variability. In the following analysis, the full dataset has therefore been left intact. Due to the remnant spatial autocorrelation there is a degree of uncertainty related to the number of independent samples cited in the foregoing statistical tests of significance.
a. Redistribution parameters
All of the tested Sx parameters were significant predictors of snow depth (Table 1). The seemingly low, yet highly significant, linear R2 values are reflective of the potential for complex nonlinear interaction of the controlling processes in high-elevation basins (Elder et al. 1991) and may also be due in part to scaling effects caused by differences in the model and process scales (Blöschl 1999). These factors were accentuated in this rugged, steep basin where sampled snow depths were observed to vary by an order of magnitude over very small distances.
The level of significance of the Sx parameter slightly increased as dmax increased from 50 to 100 m and then slowly decreased as dmax increased beyond 100 m. Based on maximum upwind slopes, the Sx parameter was only affected by the inclusion of additional upwind terrain if the distant terrain provided a higher degree of shading than the near terrain previously analyzed. Therefore, as dmax increased, a sample's Sx value could either increase or remain unchanged. For example, a sample located on a perfectly flat area with an abrupt obstacle located 1000 m upwind would have Sx values of 0 for all Sx with dmax < 1000 m, and Sx values greater than 0 for dmax ≥ 1000 m. Analysis of the coregistered Sx100 data indicated that the samples most affected by increases in dmax were those with low Sx100 values, so that areas with low proximal shelter tended to be classified with greater relative shelter as the extent of upwind terrain analyzed increased (Fig. 5). The increase in the level of significance of Sx as dmax increased from 50 to 100 m suggested that, based on this dataset, the sheltering effects of upwind features in this region provided additional information on snow accumulation. The decrease in significance levels as dmax increased beyond 100 m suggested that the sheltering effects of terrain in this region may not affect snow accumulation as much as in the near terrain. Based on these results, only the Sx100 parameter is applied in the succeeding analyses.
As a function of Sx100 the data were arbitrarily classified into exposed (Sx100 < 0°), slightly sheltered (0° < Sx100 < 15°), moderately sheltered (15° < Sx100 < 30°), and extremely sheltered samples (Sx100 > 30°) to observe the ability of Sx100 to define zones of varying accumulation. Box plots of depth for each of the Sx100 groups are shown in Fig. 6 and distinctly indicate the ability of the Sx parameter to segregate regions of high and low accumulations.
Comparisons were also drawn between the amount of observed variance in depth explained by Sx100 and that explained by elevation, net potential radiation, and slope. In linear models of snow depth, Sx100 (p < 0.0001, R2 = 0.19) had a substantially stronger relationship to the dependent variable than elevation (p = 0.88), net potential radiation (p = 0.38), and slope (p = 0.0002, R2 = 0.03). The implication was that Sx100 characterized a significant snow distribution process independent of the processes characterized by the other parameters.
Though elevation and radiation input have been shown to be important controlling processes on snow distribution (e.g., U.S. Army Corp of Engineers 1956; Meiman 1970; Elder et al. 1995), Balk et al. (1998) similarly found low levels of significance for elevation and radiation as predictors of snow depth in another Colorado Rockies Front Range headwater basin. While issues of scale and leeward locations may be affecting the significance of elevation as a factor in these Front Range basins, the low significance of radiation input, which exhibits tremendous heterogeneity over the study areas, highlights how complex process interactions can diminish measures of statistical significance taken from individual parameter tests.
The grid of Sx100 (Fig. 7) highlights the areas of the basin located downwind of sheltering topography and those areas that lack a proximal sheltering feature. These highlights are particularly evident as contour lines are traced beneath the watershed-defining northern and southern ridges and in the relatively flat area approximately 350 m east of the Continental Divide. The distribution of Sx100 appears to closely follow the pattern of sampled depths, particularly in the mid- to lower basin, where the dominant wind pattern may be more closely associated with that observed at the D-1 meteorological station. Due to complex wind patterns in the upper basin, the D-1 data may not be representative of localized conditions and may be the source of some error.
Drift zones were modeled using a separation distance of 300 m and threshold values for Sb and Sxo of 7° and 5°, respectively. Detailed analysis of the Sb parameter using an earlier dataset that encompassed the current study area indicated that increasing the separation distance effectively increased the extent of downwind terrain affected by each topographic break, and that a 300-m separation distance was a best fit to conditions in this basin (Winstral 1999). The best fit of the 300-m separation distance to these data may reflect coincidental downwind and downslope expansion of the accumulation zones due to avalanching and sloughing. Additional data would be required to separate these redeposition effects, while application of Sb to other basins would require calibrating sepdist to local conditions. Since grid-determined slopes and the process of deriving averaged parameters across an upwind window both tend to smooth outstanding features, Sb-calculated slopes would be less than those obtained from field measurements. Guided by field studies citing minimum requirements for flow separation of 7.5°–12° (Berg 1977; Tabler 1994), the 7° threshold was chosen to determine lee-slope drift zones. The threshold for Sxo applied here states that the Sb-defined upwind slope break must have an average maximum exposure of less than 5° to be considered sufficiently susceptible to the large transport fluxes necessary to produce significant downwind drifts.
The general linear model
was fit to the dataset. Results indicated that the p values for the coefficients associated with both Sx and D0 from the fitting of Eq. (10) to these data were significant (Table 2). Given similar conditions of near shelter (100 m), points located within defined drift zones had significantly higher observed snow depths than points located outside of defined flow separation zones.
The areas of the basin modeled as drift zones can be seen in Fig. 7. The largest modeled drift zone in the lee of the Continental Divide corresponds well with the deeper depth samples taken in the central and southern regions of this zone. The width of the northern extension of this modeled drift zone is, however, questionable. It is possible that winds are being funneled down the Continental Divide, producing a more northerly flow pattern at Navajo Peak not evident in the wind data garnered from the D-1 meteorological station. Another possibility is that lateral convergence of flow, not accounted for in the presented modeling techniques, may be occurring downwind of Navajo Peak, producing higher wind gusts in this area. From the available data, some questions also arise as to the accuracy of the smaller drift zones modeled to the south of the spine between Navajo Peak and the D-1 observation site. Though there were insufficient data available to verify the westernmost of these drift zones, the pattern of deeper samples downslope of it suggest that an accumulation zone may exist upslope of the samples. The other drift zones along this spine have more southerly aspects, and their delineation as drifts is very sensitive to any southerly shift in the prevailing wind direction. The use of a static flow pattern as employed in this study does not account for the natural variance of winds seen through the season. These sites may in fact be subject to alternating periods of deposition and scour based on changes in wind flow, and perhaps a modeling approach that accounts for wind dynamics at a smaller timescale would be more appropriate (see Marks and Winstral 2001).
b. Snow distribution models
Two regression trees were grown to distribute snow depth. In the redistribution tree, the set of predictors contained the variables Sx100, D0, elevation, net potential radiation, and slope, while the nonredistribution tree was based on the latter three predictors. Standard diagnostic checking did not reveal any gross heteroscedasticity or non-normality of the residuals in the regression tree modeling that would have necessitated transformation of the dependent variable, snow depth. The plot of model R2 versus tree size (Fig. 8) indicated that the redistribution tree was considerably more accurate than the nonredistribution tree. On a node-by-node basis, the increase in accuracy was approximately 12 percentage points through tree sizes of 30. The cross-validation tests (Fig. 9) for the redistribution tree produced a distinct minimum, with a suggested tree size of 16 (R2 = 0.50). Cross-validation results for the nonredistribution tree, however, were inconclusive, producing a range of optimal tree sizes. The nonredistribution tree's cross-validation plot had a minimum at 9 nodes (R2 = 0.27) and exhibited slightly increasing deviances as tree size grew to 20 nodes (R2 = 0.42), before rising at a rate clearly indicative of model overfit. Based on these results, dependent on the selection of an optimal tree size for the nonredistribution tree, the redistribution tree explained 8%–23% more variance in the observed snow distribution than the tree that did not include a redistribution parameter.
The redistribution regression tree with 16 terminal nodes is shown in Fig. 10. The first three splits of the tree, all based on the Sx100 variable, divided the dataset into four groups of increasing depth correlated with increasing wind shelter. As the tree evolved, splitting decisions were increasingly based on elevation, radiation input, and the drift delineator (D0). None of the splits were based on slope. It was encouraging that all of the splits based on the redistribution parameters were physically meaningful (i.e., increasing shelter contributed to greater modeled depth). The structure of the tree is indicative of the Sx100 parameter, explaining the first-order variability of the observations, with elevation, radiation input, and the drift delineator (D0) having a greater influence on the second-order variability. Evaluation of D0 in this regard is, however, biased. Though regression trees are capable of handling both continuous and categorical data in the same tree, continuous variables have a distinct competitive advantage over binary classified variables such as D0 in the determination of node splits (Winstral 1999). The D0-delineated drift samples largely followed the splits based on increasing values of Sx100. Of the original 504 samples, 56 (11%) met the D0 criteria for drift classification, while, after the first tree split, the 3.75-m node contained 121 sample points, of which 34 (28%) were D0-drift samples, with the following terminal node of 6.03 m containing 11 points, of which 6 (55%) were D0-drift samples. Interestingly, while overall results indicated the need for the slope-break exposure qualification (Sxo) on drift delineation, the five data points at the 6.03-m node that failed to meet the D0 requirements had sufficient upwind slope breaks yet failed to meet the exposure requirement. Further, more detailed analyses of delineating drifts using the available parameters appears to be warranted.
Extensive snow surveys, like the one applied in this study, are most often limited to small basins, due to cost, accessibility, and manpower concerns. This limitation necessitates some means of extrapolating these basin-scale findings for regional applications. Examination of the regression tree images of distributed snow depth (Fig. 11) offered insight into model performance outside the immediately sampled area. Scaling issues between the broader tree-based model scale and the finer scales of the snow distribution process and modeling of Sx are evident in comparisons with Fig. 7. Both regression tree images were based on models with 16 terminal nodes, with slopes greater than 50° assumed to be too steep to accumulate snow and assigned depth values of zero. The most striking discrepancies between the redistribution and nonredistribution trees occurred along the north-northwest-facing ridge to the west of Kiowa Peak. In this high-elevation region with high wind exposure and low radiation inputs, the nonredistribution tree, lacking a direct means of accounting for wind exposure, assigned the greatest depths in the basin, 4–8 m. The redistribution tree on the other hand, designated this region as a low-accumulation zone having a predominance of 1–2-m depths, with slightly greater depths behind northerly trending transverse ridges. Visual observations made during the study period indicated that this area is in fact a region of very low snow accumulation. Aerial assessments of snow-covered area were not available for the period of the 1999 snow survey; however, there was a classified photograph from 22 May 1996 [1 May snowpack 167% of average (NRCS 1996)] available to compare with the May 1999 [1 May snowpack 112% of average (NRCS 1999)] modeled snow distributions (Fig. 12). Assuming that the locations of low- and high-accumulation zones are annually consistent regardless of snowfall totals, the classified photograph further supported the redistribution model.
Due to the spatial correlation issues mentioned earlier, there was some degree of uncertainty attached to the statistical tests of significance in the first step of data analysis. It is therefore critical to assess whether the parameters are capturing a physical process or picking up idiosyncrasies in either the data or physical characteristics of the basin. Aside from the clearly evident relationship between modeled wind exposure and snow distribution described in the preceding paragraph, there exists other visible correlations between the parameters and snow distribution that cross radiation loadings, elevation bands, and basin characteristics. An example of this is the high snow-accumulation band that extends roughly from the southernmost point in the basin in a northeasterly direction down to the valley floor and part way up the north wall of the basin (Figs. 7 and 12). Radiation inputs along this band range from low along the northeast-facing slopes in the south to high along the southeast-facing slopes in the north. There is an approximate elevation range of 350 m. While radiation loadings, elevation, and slope vary across this region, Sx100values remain consistently high (Fig. 7), suggestive of a wind-sheltered region consistent with the observed higher accumulations. A broad visual analysis of Fig. 12 also corroborates some of the earlier statistical findings—that correlations between solar radiation input, elevation, slope, and snow distribution in this basin are very weak. There are, however, striking similarities between the redistribution parameters (Fig. 7), the redistribution regression tree model (Fig. 11b), and the observed snow patterns (Figs. 7 and 12). In this basin, where wind is widely acknowledged as the dominant control on snow distribution, the terrain-based redistribution parameters have seemingly captured snow distribution patterns unrelated to radiation input, elevation, slope, or any apparent basin physiographic features other than those specifically defined by the parameters.
Digital terrain analysis was employed to develop parameters that characterized the effects of wind on snow distribution. The Sx parameter was designed to quantify the degree of shelter/exposure provided by upwind terrain, thus characterizing the wind scalar. Derivations of Sx that sequentially analyzed the nearest 50, 100, 300, 500, 1000, and 2000 m of terrain were all significant predictors of snow depth. Based on these varying scales of analysis, it was found that the upwind terrain within 100 m (Sx100) of a sample was the strongest predictor of the snow distribution sampled in this basin. Furthermore, it was found that Sx100 explained variance in snow depth that was not accounted for by differences in elevation, net potential radiation, or slope. A drift delineator, D0, was formulated based on a measure of upwind breaks in slope, Sb, and an application of Sx that determined the exposure of the upwind slope break (Sxo). In regions of similar upwind shelter, significantly deeper snow accumulations were found within the modeled drift zones.
A regression tree model of snow distribution that included Sx100, D0, elevation, radiation input, and slope explained 8%–23% more of the observed variance than the tree model based on just the latter three predictors. The model that included the redistribution parameters coincided with the general field observation that snow primarily accumulates in areas sheltered from the consistently strong winds present in this basin, with very little snow accumulating in exposed areas. Acute differences that would have a significant effect upon runoff forecasts based on these models were observed in high-elevation areas that were shaded from the sun yet lay exposed to prevailing winds. Ancillary data supported the redistribution model and lent further credence to the use of the wind redistribution parameterizations.
Many researchers have noted the importance of wind redistribution on snow accumulation patterns in alpine environs, yet the ability to include the effects of this process within spatial snow models has been limited. This analysis demonstrated how just a basic application of the redistribution parameters could substantially improve an empirical model of snow distribution while also supporting previous results obtained in other alpine basins (cf. Winstral et al. 1998; Winstral 1999). All indications have been that the presented procedures are transferable to other basins with predominant storm-flow patterns. We believe the presented findings represent a significant first step in accounting for wind-induced snow redistribution with terrain parameterizations in alpine basins.
From a modeling perspective, the strong point of these redistribution parameters is their computational efficiency, which readily allows for incorporation of additional detail into the modeling procedure. The parameters are ideally suited for distributed modeling. Currently, we are working on linking the terrain parameters to a distributed snow accumulation and melt model in which the redistribution parameters evolve through time relative to the developing snowpack, and snow accumulation is modeled at each time step as a function of snowfall, winds, and the redistribution parameters. We anticipate that additional climatic inputs such as these will lead to the development of a computationally attractive snow redistribution model that is robust to intra- and interannual variability.
Support for this project was provided by the USACE Cold Regions Research and Engineering Laboratory. Logistical support and/or data were provided by Mark Williams, Don Cline, the NSF-supported Niwot Ridge LTER project, and the University of Colorado Mountain Research Station. We would like to express our gratitude to G. Blöschl, C. Luce, D. Marks, D. Tarboton, and two anonymous JHM reviewers for their thoughtful and valuable comments on the manuscript.
Mention of trade names or commercial products does not constitute endorsement or recommendation for use.
Corresponding author address: Adam Winstral, Northwest Watershed Research Center, 800 Park Blvd., Suite 105, Boise, ID 83712. Email: firstname.lastname@example.org