## Abstract

Snowpack properties vary dramatically over a wide range of spatial scales, from crystal microstructure to regional snow climates. The driving forces of wind, energy balance, and precipitation interact with topography and vegetation to dominate snow depth variability at horizontal scales from 1 to 1000 m. This study uses land surface elevation, vegetation surface elevation, and snow depth data measured using airborne lidar at three sites in north-central Colorado. Fractal dimensions are estimated from the slope of a log-transformed variogram and demonstrate scale-invariant, fractal behavior in the elevation, vegetation, and snow depth datasets. Snow depth and vegetation topography each show two distinct fractal distributions over different scale ranges (multifractal behavior), with short-range fractal dimensions near 2.5 and long-range fractal dimensions around 2.9 at all locations. These fractal ranges are separated by a scale break at 15–40 m, depending on the site, which indicates a process change at that scale. Terrain has a fractal distribution over nearly the entire range of scales available in the data. Directional differences in the fractal dimensions for each parameter are also present at multiple scales, and are related to the wind direction frequency distributions at each site. The results indicate that different sampling resolutions may yield different results and allow rescaling in specific scale ranges. Resolutions of 10 m and finer are consistently self-similar, as are resolutions greater than 30 m, though the coarser resolutions show nearly random distributions.

## 1. Introduction

Knowledge of the spatial distribution of snow in mountain regions is critical for accurate assessment and forecasting of snowmelt timing (Luce et al. 1998), snowmelt volume (Elder et al. 1991), and avalanche hazard (Conway and Abrahamson 1984; Birkeland et al. 1995), for initialization of synoptic- and global-scale weather and climate models (Liston 1999; Groisman and Davies 2001), and for investigations of ecologic dynamics and biogeochemical cycling (Jones 1999; Brooks and Williams 1999). The interactions of snowfall and wind with terrain and vegetation create highly variable patterns of snow accumulation. These complex interactions present formidable sampling and modeling problems for characterizing spatial snow properties (Elder et al. 1991). The seasonal snow system 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, soil moisture distributions, and growing season duration (Jones et al. 2001), greatly influencing processes governing energy, water, and biogeochemical cycling in mountain and earth surface systems.

The nature of the process interactions creating spatial snow distributions are complex and the observed variability changes with the scale of observation (Blöschl 1999); scale is thus central to any assessment of the spatial distribution of snow. Investigations have shown that spatial snow distributions in a variety of environments exhibit fractal characteristics (Shook et al. 1993; Shook and Gray 1996; Shook and Gray 1997; Kuchment and Gelfan 1997; Granger et al. 2002; Litaor et al. 2002). Fractal distributions indicate self-similar properties over multiple scales and may provide a theoretical basis for sampling, modeling, and rescaling spatial snow data and for understanding the underlying process interactions. In this study, snow depth, topography, and vegetation topography data are examined for fractal characteristics and consistencies in scaling characteristics among the three datasets.

## 2. Background

### a. Scale and snow measurement

It is difficult to obtain measurements of snow depth distribution over a large area at a resolution that approximates the true variability present (Blöschl 1999). Manual data collection, via snow stakes or probes, is labor intensive, expensive, and potentially dangerous in steep mountain environments. Airborne and satellite remote sensing methods have been used effectively for estimation of snowpack properties such as snow-covered area and albedo (Dozier and Painter 2004), but are not yet to the point where they provide real-time operational estimates of snow depth or water equivalent on the scale of individual catchments (Marchand and Killingtveit 2004). Recently, the use of high-resolution airborne laser altimetry [light detection and ranging (lidar)] has been explored for fractal analysis of terrain and vegetation patterns (Pachepsky et al. 1997; Pachepsky and Ritchie 1998) and for gathering spatial snow depth data (Hopkinson et al. 2001; Miller et al. 2003). Pachepsky et al. (1997) and Pachepsky and Ritchie (1998) used fractal analysis techniques to study spatial patterns in terrain and vegetation distributions derived from lidar data, but did not examine snow distributions. Hopkinson et al. (2001) and Miller et al. (2003) focused on the utility of lidar altimetry for measuring snow depth, but did not analyze spatial snow distributions. This study builds on the previous work, using fractal analysis to examine spatial patterns and scaling in lidar-derived snow depth distributions.

