Fractal dimensions derived from log–log variograms are useful for characterizing spatial structure and scaling behavior in snow depth distributions. This study examines the temporal consistency of snow depth scaling features at two sites using snow depth distributions derived from lidar datasets collected in 2003 and 2005. The temporal snow accumulation patterns in these two years were substantially different, but both years represent nearly average 1 April accumulation depths for these sites, with consistent statistical distributions. Two distinct fractal regions are observed in each log–log variogram, separated by a scale break, which indicates a length scale at which a substantial change in the driving processes exists. The lag distance of the scale break is 15 m at the Walton Creek site and 40 m at the Alpine site. The datasets show consistent fractal dimensions and scale break distances between the two years, suggesting that the scaling features observed in spatial snow depth distributions are largely determined by physiography and vegetation characteristics and are relatively insensitive to annual variations in snowfall. Directional variograms also show consistent patterns between years, with smaller fractal dimensions aligned with the dominant wind direction at each site.
Many hydrologic (Luce et al. 1998; Elder et al. 1991), ecologic (Jones 1999; Brooks and Williams 1999), and climatic (Liston 1999; Groisman and Davies 2001) investigations require that spatial distributions of snow amounts be known or estimated. In complex and mountainous terrain, precipitation and snow redistribution amounts and patterns are often highly variable and difficult to measure. The complex interactions of snowfall and wind, with terrain and vegetation, introduce formidable sampling and modeling problems for characterizing spatial snow properties (Elder et al. 1991; Winstral et al. 2002; Erickson et al. 2005). However, because seasonal snow and its spatial distribution at many scales is dynamically linked to hydrologic, atmospheric, and biologic systems through forcing of runoff characteristics, heat and energy fluxes (Liston 1999), soil moisture distributions (Sturm et al. 2001), and growing season duration (Jones et al. 2001), the ability to measure and model snow distributions is critical to understanding and representing the processes governing energy, water, and biogeochemical cycling in mountain and earth surface systems.
The nature of the process interactions creating spatial snow distributions is complex, and the observed variability changes with the scale of observation (Blöschl 1999). Thus, scale is central to any assessment of the spatial distribution of snow. Recent work has shown that spatial snow distributions in a variety of environments exhibit fractal characteristics (Shook et al. 1993; Shook and Gray 1996, 1997; Kuchment and Gelfan 1997; Granger et al. 2002; Litaor et al. 2002; Deems et al. 2006). Fractal distributions indicate consistent driving process relationships over the range of applicable scales and provide a theoretical basis for sampling, modeling, and rescaling spatial snow data, and for understanding the underlying process interactions.
Deems et al. (2006) examine snow depth, topography, and vegetation topography for fractal scaling characteristics at three locations. Snow depth at all sites is found to show two distinct regions of fractal scaling separated by a scale break where process dynamics appear to change character. The lag distance of the scale break varies among the sites but is consistently longer where the overall relief of the study area is higher, as found by Shook and Gray (1996). The fractal dimensions of each scale range are very consistent between sites, suggesting that the scaling properties in the snow depth distributions are consistent among different physiographic and vegetation covers within the same snow climate. Additionally, the fractal distributions are shown to be anisotropic, with directional patterns related to the dominant wind direction at each site.
In this study, building on the work presented in Deems et al. (2006), a second lidar data acquisition allows an investigation of interannual consistency into the observed fractal distributions and scaling features at two of the sites used in the prior study. Consistency between different snow seasons with differing accumulation characteristics would indicate that the scaling characteristics are intrinsic to the specific site and are relatively insensitive to variations in weather patterns. A robust estimation of scaling characteristics for a single site would allow for more efficient and accurate sampling design, data interpolation, and rescaling of existing data.
2. Study areas and snow seasons
a. Study areas
Two sites from the National Aeronautics and Space Administration’s (NASA) Cold Land Processes Experiment (CLPX) in Colorado were used in this study: the Walton Creek and Alpine intensive study areas (ISAs; Cline et al. 2001; Figs. 1 and 2). The moderate-elevation Walton Creek site, in the Yampa River Basin, is characterized by a broad meadow environment interspersed with small, dense stands of coniferous forest; low rolling topography; and deep snowpacks. The high-elevation Alpine study site, in the Upper Colorado River basin, has alpine tundra with some subalpine coniferous forest and is generally north facing, with moderate relief (Cline et al. 2003). Walton Creek receives high annual snowfall and calm-to-moderate winds, whereas the Alpine site receives lower annual precipitation totals and has greater wind exposure above the treeline.
The two study areas provide a contrast of terrain and vegetation environments, from low to high relief and from alpine tundra to subalpine forest of varying density. Field observations indicate that snow distributions at both sites are dominated by wind redistribution and interaction, with terrain and vegetation patterns. Wind direction frequency distributions, calculated from meteorological data collected at each of the study sites for both snow seasons in this study (Elder et al. 2008), indicate that winds of speeds greater than 5 m s−1 (at 10 m elevation) are confined to narrow direction bands (Fig. 3).
b. 2003 and 2005 snow seasons
The snow season evolution for 2003 and 2005 can be characterized through micrometeorological data collected at each CLPX study site as well as data from nearby U.S. Department of Agriculture (USDA) Natural Resource Conservation Service snowpack telemetry (SNOTEL) sites. These data allow an interannual comparison of precipitation amounts and wind redistribution, two major controls on spatial snow accumulation patterns at slope-to-basin scales.
The Yampa River basin as a whole experienced very similar seasonal total snow water equivalent (SWE) and precipitation accumulation in 2003 and 2005, as did the Upper Colorado basin. As of 1 April in both 2003 and 2005, the Yampa River basin was at 90% and 89% of average SWE, respectively, whereas the Upper Colorado River basin was at 98% of average in both years. Total precipitation for the Yampa River basin (Upper Colorado basin) was reported as 92% (97%) of average in 2003 and 88% (94%) in 2005. The basin-wide percent of average values show that at the scale of synoptic weather patterns, both 2003 and 2005 snow years were near the 30-yr average. This scale of comparison gives a qualitative indication that the storm patterns in the region, which exert a control on the dominant wind direction during snowfall events, were substantially similar in both years.
To compare the interannual variation in snow accumulation at nearby SNOTEL stations, we examined daily SWE totals because a snow depth sensor was not installed at the Rabbit Ears SNOTEL until 2006. As with the basin SWE summary data, 1 April SWE values from the individual SNOTEL stations show similar values for both years, though with greater variation than the basins as a whole. Relative to the entire basins, the Rabbit Ears and Berthoud Pass Summit SNOTEL sites showed higher SWE accumulations in 2003 and lower SWE in 2005 relative to the 30-yr mean for each site. As might be expected, there is more year-to-year variation in individual station data than in the basin-wide aggregate data.
Figure 4 shows daily SWE observations from the nearby SNOTEL sites for both 2003 and 2005. The curves for both years are similar at both sites, with lower SWE totals in 2005 and complete depletion by early to mid-June. The difference in SWE amounts between years at Berthoud Pass Summit appears to be mainly a result of a precipitation event in mid-March 2003—a large easterly, upslope storm that did not reach far enough to the west to affect the Rabbit Ears site.
Individual storms assume a greater importance in areas that routinely experience shallow snowpacks, like the above-treeline region at the Alpine site, because a single storm can multiply the existing snow totals several fold (Serreze et al. 2001). In deeper snowpack environments, such as the Walton Creek site, a single storm snowfall total is usually a much smaller percentage of the snowpack (Fig. 5a). This factor is critical when examining the measurements taken in two different years because a change in measurement date by only a few days can dramatically affect the observed distribution. This fact is well demonstrated by the snow depth time series from the Alpine study site (Fig. 5b). Both years show an extremely temporally variable snowpack, though for any given day of the water year, the 2005 snowpack tended to be deeper. Notably, both series show a large storm occurring late in the season, which significantly changed the snow depths for several days. However, each lidar dataset was acquired several days poststorm when snow depths were much lower and of similar magnitude. Therefore, though the snow distributions at the Alpine site appear similar in magnitude between years, different measurement dates may have produced very different results. Measured distributions at the Walton Creek site, by contrast, are not as sensitive to acquisition date.
The other major weather variable that drives snow distributions is wind. Wind speeds above a transport threshold value can redistribute significant volumes of snow. Figures 3a–d show the frequency distributions of wind direction for both years at both sites for wind speeds greater than 5 m s−1 and for wind during or just after precipitation events. The distributions show little difference in the direction of peak wind frequency for the two years at each site. The 2005 data show an increase in the frequency of southwesterly winds associated with precipitation events at both sites, making the distributions slightly bimodal. However, the dominant precipitation wind direction is consistent between the years. Overall, the two snow depth datasets obtained for each site represent spatial snow distributions near maximum accumulation in the two very similar snow years.
a. Lidar altimetry
Lidar altimetry is used for many applications, including terrain modeling and forest structure mapping. Lidar altimetry offers the ability to process multiple laser returns to get bare ground or snow, and canopy height information. The combination of high accuracy and areal coverage makes airborne lidar datasets appealing for mountain hydrology studies. The use of lidar is becoming more common for volume estimation and change detection, including snow depth retrieval (Hopkinson et al. 2004; Deems et al. 2006) and glacier mass balance estimation (e.g., Hopkinson and Demuth 2006). Wehr and Lohr (1999) and Baltsavias (1999) provide a comprehensive overview of airborne lidar systems and data collection procedures.
Accuracy of airborne lidar data is a function of system and flight parameters and ground surface conditions [see Baltsavias (1999) and Hodgson and Bresnahan (2003) for comprehensive reviews]. Errors due to global positioning and inertial navigation systems are assumed to be random and typically contribute on the order of 10 cm to the range error (Baltsavias 1999). Flight parameters such as altitude, flight line orientation, swath overlap, scan angle, pulse rate, and beam divergence control the spacing of return pulse locations and, therefore, the accuracy of the filtered data product.
Vegetation cover and terrain geometry (mostly slope) influence the strength of ground returns and vertical accuracy, respectively. Steep slopes can reduce accuracy through “time walk” and horizontal error (Baltsavias 1999; Hodgson and Bresnahan, 2004). Time walk refers to the increased spread of laser pulse energy on a slope and the corresponding delay in rise time, or the time required for return energy to exceed the recording threshold. The magnitude of time walk errors depends on the beam divergence and pulse length parameters. Vertical error can also be induced by errors in horizontal position, whereby an otherwise accurate range measurement subject to (x, y) positional error has an apparent position above or below the ground surface. Because the horizontal error distribution is assumed to be random in magnitude and direction, the vertical error induced by this mechanism will vary randomly from 0 (for horizontal errors of zero or perpendicular to the slope direction) to a maximum determined by
where ɛz and ɛx,y are the vertical and horizontal error, respectively, and α is the terrain slope (Hodgson and Bresnahan 2004). Error contribution from the terrain slope can be potentially large though in practice, extra swath overlap is added or flight line orientation is changed to increase laser shot density in areas of complex terrain. A larger number of laser shots provides more candidate ground returns from which the filtering algorithm can estimate the ground surface location. Sufficiently large slope-induced errors would be subject to removal by filtering algorithms.
Three lidar datasets from each site were used for this study; they were acquired on 9 April 2003, 19 September 2003, and 1 April 2005. The raw data were normalized using ground control points, postprocessed to remove redundant data points and noise, and then classified as ground or vegetation points by the lidar contractor using a proprietary morphological filter, producing a dataset with 1.5-m nominal horizontal spacing between data points and 0.05-m vertical measurement tolerance according to contractor specifications (Miller 2003). The last-return signal from the September 2003 mission provides ground surface (“bare earth”) elevations, the September first-return data measures the terrain-plus-vegetation (vegetation topography) elevations, and the April 2003 and 2005 last-return data provide snow surface elevations. A 1-m resolution digital elevation model (DEM) was produced from the bare earth point data using inverse–distance–weighting interpolation. The DEM elevations were then subtracted from the snow surface elevation points, producing datasets of snow depth point estimates. The point datasets of snow depth, bare earth elevations, and vegetation topography were used for the variogram fractal analysis.
Figure 6 shows the histograms for the derived snow depth point datasets. The deep snow environment at Walton Creek produces a relatively normal distribution of snow depths, with slightly deeper overall snow accumulation in 2003. The Alpine site shows a much higher proportion of low and zero snow depths, resulting in a more complex distribution, similar to the mixed lognormal depth distributions reported by Marchand and Killingtveit (2004). Of interest are the substantial number of zero and negative snow depths calculated for the Alpine site. Some negative depths are to be expected because of random measurement errors; however, some unexpectedly large negative depths appear in the 2003 Alpine dataset. Though the source of these erroneous numbers is unknown, they are possibly the result of a georegistration or filtering error on the part of the lidar contractor. A visual comparison of the location of the negative values indicates that the vast majority are located in snow-free areas above treeline, as seen in the orthophotos (Fig. 2b). Snow depth maps (Fig. 3) show that the spatial depth patterns are realistic and exhibit textures expected of forested and tundra areas. Therefore, we decided to set the negative values to zero for this analysis. Table 1 contains summary statistics for the snow depth data used.
b. Variogram fractal analysis
Fractal dimensions are calculated using variograms, which show the amount of variance between samples as a function of their distance scale of separation. As referenced in Deems et al (2006) omnidirectional and directional semivariograms, γ(r), are estimated using 50 log-width bins:
where r is the lag distance of bin k, N is the total number of pairs of points in the kth bin, and zi and zj are the snow depth values at two different point locations i and j (Webster and Oliver 2001). Log-width distance bins are used to provide equal bin widths when the variograms are transformed to log–log space. Log-width bins also allow for a greater bin density at short lag distances than would linear-width bins, which aids in resolving the variogram structure at short length scales; therefore log-width bins allow a more precise power-law fit than would standard linear distance bins (Klinkenberg and Goodchild 1992). The variograms were restricted to a maximum distance of 1100 m, the diameter of the largest circle that can be fit within the mapped area, which helps avoid directional bias (Mark and Aronson 1984). Additionally, any nonstationarity in the data was not removed, though that is standard practice in a geostatistical analysis. In fractal analysis, the fractal dimension is an index of the relative balance of long- and short-range processes, and removing large-scale trends would bias the calculated fractal dimensions toward short-range variability (Klinkenberg and Goodchild 1992).
The variograms are log-transformed to allow identification of regions that can be described by a power law. Linear regions in each log–log variogram are identified visually, and each straight section is fit with a power function of the form
by varying coefficients a and b to minimize the squared residuals.
The fractal dimension (D) is estimated from the slope (power) of the log–log variogram by Gao and Xia (1996)
When breaks in slope are observed in the log variogram plots, the scale break length is determined from the intersection of the two fitted power-law curves. It is recognized that, rather than a discrete break, the slope changes continuously from one linear segment to another. It should be noted, however, that this change in slope occurs over a relatively short scale range. The power-law relationships display very good fits with the transition region included in the power-law segments (R2 based on the log data greater than 0.90). Solving the two power-law equations to derive a single break point provides a consistent measure that is comparable between datasets. Shook and Gray (1996) used a similar methodology, though they assumed that the longer-range sections approached a completely random spatial distribution (flat slope, D = 3). Their technique results in slightly larger estimates of the scale break distance than the technique used here. Our methods likely produce a slightly smaller estimate of the scale break location.
The 2003 and 2005 snow depth datasets were compared to examine the consistency of the spatial snow depth distributions from year to year. Because of a problem with aircraft equipment during the 2005 mission, full coverage of the Alpine site was not obtained. Therefore, the 2005 Alpine variograms were compared with 2003 Alpine variograms computed from both the full dataset and a dataset clipped to match the extent of the 2005 data.
4. Results and discussion
a. Omnidirectional variograms
Three variogram features were examined for interannual consistency in the spatial snow depth distributions: overall semivariance, the lag distance of the scale break, and the slopes of the power-law fractal segments. Results are consistent between years at both sites, showing only a change in the overall magnitude of the variability at all lag distances, with essentially no changes in the spatial structure or fractal scaling properties. This interpretation is supported by the snow depth difference maps (Fig. 7, bottom panels). Though the overall patterns and scales of variability appear similar in both years and the actual depth values are comparable, with correlation coefficients of 0.88 for Walton Creek and 0.92 for Alpine, there are substantial differences in depth values between years at both sites. The difference maps highlight changes at the scale of individual grid cells, but the fractal analysis helps to quantify the similarity in the overall pattern by identifying scale features that are difficult to differentiate visually.
The Walton Creek omnidirectional variograms are nearly identical between seasons (Fig. 8a). The scale break lag distance and the fractal slopes are the same, indicating that the spatial distribution of snow depths had the same structure in both years, with virtually identical scaling properties and relative amounts of long- and short-range variation (Table 2). The only difference is a slightly higher overall semivariance in 2005. This difference could be explained by the lower overall snow depth in 2005 (see Figs. 5 and 6), which would have masked less of the terrain and vegetation roughness and produced a larger variance magnitude, though with the same spatial structure. The accumulation and wind redistribution history certainly differ between years, as suggested by the slightly higher proportion of storm winds from the southwest during 2005. The interannual consistency in fractal dimensions and scale break distance demonstrates that the process interactions that create the spatial snow depth pattern are identical in each year. Any differences in the spatial snow depth patterns are not sufficient to disrupt the overall scaling relationships created by the snowfall–wind–terrain–vegetation interactions over the 1– 1000-m scale range.
When the Alpine variograms from 2003 and 2005 are compared, the initial result indicates that there is a change in the spatial distribution and fractal dimensions (Fig. 8b). However, because an equipment issue precluded data collection for the entire site in 2005, the spatial extent of the 2003 data must be clipped to match that of the 2005 data for the comparison to be valid. Areas of different size will contain different elevation and vegetation distributions—different distributions of the roughness elements that control the pattern of wind redistribution. Therefore, unless the roughness elements have similar patterns, areas of different extent could be expected to have different spatial snow distributions. Previous studies have shown a qualitative relationship between scale break distance and the overall relief of the study area (Shook and Gray 1996; Deems et al. 2006). The clipped Alpine extent has a smaller overall relief than the full study area (298 versus 339 m) and also a shorter scale break length, which is consistent with the prior findings. Using the same variogram search distance parameters to analyze the clipped area could induce edge effect errors, though these errors are confined to the longest lag bins (as seen by Fassnacht and Deems 2006), appearing as an anomalous change in semivariance (Fig. 8). The longest lags are excluded from the power-law regressions, thereby limiting edge effect influences.
When the 2005 Alpine variogram is compared to the variogram from the clipped 2003 dataset, the distributions are very similar (Fig. 8b; Table 2); both long- and short-range fractal dimensions vary within 0.02, the scale break distance is within 5 m, and the overall variance magnitude is effectively equivalent. Thus, the scaling characteristics of the snow distributions in the two years were consistent on the dates sampled, indicating consistent process relationships despite differences in the snowfall and wind history between the two years. Low snow depths and wind exposure make the Alpine site more sensitive to the date of sampling than the Walton Creek site. It is possible that if either sampling date (2003 or 2005) had been closer to the large spring snowstorm, the resulting spatial scaling patterns might be different. However, based on the snow depth time series, the sampling dates appear to be a good representation of “normal” conditions at the Alpine site.
b. Directional variograms
Snow depth fractal dimensions at the Walton Creek site are virtually identical between years for all directions and for both the long-range (15–550 m) and short-range (1.5–15 m) and fractal segments (Fig. 9a). Two conclusions can be drawn from this comparison. First, there appears to be no difference between years in the anisotropy of the process relationships that determine the snow depth pattern and pattern scaling relationships. This conclusion is supported by the similarity between the wind direction frequencies (Figs. 3a,b). Assuming that any changes in the vegetation (and topography) are insignificant over two years, interannual differences in the wind history will exert a dominant influence on the differences in snow redistribution. Therefore, any differences in the wind direction frequencies, such as the slightly higher proportion of southwest winds associated with the 2005 precipitation events, are not sufficient to substantially alter the overall scaling pattern.
Second, in practical application, any sampling or interpolation scheme designed to capture the anisotropy evident in the short lag distances would be applicable to both snow seasons. Measured or modeled spatial distributions could be rescaled using anisotropic scaling factors, and the results indicate that the same anisotropic rescaling methods might be applicable to the data from both years.
At the Alpine site, the anisotropy is also similar between years, though slightly more variable than at Walton Creek (Fig. 9b). Two directional features are notable. First, short-range fractal dimensions are much higher in the south-southwest–north-northeast direction, which is perpendicular to the dominant wind direction. Second, the short- and long-range fractal dimensions are the same in the south-southwest–north-northeast direction. In other words, there is no scale break perpendicular to the dominant wind direction, and the snow depth distribution scales continuously from 1.5 to about 270 m. In this direction, there is virtually no change in elevation (see Fig. 2b). Because the only variation along an elevation contour is in the vegetation pattern, there is no process change at any lag distance within the range of scales represented in the data, and hence no scale break.
The wind direction frequencies at Alpine are consistent between years (Figs. 3c,d), such that the process relationships that determine snow distributions, vegetation, topography, and wind are equivalent. As with the Walton Creek site, there is a slight increase in the proportion of southwesterly winds associated with precipitation events in 2005, which seems to have had little impact on both the overall and the directional spatial distributions.
It is instructive to consider the time of year that data were collected. Both years were sampled near 1 April, the date often used to represent maximum accumulation. By the latter part of the accumulation season, the storm sequence is likely to have produced a spatial snow distribution that is consistent with the dominant storm track and wind direction. Were this same analysis conducted for 1 January data acquisitions, for example, very different distributions might be expected, with spatial patterns more sensitive to variations in early season snowfall. It is possible, therefore, that late season spatial scaling relationships are more interannually robust than those observed earlier in the year.
The Alpine site warrants additional considerations. Unlike the Walton Creek site, the Alpine site contains some steep slope regions. Lidar accuracy is expected to decrease with slope angle so some error must be assumed to propagate through the analysis. However, the large number of points in the dataset and, consequently the large number of point pairs in the variogram analysis, will tend to reduce the influence of errors from the higher slope regions. Therefore, it is unlikely that the random slope errors induce any substantial bias in the fractal variograms.
Additionally, the Alpine site contains two distinct subregions, with dramatically different accumulation and ablation histories through a single season. The below-treeline region experiences a typical subalpine forest snow season, with relatively well-defined accumulation and ablation periods. The above-treeline alpine tundra region, in contrast, is subject to substantial wind redistribution, with its generally low snow depths fluctuating, and even depleting, multiple times per season. Thus, the interannual comparison of the spatial pattern can be sensitive to the date of observation. Indeed, as the snow depth history from the study site illustrates (Fig. 5b), the time between the last accumulation event and the data collection could affect the observed snow depth distribution in the tundra area.
This sensitivity is best interpreted through the fractal variograms. The scale break is interpreted as a length scale at which snow depth processes change fundamentally in the study domain. The fractal dimension in each scale range separated by the scale break is a function of the process balance over that length range because the process interactions will control the amount of variability present in the spatial snow distributions. From this perspective, changes in snow depths in the tundra areas of the Alpine site should affect the fractal dimension of the two scale ranges by altering the measured semivariance at different lag distances. The scale break distance should remain unchanged because the scale at which the snow depth processes change is driven primarily by the topographic and vegetation distributions within the study area. Future datasets with substantially different snow depths in the Alpine study site will enable further analysis of the effects of depth variation on the fractal dimension.
Snow depth spatial patterns and scaling behavior are compared for 2003 and 2005, using snow depth datasets derived from airborne lidar measurements at two midlatitude, mountainous sites. Both omnidirectional and directional analyses show strong interannual consistency for the two years, with substantially different accumulation histories. The spatial distributions differ slightly in overall variance, but the fractal dimensions over the two distinct scale regions and the length of the scale break that separates the scale regions are nearly identical.
These results demonstrate interannually consistent process relationships among the major driving factors controlling snow accumulation and redistribution, namely, precipitation, wind, vegetation, and topography. The snow seasons that produced the observed snow distributions were similar at the regional, basin, and individual site scales, according to data from nearby SNOTEL sites and micrometeorological data collected onsite. However, considerable differences in the temporal accumulation patterns exist between years. The results suggest that the observed scaling properties of the snow depth distributions are characteristic of the site locations and are relatively insensitive to interannual variation in snow accumulation. The interannual consistency in the scaling characteristics indicates that the fundamental process relationships producing spatial snow depth patterns at these sites are consistent from year to year. The scaling properties of snow distributions produced by substantially deeper or shallower snow years, or at times earlier in the accumulation sequence, remain to be investigated.
This project was supported by NASA headquarters under Earth System Science Fellowship Grant NNG04GQ40H. The lidar and meteorological data were collected as a part of the NASA CLPX. Two anonymous reviews greatly improved this manuscript.
Corresponding author address: Jeffrey S. Deems, Department of Civil and Environmental Engineering, University of Washington, Box 352700, Seattle, WA 98195-4325. Email: email@example.com
This article included in the The Cold Land Processes Experiment (CLPX) special collection.