Mark W. Williams¹, Richard Sommerfeld², Sam Massman¹, and Mark Rikkers¹
¹Institute of Arctic and Alpine Research and
Department of Geography, CB 450,
University of Colorado, Boulder, CO
²Rocky Mountain Forest and Range Experimental Station
US Forest Service, Department of Agriculture
Fort Collins, CO
The melting of snow is known to be an inhomogeneous process. The spatial distribution of meltwater flowing from the bottom of melting snowpacks is the result of horizontal and vertical flow paths within the snowpack. The ability to characterize the spatial distribution of these meltwater flowpaths would be useful in developing snowmelt runoff models which could better characterize snowmelt hydrographs. Near infrared aerial photos of melting snow have been analyzed using a moving window analysis which can characterize correlation lengths in the reflectance of the snow surface. Near infrared is sensitive to snow grain size which indicates the concentration of melt water; the grains grow faster if the liquid water content is higher. The probability of finding such correlation lengths was about 0.22 in May, 1997 when the melt had just started and rose to 0.68 by June when the melt was well established. Correlation lengths for all sampling dates ranged from 5 to 7 m. Liquid water content at the snow surface was sampled with a dielectric sensor at 0.5 m intervals on two 100-m² grids. Semi-variograms showed a sill at 5 to 6 meters. The liquid water measurements at the snow surface suggest that the correlation lengths derived from the infrared aerial photos represent surface expressions of vertical flowpaths through the melting snowpack. A circular array of 16 snowmelt lysimeters each with areas of 0.2 m² was operated for two years at Niwot Ridge in the Colorado Front Range. Variograms indicated that flows were correlated over a distance of 5 to 7 m. These three independent methods all suggest a correlation length of 5-7 m for vertical flowpaths draining ripe snowpacks in the Rocky Mountains.
The melting of snow is known to be an inhomogeneous process. The spatial distribution of meltwater flowing from the bottom of melting snowpacks is the result of horizontal and vertical flow paths within the snowpack [Marsh and Woo, 1985; Kattelmann, 1985; Marsh, 1988]. Simple tracer experiments using dye applied to the snow surface have shown that liquid water often flows through a snowpack in preferential pathways, becoming concentrated at the interfaces of horizontal layers and then breaking through in preferential channels [Gerdel, 1949]. However, attempts to characterize the spatial distribution of preferential flowpaths have had only limited success [e.g. Marsh and Woo, 1985, Kattelmann 1989]. Attempts to understand meltwater flow through snow from first principles have also had only limited success [e.g. Colbeck, 1979, 1991]. The ability to characterize the spatial distribution of these meltwater flowpaths would be useful in developing snowmelt runoff models which could better characterize snowmelt hydrographs.
The grain size of wet snow increases rapidly compared to drier snow [Wakahama, 1968, Colbeck 1979]. The reflectivity of snow decreases with increasing grain size in the near infrared because ice is moderately absorptive in those wavelengths [Nolan and Dozier, 1993]. Sommerfeld et al.  hypothesized that vertical flowpaths would cause larger snow grains on the snow surface compared to surrounding areas of the snow surface that did not have vertical flowpaths. In turn, the larger grain sizes caused by the concentration of liquid water at the top of these vertical flowpaths should cause a bimodal brightness distribution in infrared aerial photos of melting snow. They further hypothesized that the separation of the modes would be maximized in images spatially filtered at a spacing characteristic of the spacing of the melt channels. They supported this hypothesis with some data in both visible and infrared aerial photos. The resulting characteristic distance was about 6 m, consistent with results from a lysimeter array [Sommerfeld et al., 1994]. However, the quantity of their data was small and the question of the consistency of the characteristic distance over wider areas remained open. Whether the observations reported by Sommerfeld et al.  actually reflected differences in surface wetness was not answered.
Our objective here is a more robust test of the hypothesis first proposed by Sommerfeld et al.  that vertical flowpaths in melting snowpacks have a characteristic length of about 6 m. We test these hypotheses:  The characteristic distance separating bimodal reflectance peaks from aerial photographs of a draining snowpack in Colorado would be the same as that reported by Sommerfeld et al.  for a Wyoming snopack.  The frequency of bimodal distribution of brightness values from aerial photographs of snow covered area would increase from early in the melt season to later in the melt season.  Direct measurements of the liquid water content of the snow surface will have the same characteristic distance as bimodal distributions of brightness values.  The characteristic length of meltwater discharge measured in snow lysimeters at the base of the snowpack will have the same characteristic length as the bimodal brightness distribution of aerial photographs.
All experiments were conducted on Niwot Ridge, located in the Colorado Front Range of the Rocky Mountains about 5 km east of the Continental Divide (40䓃'N, 105䓣'W). This site is an UNESCO Biosphere Reserve and a Long-Term Ecological Research (LTER) network site. Niwot Ridge is an alpine interflueve and was not glaciated during the Pleistocene; and extends from the continental divide at an elevation of 3,700 m east about 6 km to an elevation of 3,400 m. Climate is characterized by long, cool winters and a short growing season (1-3 months). Since 1951, mean annual temperature is -3.8°C and annual precipitation is 1,000 mm [Williams et al., 1996a].
Figure 1. Oblique aerial photograph taken in April near maximum snow accumulation. Niwot Ridge is the broadening ridge in the center of the photo and running left to right. There are meteorological stations at D1 and the Saddle site.
Snowpack conditions are monitored about weekly at the Saddle site on Niwot Ridge (Figure 1) following the protocol of Williams et al. [1996b]. Measurements include snow depth, snow density at ten-cm increments, snow temperature every ten cm, snowpack stratigraphy, and grain size. Depth-weighted values are then calculated for snowpack density and temperature.
Aerial photographs were acquired using a 70 mm Kodak Infrared Aerographic type 2424. The camera was a motor drive Hasselblad, with a 250 mm lense. Visible light sensitivity was attenuated with a yellow filter (Wratten No. 12) resulting in a spectral response of approximately 900 to 1200 nm. Aerial photographs were taken at an elevation of 5,000 m on three dates: 14 June 1995, 13 May 1997, and 24 June 1997.
Black and white negatives from the three flights were digitized using a Sony 3 CCD color video camera with a resolution of 640 X 780 pixels per frame. The resulting scale factor, 570 mm pixel, was determined by ground and map measurements of identifiable features. Camera images were acquired in ImagePro Plus (Media Cybernetics, 1995) as grey scale images. The software includes a fast Fourier transform capability. The resulting two dimensional spectrum was edited in the low pass mode to selectively remove higher spatial frequencies. Each frame included nine moving window analyses at three different pixel resolutions: 12.2, 25.6, and 51.2.
No formal method was used to randomize this survey. However, the aircraft pilot made no attempt to acquire images of specific areas and the only criteria used to select frames for analysis was (1) the frame was largely snowcovered and (2) it did not significantly overlap other analyzed frames. It is thus unlikely that there was any bias in the image selection. In all, about 50 frames were examined.
The liquid water content of the snow surface was measured with a portable dielectric sensor, a Finnish Snow Fork (FSF) [Shivola and Tiurri, 1986]. The FSF measures the liquid water content at the snow surface to a depth of 6 cm. The FSF was used to measure the liquid water content of near surface snow because it is fully automated with a fast response time. Comparisons of liquid water measurements in snow between the FSF and the more standard Denoth Wetness Sensor [Denoth et al., 1984] showed no significant difference in either mean values or range [Rikkers, 1997]. Measurements of liquid water content using the FSF were made every 0.5 m on a 10 m x 10 m grid and a 10 m x 7.5 m grid in July of 1996.
The spatial autocovariance of the liquid water content was determined using a semivariogram as suggested by Carr :
During the summer of 1994, a circular array of sixteen snow lysimeters, each 0.2 m2 in area with 25-cm sidewalls, was installed with a 5 m radius and approximately 2 m spacing along the arc of the circle, as reported in Rikkers et al. . The lysimeters were individually plumbed and output routed by gravity into dedicated tipping buckets (Davis Rain Collector II) in a subnivean laboratory. All buckets were mounted and leveled inside the subnivean laboratory and drain funnels directed all meltwater through the floor to a gravel bed below.
Individual tipping buckets were wired to a Campbell Scientific CR10 data logger and the data storage palm-card was replaced and downloaded weekly. Meltwater flux was recorded as the total number of tips at ten-minute intervals. Tipping buckets were calibrated to volume by slowly pouring a known amount of water into the buckets and recording the number of subsequent tips; one tip was 4.6 ml (precision of about 10%).
The spatial statistics of meltwater discharge were calculated using a spatial variogram which compares the population variance with the covariance of fluxes at different linear lag distances, following the protocol of Sommerfeld et al. (1994):
where h is the lag distance between pairs of lysimeters, n is the total number of points, n h is the number of pairs, Z(xi) are the measured flows and Z--- is the mean of the measurements. If the covariance (2nd term) is smaller than the population variance (1st term), gamma(h) will be positive. We take the lag where gamma(h) crosses 0 to be the correlation length. The spatial variability of meltwater flux was calculated for linear lag distances of 0, 1.9, 3.8, 5.5, 7.1, 8.3, 9.2, 9.8, and 10.0 m. Computations were made using a computer script (awk) which ran the spatial variogram on each possible combination of pairs at each lag distance at ten minute intervals.
The 1997 snow year was about average, with a maximum snow depth of 2.57 m occurring on 7 May at the index site on the Niwot Ridge saddle. The snowpack became isothermal at 0°C on 10 May, with complete meltout in the second week of July. Maximum snow water equivalence (SWE) was 1,016 mm. The 1995 snow year was characterized with spring snow and a later melt date, with snowmelt starting about 1 June [Cline, 1997].
Figure 2 illustrates the process of finding the spatial distribution of brightness values in the snow covered images from 1995. North is to the top of the scene and the photos were taken about 11 AM. The smoothed area was spatially filtered using a Fourier transform of 25.6 X 25.6 pixels. The Fourier transform window contains both the light and the dark regions so that the correlation length is one-half the length of the Fourier window, or 7.5 m. There is a tri-modal distribution of brightness values, shown in the inset histogram. Three different classes of snow crystal sizes are probably represented in the scene: (1) large crystals in the wetter areas; (2) intermediate size crystals in the dryer areas and (3) fine grain crystals from a light dusting of new snow.
Figure 2. Aerial photographs of snow covered area on Niwot Ridge after a fast-fourier transform. The frequency distribution of brightness values for the windowed area of the image is shown on the inset histogram.
The aerial photographs on 13 May 1997 were taken slightly after the snowpack had become isothermal. The aerial photographs on 24 June 1997 were taken well into the melt season. The aerial photographs on 14 June 1995 were taken about two weeks into the melt season, after a light snowfall of several cm. Bimodal distributions of brightness values were found from all three flights. The probability of finding multimodal brightness distributions on the 13 May 1997 flight was about 0.22; the lowest of the three flights. For the flight on 24 June 1997, the melt season was the most advanced of the three flights. The probability of finding multimodal distributions for this flight was 0.68. The snow for the 14 June 1995 flight was well into the beginning of the melt season but had a small amount of new snow on the surface from a snowstorm two days before the flight. The probability of finding a bi or multimodal distribution is these images was 0.60, intermediate between the two 1997 flights. The correlation between the stage of snowmelt and the probability of finding a multimodal distribution was good, supporting the hypothesis that the character of the brightness distribution was a function of snow wetness and grain size.
For all three flights the best separation of the modes in the brightness distribution was associated with a lag distance of 5 to 7 meters (Figure 3). However, in most frames there was a mult-modal distribution of brightness values. To illustrate, in both Figure 2 and Figure 3 there are multi-modal peaks in the frequency of brightness values. In most cases the lower-magnitude peaks in brightness values were nested within higher-magnitude brightness values with the characteristic distance of 5-7 m (eg Figures 2 and 3). These results suggest that the characteristic distance of 5-7 meters in the bimodal distribution of brightness values may be composed of higher-frequency patterns in vertical flowpaths.
Figure 3. A histogram of brightness values from aerial photographs of snow-covered area, illustrating the bimodal distribution image brightness values. The distance between peaks in the bimodal distribution was 6 m.
Near-surface measurements of snow wetness covered a 10 x 10 m grid (n = 441) and a 10 x 7.5 m grid (n = 336) (Figure 4). Both data sets were normally distributed at alpha = 0.05, with p = 0.08 for the 100 m° grid and p = 0.10 for the 75 m°grid. The average liquid water content was similar for both experiments, 9.1% for the larger grid and 8.8% for the smaller grid.
Figure 4. A 3-dimensional surface plot of liquid water content of the snow surface from a 100 m° grid with 50-cm sampling resolution. The x and y axes are length along the grid and the z-axis is percent liquid water content of the snow surface by volume, as measured with the Finnish Snow Fork.
Semivariogram analysis was similar for both experiments; we illustrate with results from the larger grid (Figure 5). The sill of both variograms occurred at a lag distance of about 5 m, indicating that the liquid water content of the snow surface was autocorrelated from 0 to 5 m. In both cases, there was a small change in slope at about 2.5 m, suggesting that the sill at 5 m may in fact be composed of correlation lengths with a higher spatial frequency. The geostatistical analysis suggests that the liquid water content of the snow surface at this time had a spatial organization of about 5 m, similar to the correlation lengths derived from aerial photographs of the snow surface during snow melt. The variograms appear to indicate an anticorrelation at about 10 m, though the data are limited at this lag distance.
Figure 5. A semivariogram of lag distance in m versus semivariance (gamma (h)) for measurements of the liquid water of the snow surface from a 10 m X 10 m grid with measurements every 50 cm. The semivariogram from a 10 m X 7.5 m grid was virtually identical.
Meltwater discharge was recorded for all 16 lysimeters in 1995, 1996, and 1997. However, meltwater discharge was inconsistent at times in 1996 because of icing problems and these data are not reported. Illustrating with results from 1997 (some of the 1995 data was reported in Rikkers et al. 1996), meltwater was collected in all lysimeters from 11 May to 5 July. Total measured flow in each of the 16 lysimeters ranged from a low of 420 mm to a high of 16,000 mm (Figure 6). Visual inspection of the cumulative discharge in each lysimeter suggests a repeating pattern about every four lysimeters (Figure 6). There was relatively low flow in lysimeters 1-4, much higher flow in lysimeters 5-8, lower flow in lysimeters 9-12, and higher flows in lysimeters 13-16.
Figure 6. Total measured discharge (mm) in each of the sixteen lysimeters for 1997. Maximum snow water equivalence (SWE) measured at the index snow pit was 1,016 mm. Discharge in lysimeters ranged from 42% of snowpit SWE to 16 times snowpit SWE.
Figure 7 is a plot of gamma*(h) versus lag distance in meters in ten-minute time steps for the entire 1997 and 1995 melt seasons. For 1997, the results indicate a crossover at about 5 m where the gamma*(h) values become negative. The gamma*(h) values then become positive at about 9 m. For the 1995 season, the gamma*(h) values become negative at about 6.5 m. Similar to 1997, the gamma*(h) values in 1995 approach zero at a lag distance of about 9 meters but do not become positive. These lag distances from both 1997 and 1995 suggest a repeating correlation length of 5-7 m, similar to the correlation lengths from aerial photographs and from measurements of the liquid water content of the snow surface.
Figure 7. Variogram for measured discharge for the 1997 and 1995 snowmelt seasons, normalized by dividing by the population variance. Measured lag distances ranged from 0 to 10 m.
The consistency of lag distances of 5-7 m from three independent measurements of vertical flowpaths in a draining snowpack suggest a conceptual model of meltwater flow through snow (Figure 8). Infiltration starts as a wetting front descending more or less uniformly from the melting snow surface layer; surface melt itself may not be uniform because of surface irregularities which influence energy transfer from the atmosphere to the snow surface. Uniformity of vertical flow is interrupted by grain-scale variations in hydraulic conductivity and by hydraulic barriers. At a hydraulic barrier to unsaturated flow, such as a fine-grain to coarse-grain transition or a pre-existing ice layer (layers B, C in Figure 8), the wetting front stops and water content at the barrier increases. Horizontal snow layers can act as laterally connected hydraulic pathways, and allow lateral movement of water by capillary forces and in response to hydraulic gradients. Water infiltrating from above pools at the hydraulic barrier, and can move laterally along the upper surface of the barrier until a saturation level is reached which allows water to break through the barrier to snow beneath. After breakthrough, the affected region of the snowpack may be dominated by vertical water flow, with rapid advection of free water into the pipe, and for initially subfreezing snow, conduction of released heat of fusion into snow surrounding the pipe. These vertical flowpaths have a characteristic distance (wavelength) on the order of meters (Figure 1) and behave analogously to macropores in soils.
We suggest that once these vertical flowpaths become established there is a positive feedback situation such that these vertical flowpaths are reused and become larger over time. The presence of relatively large amounts of liquid water in the vertical flowpaths causes settling at the snow surface and the formation of a local depression. More liquid water becomes entrained in these localized depressions, causing an increase in grain size and a decrease in albedo, further enhancing melt at the snow surface. Within a vertical flowpath, the presence of liquid water causes an increase in grain size and an increase in the size of pore spaces. Consequently, the hydraulic conductivity is greater within the pipe compared to surrounding snow. Meltwater therefore moves vertically through these flow features at a faster rate compared to vertical flow in the surrounding snowpack. Our lysimeter results suggest that once these flowpaths become activated, they are generally re-used through the melt season.
Figure 8. Schematic diagram of water infiltration through snow. Meltwater flux is proportional to the size of arrows and is in the direction of arrows. Lateral flow dominates at hydraulic barriers to flow, such as the transition from fine-grained snow (B) to coarse-grained dry snow (C). Vertical flowpaths occur where the hydraulic barrier is breached and produce draining areas of the snowpack. Wavelength is the characteristic distance between vertical flowpaths.
Three independent techniques were used to test the hypotheses presented in the introduction: (1) low pass spatial filtering to identify multimodal brightness distributions in near infrared images of melting snowpacks; (2) variograms of the spatial distribution of meltwater from the bottom of melting snowpacks; and (3) capacitance measurements of the liquid water content of the upper surface of melting snowpacks. The three techniques support the hypotheses presented in the introduction.
 The characteristic distance separating bimodal reflectance peaks from aerial photographs of a draining snowpack in Colorado would be the same as that reported by Sommerfeld et al.  for a Wyoming snowpack. The characteristic distances from the aerial photographs reported was approximately 6m.
 The frequency of bimodal distribution of brightness values from aerial photographs of snow covered area would increase from early in the melt season to later in the melt season. The probability of finding multimodal brightness distributions increased from 0.22 to 0.60 and 0.68, corresponding to the extent of snowmelt.
 Direct measurements of the liquid water content of the snow surface will have the same characteristic distance as bimodal distributions of brightness values. Although there was a finer structure evident in the liquid water measurements, regions of higher and lower liquid water content segregated into separate areas of about 5 m in diameter. Variograms showed a sill at about 5 m and appeared to indicate an anticorrelation at about 10 m although the data at this distance is limited.
 The characteristic length of meltwater discharge measured in snow lysimeters at the base of the snowpack will have the same characteristic length as the bimodal brightness distribution of aerial photographs. Variograms indicate that correlation length of the meltwater collectors is 5 to 7 m and there is an indication that the melt is anticorrelated at about twice that distance.
The data presented here provide further evidence for the development of a correlation length in melting snow of about 5 to 7 m. These various spatial statistics are similar to those found in a preliminary study in the Snowy Range of Wyoming indicating that they may be of regional extent. An important result of this study is that infrared aerial photos can be used to investigate these spatial statistics over wide regions without labor intensive surface measurements.