Snow properties vary over multiple scales, from microscopic crystal morphology and grain bonding to landscape-scale patterns of accumulation and ablation. The amount of variability at each scale is controlled by the process interactions governing snow accumulation, metamorphism, and ablation; as the process interactions change with scale, so too does the variability in the snow cover. As with many earth-surface phenomena, ordered behavior may be found at some scales, while others show incomprehensible complexity (Phillips 1999). Identification of “natural” scales of observation or modeling can lead to better assessment of snow distributions or more effective models. It is important to describe scaling properties of snow depth variation and to identify scales where process interactions change in order to more effectively characterize and model spatial patterns in the snow cover.

Orographic effects on precipitation, dominant at the mountain-range scale, have been observed over elevation changes as small as 50 m (Barry 1992). However, local-scale (1–1000 m) snow depth variability in mountain environments is dominated by wind redistribution of snow (Elder et al. 1991; Blöschl and Kirnbauer 1992; Winstral et al. 2002). Wind redistribution coincident with or subsequent to a precipitation event can quickly alter accumulation patterns. This superposition of accumulation processes leads to a dependence on both orography (terrain) and roughness element (vegetation) distribution (Elder et al. 1991). Therefore, the scaling properties of the wind–terrain–vegetation interactions are a critical focus for understanding patterns of snow accumulation over different scales.

Efforts to model the snow depth or snow water equivalent (SWE) distributions resulting from this complex system, while undoubtedly successful, show that there is still progress to be made. Snow distributions have been modeled extensively using terrain components as predictors in diverse model types including simple linear regression, binary regression trees (Elder et al. 1995), and geostatistical methods (e.g., Balk and Elder 2000), predicated on the wind–terrain interactions that govern snow redistribution and the relative ease with which terrain can be modeled. Various methods of parameterizing terrain have been explored, including using “curvature” (the first derivative of slope or aspect) (e.g., Blöschl et al. 1991; Liston and Sturm 1998), shelter from versus exposure to wind (e.g., Purves et al. 1998), and characterization of upwind and downwind terrain features (e.g., Winstral and Marks 2002; Winstral et al. 2002). The work of Winstral et al. (2002) confirms the notion that a terrain parameter most suited for addressing snow accumulation addresses the upwind terrain features that directly modify the wind field on the 10–100-m scale. They demonstrated that a “wind exposure index” can identify areas of convergent and divergent wind flow, and consequently predict areas of snow scour and deposition.

Winstral et al. (2002) also state that the relatively low portion of the snow depth variance explained by the terrain parameters could be related to differences in process and model scales. This issue, treated in detail by Blöschl (1999), is a critical problem for the development and testing of spatial snow accumulation and melt models. The true spatial variability is never fully knowable, and our assessment and characterization of it is limited by the scales at which we observe snow properties, and the scale at which we model them.

Measurement of spatial phenomena is in essence a filtering of the actual variability on the ground, as the resolution at which measurement occurs is usually different from the actual process scale (Blöschl 1999). Depending on the true nature of the spatial phenomenon, sampled variability may be either smoother or rougher than what exists on the ground. To assess this scaling effect induced by the data collection, we must therefore know the process scale. The most common method used to estimate the process scale is calculation of the correlation length from a variogram. The correlation length is the lag distance (length) beyond which there is no correlation between adjacent points. In practice, the correlation length is defined as the scale where the semivariance reaches 95% of the spatial (sill) variance and only exists for stationary (constant mean) or locally stationary processes where the variogram has a well-defined sill (Webster and Oliver 2001).

The correlation length, while useful for discriminating autocorrelated from uncorrelated scale regions, cannot describe any scaling relationships at longer or shorter scales. The slope of a log–log variogram, and therefore its fractal dimension, may provide an important alternative to the correlation length (Blöschl 1999). Fractal analysis performed on a variogram, in addition to identifying important breaks in scale analogous to the correlation length, can also characterize the scaling properties of the spatial distribution in question. Further, if the fractal dimension changes with scale (i.e., the phenomena is multifractal), the length scale of the changes can point to important changes in the process relationships that create the spatial pattern in question.

Other studies have found self-similar, or fractal, behavior in snow properties. Tabler (1980) demonstrated consistent wind drift geometry over several orders of magnitude, from centimeter-scale models to drifts behind 4-m snow fences. This scale invariance suggests that wind-induced patterns of snow accumulation may exhibit scaling properties over a significant scale range, with important implications for analyzing the spatial distribution of snow in mountain environments. Shook et al. (1993) found that patchy snow covers show fractal scaling in the area–frequency relationships of snow and snow-free patches. They attributed this effect to a fractal distribution of SWE. Shook and Gray (1996) analyzed snow depth transects in southern Saskatchewan, Canada, and found that snow depth shows a fractal distribution over lengths from 1 to around 30 m, with spatially random behavior at longer distances. Subsequently, Shook and Gray (1997) used fractal distributions to model snow covers in an effort to supplement available field data for snowmelt model development. They also describe the fractal nature of patchy snow covers and note that the fractal patches arise from the fractal spatial distribution of SWE. Blöschl (1999) examined snow-covered area data from diverse sources across a wide range of scales and showed that power-law fractal patterns could exist over several scale regions covering multiple orders of magnitude. Kuchment and Gelfan (1997, 2001) demonstrated fractal behavior in both linear snow course and spatially aggregated snow course data and used the fractal scaling relationships to estimate spatial snow depth variability at different scales. Granger et al. (2002) also examined patchy snow covers and found the patch length–area relationship to be fractal. In contrast to previous work, the rich spatial coverage of the lidar dataset used in this study allows a more robust spatial analysis than transect data. Furthermore, the areal coverage of the data allows an analysis of directionality in the spatial patterns of snow depth.

### b. Fractals and fractal measurement

The spatial distribution of snow depth is the result of the complex interactions of many physical processes. Some of the interactions are highly nonlinear or chaotic, producing a pattern with such complexity that it appears random. Such stochastic variation is often treated with geostatistical methods, employing the variogram and correlation length, as described above. Geostatistical techniques characterize a spatially varying phenomenon as a random field, that is, the sum of a mean value and a random component (Webster and Oliver 2001) where the spatial pattern is composed of random deviations from the mean value. The random characterization neglects any spatial structure that may exist in the random component.

Treating snow depth in the above manner implicitly limits the explainable variability by defining the scale at which the mean value is computed and does not allow for multiscale effects arising from process interactions over multiple length and time scales. A fractal analysis, conversely, allows for stochastic structure within the random component of the distribution and explicitly allows multiple scales of variability. Scale invariance, or fractal geometry, has been observed in snow distribution patterns (Shook and Gray 1996; Kuchment and Gelfan 1997; Granger et al. 2002) over specific scale ranges.

A fractal is a structure that has infinite detail—that is, more detail is resolved as the scale of observation is reduced (resolution increased). The structure is therefore *self-similar*—small parts of the structure have the same shape and/or statistical properties as the whole. Mathematical fractals display perfect self-similarity. In practice, no perfect natural fractals have been found, though numerous natural systems show fractal properties over distinct scale windows. Fractal concepts have seen numerous applications in the earth sciences, including coastline geometry (Mandelbrot 1983; Phillips 1986), earthquake size distribution (Turcotte and Huang 1995), landslide magnitude (Malamud and Turcotte 1999), and precipitation events (Peters and Christensen 2002).

A fractal is quantified by its Hausdorff–Besicovich, or fractal, dimension (*D*). Fractal dimensions are consistent with Euclidean (topological) dimensions, but are noninteger. The fractal dimension of a linear feature has a value between 1 and 2, with numbers near 1 indicating a smooth line and values near 2 describing a nearly plane-filling object. Similarly, surface features have *D* values between 2 and 3, with low values being close to planar, and high values being nearly volume filling. Indeed, a fractal is defined as a feature whose Hausdorff–Besicovich dimension exceeds its topological dimension (Mandelbrot 1983; Lam and DeCola 1993); *D* therefore is an index of the complexity or roughness of the feature or distribution in question (Turcotte and Huang 1995; Gao and Xia 1996).

Natural features most often exhibit fractality only over a specific range of scales (Mark and Aronson 1984; Gao and Xia 1996). Often features will have one or more scale windows over which they exhibit fractal self-similarity, with a distinct fractal dimension characterizing each window. This multifractility is not unexpected, as different driving processes are likely to become dominant at different scales in any given system (Blöschl 1999; Phillips 1999). These “scale breaks” can be identified by a change in the fractal dimension value and can indicate important changes in process or process interactions.

Complex, nonlinear systems have been shown to produce fractal distributions, indicated by power-law behavior in size–frequency relationships (Malamud and Turcotte 1999). Turbulent wind fields behave chaotically, with vortices resulting from the nonlinearities of viscous fluid dynamics in the near-surface atmosphere (Poveda-Jaramillo and Puente 1993). Therefore, in environments where wind plays a significant role in the spatial distribution of snow, snow depth is expected to show a fractal distribution over scales from 1 cm to 1000 m, where wind redistribution dominates.

If the spatial distribution of snow depth has fractal characteristics, the spatial distribution has a structure that produces similar statistical snow depth distributions at multiple scales, allowing spatial pattern characteristics to be transferred from one scale to another using a scaling parameter (Kuchment and Gelfan 1997). Fractal analysis can also provide information about the scale, scope, and resolution of modeling and field sampling efforts, as well as reveal any particular scale ranges over which snow depth distributions can be up- or downscaled. Knowledge of the spatial scales over which a phenomenon shows a fractal distribution also exposes scales at which important changes in governing processes occur (Burrough 1981).

### c. Objectives

Previous analyses applying fractal techniques have important ramifications for up- and downscaling of snow depth distributions, sampling design, and the appropriate resolution of snow accumulation models. This study builds on these previous efforts, using high-resolution, spatially extensive data that allow a detailed analysis over a large area and in multiple directions. We examine the fractal characteristics of snow depth distributions and the associated terrain and vegetation patterns. Lidar-derived snow depth, topography, and vegetation topography datasets from three study sites are subjected to variogram analysis to look for fractal patterns and scale limits to fractal behavior. The variogram analysis is next restricted to 16 compass directions to examine anisotropy in the datasets. The fractal characteristics of snow depth, topography, and vegetation topography are then compared to look for any scaling characteristics consistent among the three datasets.

## 3. Study areas

This study uses data from the 2003 National Aeronautics and Space Administration (NASA) Cold Land Processes Experiment (CLPX) in Colorado, at the Buffalo Pass, Walton Creek, and Alpine Intensive Study Areas (ISAs; Cline et al. 2001; Figs. 1 and 2). The Buffalo Pass site is characterized by dense coniferous forest interspersed with open meadows, low rolling topography, and deep snowpacks. The Walton Creek site provides a broad meadow environment, interspersed with small, dense stands of coniferous forest, low rolling topography, and deep snowpacks. The Alpine study site contains alpine tundra, with some subalpine coniferous forest, and is generally north-facing with moderate relief (Cline et al. 2003). Buffalo Pass and Walton Creek are in a similar synoptic-scale terrain position, receiving high annual snowfall, while the Alpine site receives lower annual precipitation totals and has greater wind exposure above tree line.

The three study areas provide a variety 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 all sites are dominated by wind redistribution and wind interaction with terrain and vegetation patterns. Wind direction frequency distributions, calculated from meteorologic data collected at each of the study sites for a period of 1 October 2002 through the lidar flight date of 9 April 2003 (Elder and Goodbody 2004), indicate that winds of speeds greater than 5 m s^{−1} (at 10-m height) are confined to narrow direction bands (Fig. 3).

## 4. Methods

### 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 data used here have 1.5-m nominal horizontal spacing between data points and a 0.15-m vertical measurement tolerance. Positioning data were collected using a combination of global positioning system and inertial navigation system technologies. The raw data were normalized using ground control points and then postprocessed to remove redundant data points and noise (Miller 2003).

Two lidar datasets are used for this study—one acquired on 9 April 2003 and the other on 19 September 2003. The last-return signal from the September mission provides ground surface (“bare earth”) elevations, the September first-return data measure the terrain-plus-vegetation (vegetation topography) elevations, and the April 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 a dataset of snow depth point estimates. The point datasets of snow depth, bare earth elevations, and vegetation topography were used for the variogram fractal analysis.

### b. Variogram fractal analysis

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 *k*th bin, and *z _{i}* and

*z*are the snow depth values at two different point locations

_{j}*i*and

*j*. 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. The variograms were restricted to a maximum distance of 1100 m, the diameter of the largest circle that could 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 (*R*^{2} 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 more conservative estimate of the scale break location.

## 5. Results

### a. Snow depth

All three datasets show strong power-law distributions over two segments, separated by a relatively distinct scale break (Fig. 4). The short-range window is characterized by fractal dimensions of approximately 2.5, while at longer distances the *D* values approach 3, indicating nearly spatially random distributions (Table 1).

The directional variograms (Figs. 5a and 5b) show short-range *D* values that are smaller parallel to the dominant wind direction, showing smoother (i.e., autocorrelated) patterns with the wind and rougher patterns normal to the wind (see Fig. 3). The long-range, directional *D* values (Fig. 5b) show the opposite pattern of the short-range values, namely, that the smallest fractal dimensions are found perpendicular to the dominant wind direction.

### b. Bare earth topography

The topography datasets show power-law behavior over several orders of magnitude, from around 10 m to over 1000 m (Fig. 6; Table 2). The power-law fits are strong, with *R*^{2} values greater than 0.90. The Buffalo Pass and Alpine variogram plots are nearly coincident, with very similar *D* values. The Walton Creek terrain data show a much higher fractal dimension than the other two sites, indicating a much higher relative roughness. The small-scale terrain variations are large relative to the overall relief of the area. In contrast, the Buffalo Pass and Alpine sites have a large-scale relief that is large relative to the small-scale terrain roughness.

The directional terrain variograms show strong anisotropy (Fig. 7). The smallest *D* values, and thus the smoothest relative roughness, are found parallel to the direction of maximum relief.

### c. Vegetation topography

The Buffalo Pass and Walton Creek vegetation topography data show similar scale breaks, at 50 and 60 m, respectively, while Alpine data have a scale break at 30 m (Fig. 8; Table 3). At scales below the break, the fractal dimensions are consistent between sites, varying from 2.81 to 2.86. Beyond the scale break, the distributions are virtually identical to the pattern in the bare earth topography data.

The directional variograms show that the short-range *D* values have little variation with direction, with the exception of a decrease in *D* values in the north–south direction at Buffalo Pass, likely the effect of the ribbon forest orientation at that site (Fig. 9a; Fig. 2a). The long-range fractal dimensions show strong anisotropy, with directional patterns very similar to the bare earth terrain data (Fig. 9b).

## 6. Discussion

The results of the fractal analysis show that the snow depth and terrain datasets are well described by the fractal model over specific scale intervals. The fractal dimensions obtained are consistent with previous studies, as are the scale windows over which snow depth exhibits fractal behavior. There is distinct anisotropy in the fractal dimensions. Topography and vegetation topography also fit the fractal model well; vegetation topography shows a scale break at a scale magnitude compatible with that observed in the snow depth data, while bare earth topography shows no scale break within the range of scales available in the data.

### a. Omnidirectional distributions

The scale break observed in the snow depth distributions indicates a change in the process dynamics controlling snow accumulation. The spatial and/or functional relationships between wind, terrain, and vegetation, which combine to produce the snow distribution patterns observed at a particular scale, are substantially different when observed at resolutions above and below the scale break. The values of the snow depth fractal dimensions on either side of the scale break lend insight into how the resulting patterns are different and may guide investigation into the actual process dynamics.

The short-range snow depth *D* values are nearly identical at all three sites, indicating that the process interactions producing the snow depth patterns are essentially equivalent in all three areas. This similarity in *D* values also implies that the physical characteristics that differ between sites are insignificant in determining the snow depth patterns over short distances. Fractal dimensions have been used to differentiate features that are formed by different process interactions (Burrough 1993); the spatial organization indicated by the fractal distribution may be characteristic of snow accumulation patterns in wind-dominated areas when observed over a square-kilometer extent. The snow depth *D* values of 2.5 indicate that at length scales below the scale break, variation due to short-range processes, such as snow drifting around individual trees, is balanced by variation caused by longer-range process influences such as the direction of wind flow above the canopy (Burrough 1993; Kuchment and Gelfan 1997; Blöschl 1999). At distances longer than the scale break, the snow depth fractal dimensions vary from 2.91 (Alpine) to 2.97 (Buffalo Pass). The high *D* values indicate that the distribution in this upper fractal window is dominated by high-frequency variations (relative to the power-law region). Though the patterns approach spatially random distributions, the spatial structures indicate fractal distributions. The high values suggest, however, that spatial interpolations from sparsely distributed data points may not be appropriate.

The fractal dimensions of the different scale windows have important ramifications for spatial modeling and mapping. To interpolate a surface from point observations, the parameter variability as a function of the distance from a data point must be known or estimated. If the spatial distribution is dominated by short-range variation (high *D* value), then the point values are essentially independent, and the variable is probably better characterized as a random distribution. The fractal dimensions calculated for snow depth suggest that this is the case at length scales longer than the scale break. However, at length scales shorter than the scale break, the *D* values indicate a near-equal balance of long- and short-range variation—more spatial structure (Burrough 1993). Spatial interpolation may therefore be reasonable over the scale range below the break point.

The vegetation topography dataset represents the surface roughness elements that wind and blowing snow interact with directly. The bare earth topography certainly influences snow accumulation, predominately through flow concentration and separation, and orographic precipitation mechanisms. However, snow depth variability is dominated by blowing snow–vegetation interactions across a wide range of scales dependent on vegetation characteristics interacting with local topography. Therefore, the scaling properties of the vegetation topography dataset are critical to investigating scale characteristics of snow depth.

The vegetation topography data show a scale break at the same order of magnitude as snow depth, though at a slightly longer absolute range. This scale break indicates a switch from a pattern dominated by short-range vegetation variability to one predominately influenced by longer-range variations in terrain. A pattern change in the interaction of terrain and vegetation processes would no doubt influence the patterns of snow deposition and redistribution observed at scales straddling the scale break. All three sites show a scale break separating a scale region of high-frequency variability from a region of terrain-influenced low-frequency variability.

The bare earth terrain shows no scale break over the scale range available in the data. Other studies have shown that terrain patterns in a wide variety of physiographic environments typically exhibit a scale break at lengths around 1.5 to 2 km (e.g., Klinkenberg and Goodchild 1992; Mark and Aronson 1984), which is larger than the size of the study area used here. An analysis of vegetation height data (vegetation topography with bare earth terrain subtracted), also does not display a scale break. That the vegetation topography data do show a sharp break in slope, while the individual terrain and vegetation height datasets do not, indicates that a process interaction changes at that particular length scale.

Kuchment and Gelfan (2001) analyzed spatial patterns of snow depth from transect data at slightly longer scales than presented here, and so were unable to resolve any scale breaks that may be present in the 10–50-m range. However, for longer length scales they produced fractal dimensions of 1.82–1.93, similar to those found in this study (linear transects will have a fractal dimension between 1 and 2, and adding 1 to the *D* value will give the equivalent fractal dimension for a surface). Their study sites were diverse, varying from low to high relief, and steppe to forest to tundra environments. Shook and Gray (1996) also used transect data to compute fractal dimensions and scale breaks (cutoff lengths, in their study), with results (*D* = 1.53, scale break at 30–500 m) similar to those obtained here. The fractal slope at distances longer than the scale break was also nearly random (*D* = 1.94). In both of the above studies, the terrain and vegetation in the study areas were significantly different from the sites used here. The consistency of results between all three studies suggests that the spatial properties of snow depth can be consistent across different environments, given a strong dependence on wind redistribution.

Shook and Gray (1996) demonstrated differences in the scale break length between sites and showed a dependence of scale break distance on topographic relief. Our results follow a similar scale break pattern, with the scale break distances increasing with the overall study area relief (Table 1). However, there are several other differences between the three study areas that could impact the scale break location. For instance, the forest structure is dramatically different at each site, with dense stands and open meadows at Walton Creek, sparse stands and ribbon forests at Buffalo Pass, and dense subalpine forest adjacent to treeless tundra at the Alpine site. Notably, the Buffalo Pass and Walton Creek snow depth scale breaks are essentially equivalent, while the scale break at the Alpine site is at a much longer distance. This implies a regional difference, perhaps related to the vegetation–terrain interaction. The vegetation topography scale breaks are also nearly identical at Buffalo Pass and Walton Creek, and are at the same order of magnitude as the snow depth scale breaks. These corresponding scale features indicate that the vegetation–terrain interaction at Buffalo Pass and Walton Creek exhibit similar scaling properties, and that the scale patterns of vegetation topography are qualitatively related to the regional differences in scale break length. Therefore the scale break differences cannot be solely related to relief effects, though the influence of broad-scale terrain trends is undoubtedly important.

The fractal spatial distribution of snow depth at short scales suggests that rescaling of observed distribution patterns would be possible as long as the resolutions or grid sizes of the original and rescaled patterns both fall within the same fractal window. For example, a spatial distribution with a 30-m resolution could be rescaled to a 10-m resolution at the Alpine site, but not at either the Buffalo Pass or Walton Creek sites. The scale window extent also has implications for sample design. It is likely that as sample spacing is increased toward the upper limit of the scale window, the fractal slope and the location of the scale break will become progressively less resolvable, and therefore representation of the spatial structure will become increasingly poor. On the other hand, if resolutions above the scale break are desired, the data suggest that upscaling from 1000-m resolution is possible, allowing for less expensive data collection at coarser resolutions. However, since the *D* values in the upper scale range are very high, the spatial distribution of snow depth is nearly random. Investigations conducted at resolutions in this range would provide information on the magnitude of snow depths, but will show very little spatial structure.

### b. Directional distributions

The large volume of data present in the lidar datasets enables an investigation of directional anisotropy in the magnitude of the fractal dimensions. The snow depth data show opposite patterns between the short-range *D* values and the long-range *D* values; the process change that occurs at the scale break appears to alter the relationships producing the directional fractal dimension patterns.

For distances shorter than the scale break, snow depth fractal dimensions are larger in directions normal to the prevailing wind. These length scales are of a magnitude consistent with snowdrifts behind trees and small terrain features. A depth transect at the 1–15-m scale measured across a series of drifts would show a higher roughness (negative correlation) than a transect lengthwise along an individual drift, and thus a higher fractal dimension (e.g., see the drift patterns in Fig. 2a). At longer length scales, the largest snow depth *D* values occur parallel to the dominant wind direction. The large-scale trends are influenced by large terrain trends, such as those produced by orographic precipitation or lee slope drifting, and thus would exhibit maximum roughness in the direction consistent with the dominant wind-steering terrain feature.

The terrain data show anisotropy that is also similar to, and likely induces, the dominant wind direction at each site. The lowest fractal dimensions are found parallel to the direction of maximum relief at each site, implying a lower roughness in this direction. Despite this directional variation, the *D* values are relatively small in all directions, indicating that the high persistence and the dominance of low-frequency variability in the terrain data are relatively consistent in all directions.

The vegetation topography data are essentially isotropic at scale ranges shorter than the scale break. The exception is Buffalo Pass, which shows a distinct decrease in fractal dimension in the north–south direction. This anisotropy is likely due to the presence of north–south-oriented ribbon forests in the northeast quadrant of the study site (Fig. 2a). For scale lengths longer than the scale break, vegetation topography shows anisotropy very similar to the terrain data. This similarity is not surprising, as the terrain variability is dominating the vegetation topography fractal signal in this scale window.

As with the omnidirectional results, the anisotropic behavior of the snow depth spatial distributions could be important for measuring and modeling snowpacks. Sampling schemes could have varying spatial resolutions dependent on orientation to dominant wind patterns. Models representing accumulation processes, especially wind, may also require parameters that change with direction. A better understanding of the actual process dynamics that cause the scaling features and directional anisotropy is required if these concepts are to aid practical application of measurement and modeling of snowpack processes.

## 7. Conclusions

The complex process interactions that produce the spatial distribution of snow in mountain environments change with scale. Sampling, modeling, and prediction efforts therefore require an assessment of the scaling properties of snow distributions in order to capture and represent natural variability in the snow cover. Fractal analysis techniques applied at the local scale (1–1000 m) show promise for identifying scales at which snow-related processes or process interactions change, and for finding scale regions over which rescaling of measured or modeled distributions may be reasonable.

Snow depths, as derived from lidar altimetry surveys, show two distinct scale regions with fractal distributions, separated by a scale break. At small distances, the fractal dimensions are close to 2.5 for all study sites, indicating rich spatial complexity and a balance between high- and low-frequency variations. At larger distances, snow depth shows a very high fractal dimension (>2.9), approaching a spatially random distribution. The length of the scale break between fractal windows appears to vary with overall terrain relief; longer break distances are found in higher-relief terrain. These results are consistent with previous studies in diverse environments; however, it is important that relief is not the only characteristic that varies between sites, and other factors, such as local snowfall characteristics, could also be influencing the scale break length.

Terrain and vegetation interactions also show fractal distributions over distinct scale regions, with vegetation effects dominating at small distances and terrain effects dominating at large distances. The length of the scale break separating the two fractal regions in the terrain–vegetation distributions is of the same order of magnitude as the scale break observed in the snow depth data, which indicates that the process change revealed in the vegetation–terrain data potentially influences the scaling behavior of snow depth patterns.

The snow depth fractal dimension varies with direction and shows a strong qualitative relationship to prevailing winds and large-scale topographic orientation. The anisotropy indicates that directionality of controlling processes must be considered in pattern-analysis application or in sampling design. Terrain and vegetation topography fractal dimensions are also anisotropic; the directional variation appears to be related mainly to the terrain orientation, though the Buffalo Pass site shows some variation in the vegetation distribution associated with the presence of ribbon forests.

Efforts to assess, model, and predict spatial snow distributions are increasingly needed for water supply estimation, river forecasting, and avalanche danger evaluation. To improve these efforts, the behavior of controlling process interactions at different scales must be addressed. The fractal techniques presented here provide important measures of process change scales and scaling behavior in snow distributions. Investigations of multiscale spatial patterns as applied in this study offer a necessary perspective on snow system behavior across multiple process regimes and can distinguish scales where an ordered pattern exists from scales where disorder dominates. Quantifying system interactions at multiple scales is critical for understanding the complex relationships between wind, terrain, vegetation, and snow distributions.

## Acknowledgments

This project was supported by NASA Headquarters under Earth System Science Fellowship Grant NNG04GQ40H. The lidar and meteorologic data were collected as a part of the NASA Cold Land Processes Experiment (CLPX), coordinated by Don Cline (NOAA/NOHRSC) with additional support from Robert Davis (U.S. Army CRREL).

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Jeffrey S. Deems, Department of Geosciences, Colorado State University, Fort Collins, CO 80523-1482. Email: deems@cnr.colostate.edu