Andrew M. Fox

Department of Geography

University of Cambridge,

Cambridge, CB2 3EN,

United Kingdom

Mark W. Willliams

Department of Geography and Institute of Arctic and Alpine Research

University of Colorado, Boulder

Nel Caine

Department of Geography and Institute of Arctic and Alpine Research

University of Colorado, Boulder


Calculated surface melt and lysimeter discharge were used to establish relationships between specific surface meltwater fluxes and meltwater wavespeeds through the snowpack at a continental, alpine site in 1997. These relationships were used to determine a pore-size distribution index (e ) and to calculate equivalent permeability (k e) for the snowpack. Wavespeeds, and equivalent permeabilities derived from them, exhibited a marked decrease over the melt season and were lower then those reported in the literature. Although not a significant trend, e values seemed to increase over the melt season. The relationship between meltwater fluxes and wavespeeds was found to be dependent on the range of meltwater fluxes used. At high values the relationship is more similar to that predicted by a gravity-flow dominant regime, but breaks down for lower meltwater fluxes when capillary-flow may become more important. These results suggest that snow hydrology models need to incorporate capillary-flow processes when melt rates are less than 1 mm hr-1. However, modeling meltwater flow as gravity-flow at higher melt rates appears to be correct. The switch between these two flow regimes will occur diurnally, with capillary-flow during the morning, and gravity-flow dominant later in the afternoon. This would also partially explain daily patterns of chemical concentration in the meltwaters.


Movement of water through snowpacks is one of the least understood aspects of snow hydrology (Richter-Menge and Colbeck, 1991). It has an important influence on the timing and magnitude of snowmelt hydrographs (Caine, 1992) and on biogeochemical and geomorphological processes (Williams and Melack, 1989; Caine, 1995). Much previous research has been undertaken to investigate the physical factors controlling water movement through snowpacks (Colbeck, 1972, 1973; Ambach et al., 1981; Jordan, 1983a; McGurk and Kattelmann, 1988). Results have revealed a complex series of factors involving grain growth, snowpack thinning, water flow through variably saturated porous media, flow fingering and refreezing of meltwaters within the snowpack, which in turn cause changes in its hydraulic and thermal properties. These studies have also shown that natural snowpacks are highly heterogeneous due to microscale variations in density and ice grain structure, combined with macroscopic layering. However, few experimental studies to define factors affecting water flux through layered snowpacks under natural field conditions have been reported since Colbeck (1975) remarked on the need for more such studies in the early 1970s.

Attempts to model meltwater flux through a heterogeneous, layered snowpack include a variety of approaches: (i) modeling flow around individual ice layers (Colbeck, 1973); (ii) simplifying the effect of ice layers and preferential flowpaths by using multiple flowpath models (Marsh and Woo, 1985); and (iii) treating the snowpack as a homogeneous anisotropic media (Colbeck, 1975). Such models provide insight into flow processes within the snowpack but, without knowledge of the extent and thickness of snowpack horizons and layers and the rates at which their properties change over time, become impractical in nearly all circumstances. An outstanding question in snow hydrology modelling is how to determine an appropriate value for the saturated permeability (ks) of layered snowpacks. Point measurements of ks have little meaning if the permeability is used to predict the rate at which meltwater moves through the snowpack, encountering numerous layers with different individual permeabilities. Most studies, whether analytical or numerical, are forced to treat the snowpack as an equivalent, uniform medium in which permeability is considered as a bulk or integrated property of the entire snowpack.

Unlike soils, snow permeability is almost impossible to measure in the laboratory, because of the triple phase nature of snowpacks and changes in density, grain size, type and distribution they experience during melting. As suggested by McGurk and Kattelmann (1986), measuring an equivalent permeability (ke) of an entire snowpack in the field overcomes these problems and allows the detection of bulk changes during the melt season. ke would be a useful parameter for a variety of approaches to modeling meltwater flow through isothermal snowpacks, including gravity-flow models (Colbeck, 1972, 1979; Dunne et al., 1976; Marsh and Woo, 1985) and numerical models (Jordan, 1983b).

The only field method for determining ke uses measurements of surface meltwater flux (U) and the velocity at which this flux moves through the snowpack ([dz/dt]U) (Marsh, 1991). When these data are plotted for periods of declining flux, the intercept of the line can be used to determine a value for ke over the entire snowpack depth, and the slope of the line can be used to determine the pore size distribution index (e) of the snow. e is the exponent in the relationship between relative permeability (kw) and ks and effective saturation (S*) for unsaturated snow:

k w = k s (S*)e Eq. 1

Where kw is relative, or liquid water permeability, ke is snowpack equivalent permeability and S* is effective water saturation of snow, expressed as a fraction of pore volume occupied by mobile liquid.

e has been determined from drainage experiments, with reported values ranging between 1.5 and 4.8 (Marsh, 1991). Both e and ke depend upon physical properties of the snowpack, and may change considerably during the melt season as snowpack metamorphism affects grain size, snow density and ice layer thickness and continuity (Marsh, 1991). However, previous work has not been able to establish any relationships between either e or ke and snowpack properties or amounts of metamorphism (Denoth et al., 1979), although it has been suggested that ke is dependent upon the range of U used in calculations (McGurk and Kattelmann, 1986).

Perhaps most intriguing, McGurk and Kattelmann (1986) report that ke decreased over time in a Sierra Nevada snowpack. The decrease in ke with time and a thinning snowpack is counter-intuitive since one would expect that ke would increase with time as grain size becomes larger and more uniform under a wet metamorphism regime (Colbeck, 1973). It is unclear whether the values for ke and the change in ke over time reported by McGurk and Kattelmann (1986) were site specific or whether they apply to other sites and climatic regimes.

Here we estimate values for ke from a continental, mid-latitude site in the Colorado Rocky Mountains. Measurements of surface meltwater flux, discharge through snow lysimeters, and snowpack properties were conducted in 1997. Specific objectives of our research were to: 1) establish a range of values for e and ke at this site; 2) investigate whether they vary consistently in a way that can be used to parameterize gravity-flow and numerical models during the melt season; 3) investigate whether ke is dependent upon the range of U used in calculations; 4) evaluate whether this new information can provide insights into the physical processes that govern meltwater flow through snow.


The research site is at an elevation of 3517 m on Niwot Ridge on the eastern slope of the Front Range of Colorado (40 03' N, 105 35' W). Niwot Ridge is a 10-km interfluve extending eastward from the Continental Divide, and is characterized by low rounded hills with shallow saddles (Figure 1). Treeline in this area is at approximately 3350 m. The high elevation, exposure, and typically dry atmospheric conditions of Niwot Ridge result in long periods with clear-sky atmospheric transmissivity, high solar insolation, low magnitudes of incident longwave radiation, low air temperatures, and high wind velocities. From 1951 to 1993, mean annual temperature on Niwot Ridge was -3.8o C and annual precipitation was 1,000 mm (Williams et al., 1996a), with about 80% falling as snow (Caine, 1996). Niwot Ridge forms the northern boundary of Green Lakes Valley (Figure 1) and is an UNESCO Biosphere Reserve and a Long-Term Ecological Research (LTER) network site.

The instrument site is located on a flat bench above treeline within a broad saddle of the ridge (Figure 1). A meteorological station, subnivean (below snow) laboratory, snow lysimeter array, and index snowpit site are located on the Saddle area (Figure 1). This research area is flat and continuous with terrain that gradually increases in slope to about 15o at the top of the West Knoll about 400 m to the north west of the Saddle research site. The snow at the research area is not isolated from the surrounding snowpack, following the protocol of Schneebeli (1998) and reported in Williams et al. (1999a).


Lysimeter Discharge

During May and June 1997, meltwater draining from the base of the snowpack was measured using two unenclosed, zero tension snow lysimeters separated by a distance of 1.0 m. Each snow lysimeter was 1 m2 in area with vertical walls 0.20 m in height, following the suggestions of Wankiewicz (1978). The walls on the lysimeter edges prevented loss of water out of the lysimeters as a result of pressure-gradient differences and also prevented the influx of water from saturated zones along the snow-ground interface. The two snow lysimeters were placed out before the first winter snowfall and were located on the ground surface about 4 m from the Saddle instrument tower and subnivean laboratory, a 2.5 m x 2.5 m x 3 m structure buried to a depth of 2 m in the ground. This facility provided a temperature-controlled environment for data collection and equipment access to the seasonal snowpack and the underlying soil. The lysimeters drained by gravity into the subnivean laboratory where discharge from each lysimeter was measured using Campbell Scientific TE525 Tipping Bucket Rain Gages. Ten minute totals from each of the lysimeters were recorded on a Campbell CR21x data logger, and these were aggregated to hourly totals for this analysis. (For a full description of the lysimeter array see Rickers et al. (1996) and Williams et al. (1999a)).

Snowpack properties

Snow depth was monitored continuously beneath the instrument tower using a Campbell Scientific UDG01 Ultrasonic Depth Gauge. This sensor measures the distance from the sensor to the surface; it is currently mounted at 6.0 m above the ground, so during peak snowpack conditions it was 2-3 m above the snow surface. Snowpack conditions are monitored about weekly at the Saddle site on Niwot Ridge (Figure 1) following the protocol of Williams et al. (1996b). Snow water equivalence (SWE) measurements are made using a 1 L stainless steel cutter in vertical increments of 0.1 m. Temperature of the snowpack is measured every 0.1 m with 0.2 m long dial stem thermometers, calibrated to + 0.2oC using a one-point calibration at 0oC. Grain type, size, and snowpack stratigraphy are also recorded. Raw and calculated values, along with additional information on methods, are available at the Niwot Ridge LTER web site (http://culter.Colorado.EDU:1030/Subnivean/) (Williams et al., 1999b).

Energy balance

Surface meltwater flux (U) was calculated from direct measurements of the snow-atmosphere energy exchange. Meteorological measurements are made on a 7 m mast located approximately 15 m from the subnivean laboratory. Upward and downward looking shortwave radiation is measured with a Kipp and Zonen CM14 pyranometer; upward and downward looking longwave radiation is measured with a Kipp and Zonen CG2 pyrgeometer. Turbulent fluxes are computed using the aerodynamic profile method, with one instrument array at a fixed height of 6 m above the ground, and three instrument arrays that are mounted on a mast that is moved by hand to maintain fixed instrument heights above the snow surface of 0.5, 1.0, and 2.0 m. Air temperature and relative humidity is measured at each height with a Vaisala HMP35c probe housed in a Gill 12-plate radiation shield. Wind speed and direction are measured using RM Young 05103 wind monitors which have a range of 0-60 m s-1 and a gust survival velocity of 100 m s-1; calibration is conducted with an RM Young electric motor. Meteorological variables are measured every 10 seconds; means and standard deviations of the data are computed and stored at 10-minute intervals. Data are recorded on Campbell Scientific CR21x and CR10 Dataloggers and stored on a Campbell CSM1 storage module which is changed once a week. Instrument error was presented in Williams et al. (1999b).

Surface meltwater flux (U) was calculated using the latent heat of fusion from the residual terms in the energy balance equation using the protocol in Cline (1997b) and data presented in Williams et al. (1999b):

DQS + DQM = Q* + QH + QE + QG + QR Eq. 2

where DQS is the cold content of the snowpack, DQM is the latent heat flux during melting or freezing, Q* is the net all wave radiation flux, QH is the sensible heat flux, QE is the latent heat flux, QG is the ground heat flux, and QR is the heat advected by precipitation. Cline (1997a) found QG and QR were not significant heat exchanges until the snowpack became very thin, and so are not considered any further here. The sign convention used was that energy transferred from the snow surface to the atmosphere was negative and energy transferred from the atmosphere to the snow surface was positive (Marks et al., 1992). The accuracy of the instruments and their proximity to the lysimeters means that surface meltwater flux could be modeled with good precision over a wide range of fluxes and short time intervals, a considerable improvement over previous studies such as Jordan (1983a) and McGurk and Kattelmann (1986).

Wavespeeds ([dz/dt]U)

For the downward, non-saturated flux of water (ignoring water pressure gradients), Colbeck (1974) found any surface melt flux (U) propagates through the snowpack at speed :

Eq. 3


Where e is the pore size distribution index; ks is the saturated permeability; f* is the effective porosity, a is the flow factor (5.47 x 106 m-1 s-1 at 0oC).

The speed ([dz/dt]U)at which a certain surface meltwater flux (U) moves down through the snowpack was determined in a non-destructive way by timing how long [dt] it took U to travel the entire depth of the snowpack [dz). Surface meltwater flux was calculated for a given time step from the meteorological data, and then lysimeter discharge (Q) was used to detect when this flux reached the base of the snowpack. To identify when a particular U reached the snowpack base, calculated U and Q were converted into percent cumulative curves, to allow comparison, a protocol similar to that used by McGurk and Kattelmann (1986). We then calculated the time lag between cumulative U at a particular time step and the same amount of cumulative Q. Pressure and water flux data presented by Wankiewicz (1979) suggests that there is hysteresis in the relationship between relative permeability and water pressure on the rising and recession limbs of the diurnal hydrograph. To avoid the potential problems caused by hysteresis, [dz/dt]U was calculated on the recession limb of the diurnal lysimeter hydrograph only, when there is a minimum of saturated flow processes, and gravitational effects on flow are most dominant (McGurk and Kattelmann, 1986).

Pore-size distribution index (e) and equivalent permeability (ke)

The pore-size distribution index (e) for the snowpack was calculated for each day on which [dz/dt]U values were calculated. The slope (b) of log-log plots of [dz/dt]U against U gives the flux exponent (e-1)/e in equation (2), such that e = 3 when b is 0.67. Mean U for all values used in calculating the [dz/dt]U - U regression on an individual day was substituted into the calculated regression equation to give [dz/dt]meanU. Values of mean U and [dz/dt]meanU were then substituted into equation (3), and given a value for e, and that a is a constant, rearranging equation (3) gives ks1/ef*-1. This single value characterizes an important property of the snowpack, which Colbeck (1978) suggested should have typical values assigned, and should then be used as a parameter in operational forecasting, changing as the snowpack matures. Having calculated effective porosity (f*) for a given day using density and irreducible water content values following the procedure of McGurk and Kattelmann (1986), the saturated permeability ks for the snowpack can be determined from ks1/ef*-1. As these ks values were determined from wavespeeds calculated for the entire snowpack depth, they represent an 'averaged' value for the snowpack, an 'equivalent permeability' (ke) as defined in Marsh (1991).


The snowpack reached a maximum depth of 2,570 mm on 7 May, 1997 (figure 2) at the index sampling site, with a snow water equivalence (SWE) of 1,070 mm. The snowpack became isothermal shortly afterward on about 13 May. The end of the study period was marked by complete ablation of the on 25 June. Records from a SNOTEL site operated by the US Natural Resources Conservation Service located approximately 5 km southeast of our study site showed that maximum SWE at this location was 140% of the 1980 to 1997 average. In 1997, snow melt rates were relatively low for the first 30 days after peak accumulation. Snowpack depths were reduced by approximately 700 mm over this period (Figure 2). After this period, snow melt rates more than doubled and remained elevated until complete ablation had occurred.

Meltwater discharge records for the two lysimeters tracked each other well and were similar to those observed previously for a mid-continental area (e.g. Cline, 1997a). Discharge was first recorded in lysimeter 1 on 17 May and continued until 24 June (Figure 3). Maximum daily discharge values were commonly in the range of 2-4 mm hr-1, with a maximum value of 7.9 mm hr-1 recorded on 21 June. The discharge pattern for lysimeter 2 was similar, although between 30 May and 16 June daily maximum discharges were approximately 1mm hr-1 less than those from lysimeter 1. At the start of the melt season, hydrographs were characterized with extended rising limbs and long, gradually declining recession limbs. Little discharge was recorded from 26 to 30 May as new snow accumulated. On 1 June, maximum discharge for lysimeter one was 7.7 mm hr-1 and may have represented the sudden release of water stored in the lysimeter plumbing blocked by freezing. After 1 June, the rising limbs became more abrupt with the recession limbs decaying in an exponential fashion. Except for the very end of the melt season the net energy flux became negative during the night and discharge on the recession limb of the hydrograph approached zero in the morning hours.

A general diurnal pattern of surface melt was evident throughout the record (e.g. Figure 4a). Each night surface melt was reduced to zero for at least several hours, and each day was marked by a distinct melt curve. Once the snowpack began to gain energy, it generally took approximately 4 hours (to between 10:00 and 12:00 hrs) to satisfy the nightly energy deficit. Maximum melt rates peaked at approximately 5.5 mm hr-1. There was considerable day-to-day variation in maximum snow melt, reflecting variations in synoptic-scale weather systems affecting the Front Range, as reported by Cline (1997b) for the 1994 snowmelt season at our study site.

Diurnal melt curves were rarely a simple 'sine-type' curve as has been used in the past to model meltwater input (e.g. Colbeck and Davidson (1973)). Rather, in most cases the curves were multi-modal, with several peaks. For example, on 17 June snow melt began at 6:00 hr and continued to increase in magnitude to 2.6 mm hr-1 at 10:00 hr (Figure 4a). The rate of snow melt than decreased to about 2 mm hr-1 until 13:00 hr, at which time snow melt increased to a rate of 4.1 mm hr-1 at 14:00 hr. Radiation data show a large decline in incident shortwave radiation during the noon hours on this day. Recorded field observations show that convective clouds are a common occurrence around noon at this site. The consequent reduction in shortwave radiation received would account for many of the multi-modal melt curves observed.

The diurnal pattern of snowmelt was, in general, mirrored by that of discharge collected at the bottom of the snowpack, but also showed important differences (e.g. Figure 4a). In the morning hours on 17 June, discharge lagged surface melt by 2-3 hr. At noon, the discharge rate of 3.2mm hr-1 was greater than the preceding maximum melt rate of 2.6 mm hr-1 as a result of kinematic wave flow through the snowpack. Discharge then decreased in response to decreased surface melt at noon and then increased again at 15:00 hr in response to the 14:00 hr increase in surface melt. However, the maximum discharge rate in the afternoon of 3.4 mm hr-1 was less than the maximum surface melt of 4.1 mm hr-1.

Maximum daily discharge values are generally expected to be larger than rates of snowmelt due to the kinematic wave effect. Slow-moving, small surface meltwater fluxes in the morning are overtaken by faster-moving, larger afternoon fluxes, resulting in a kinematic wave moving through the snowpack. Maximum values of Q should be higher than those calculated for U on any given day, since peak Q will consist of the maximum value of U plus earlier meltwater entrained in the wave. This was often not the case in 1997, with peak U often greater than peak Q (Figure 4a). This appears to reflect an interference between two kinematic waves. The first, with Q peaking at approximately noon on 17 June shows the ‘classic’ pattern, whereas the second, peaking at 18:00hrs is lower than expected. This wave has been propagated through a snowpack which already has a high liquid water content due to a lack of drainage of meltwaters from earlier in the day, and is not as pronounced when it reaches the base of the snowpack.

Wavespeeds ([dz/dt]U)

The velocity of liquid water movement was calculated for days when the measured lysimeter discharge was within +30% of the calculated snowmelt volume. This value is arbitrary and was chosen as a conservative approach to ensure that vertical flows through the snowpack above the lysimeters was dominant. For a given U, the time at which it occurred was determined from daily melt curves (Figure 4a). Cumulative U at that time was then determined from the cumulative frequency curves (Figure 4b). The time difference until the equivalent cumulative Q was recorded gave the travel time (dt) for that U value to move through the full snowpack depth (dz). Snowpack depth was divided by this transit time to give a [dz/dt]U for that U value (McGurk and Kattelmann, 1986).

Wavespeeds ([dz/dt]U) were calculated using data from lysimeters 1 and 2 on 12 days in 1997 for a total of 15 sets of values (Table 1). For each day analyzed, [dz/dt]U values were calculated for a range of U on the descending limb of the hydrograph. The [dz/dt]U values [dependent variable) for a given day were then regressed on U (independent variable) using a power relationship on a log-log plot for individual days (Figure 5), generating the relationship:

[dz/dt]U = aUb Eq 4

where a is the y-intercept and b is the slope.

As predicted by equation (3). Slope b was positive in all cases (Table 1), indicating that [dz/dt]U increased with increasing U, and R2 values for these slopes generally exceeded 0.5. Values of b ranged from 0.1348 (lysimeter 1, 9 June) to 0.8209 (lysimeter 2, 14 June), with a mean of 0.4804 (SD = 0.1836) (Table 1). However, this relationship did not change consistently between days, suggesting the flow exponent b (equation 4) was not responding to a progressive change in some snowpack property such as decreasing snow depth or increasing grain size.

Supporting the findings of McGurk and Kattelmann (1986), wavespeeds decreased over time (Figure 6). To compare wavespeeds in the range of measured peak flows over the season, [dz/dt]U was calculated for each day at U = 2 mm hr-1, using the equation for the regression line developed for that day, a value we call [dz/dt]2mm hr-1. During the 30 day period between 21 May and 20 June, [dz/dt]2mm hr-1 declined from 176.6 mm hr-1 to 61.6 mm hr-1, a 65% reduction (Figure 6). A linear regression equation fitted by least-squares to these data has an R2 of 0.54 (p = 0.002), indicating that for these data there was a significant decrease in wavespeeds as the melt season progressed. Over the same time period snow depth decreased from 2,577 mm to 820 mm (Figure 2).

On individual days there appeared to be an inflection point in the [dz/dt]U – U relationship near U = 1 mm hr-1 (Figure 7). R2 for this relationship when only U values equal to and less than 1 mm hr-1 were used was 0.02 with p > 0.05 (n = 34). However, for U values equal to and greater than 1 the R2 improved to 0.36 with p << 0.001 (n = 84). This suggests that the significance and nature of the relationship between [dz/dt]U and U was dependent on the range of values for U, with a steeper regression line for higher U values.

Pore-size distribution index (e) and equivalent permeability (ke).

Fifteen values of e were calculated from the 15 flow exponents b (equation 4), which is equivalent to (e - 1)/e value (equation 3). Values of e ranged from 1.16 to 5.58, with a mean of 2.24 and standard deviation of 1.10 (Table 1). Removing the outlier of 5.58 on 17 June, e values ranged from 1.16 to only 3.32. Perhaps more important, e was found to be 2.93 when all surface meltwater fluxes > 1 mm hr-1 were used to determine a value for b (discussed above as representing an apparent inflection point in the [dz/dt]U – U relationship, Figure 7). This is similar to the value of 3 often used as a default when modeling meltwater flow through snow. A simple linear regression shows that e increased with time especially with the outlier value of 5.58 removed, but the relationship remained non-significant ( a = 0.05).

Fifteen values of ke were then calculated from the 15 values of e . A mean wave wavespeed ([dz/dt]mean U) was calculated for each day that wave speed analysis was conducted, using equation 4 and mean U for that day. The effective porosity (f*-1) was calculated from measured snowpack density, assuming an irreducible water content of 0.07, following McGurk and Kattelmann (1986). These values were then substituted into equation (3), which was then solved for ks. As ks was calculated from wavespeeds over the full snowpack depth, it represented an integrated property of the entire snowpack, and was thus an 'equivalent' permeability (ke). The calculated ke ranged over two orders of magnitude, from 3.23 x 10 -12 to 6.57 x 10 –10 m2, with a mean of 6.14 x 10 –11 m2 (SD = 1.66 x 10 –12 m2) (Table 1). Such a large range in ke was in part caused by the large range in e values calculated previously. Because of the importance of e in equation (3), any variation in e was magnified greatly when ke was calculated using these values. We also calculated ke when e was held constant at 3 in order to allow comparison with ke reported in the literature (Table 1). Holding e constant greatly reduces the variability in ke. With e = 3, the range of ke was reduced by an order of magnitude to 1.08 x 10 –11 m2 to 2.98 x 10 –10 m2, with a mean value of 1.11 x 10 –10 m2 (SD = 8.22 x 10 –11 m2).

Because snowpack properties such as grain size change over the snow melt season, ke may be expected to change as well. A regression of ke versus time using the measured values of e showed no significant change over time (R2 = 0.03; p > 0.05). However, when the pore size distribution index was held constant at e = 3, the calculated ke decrease over time was significant at the a = 0.05 significance level (R2 = 0.40; p = 0.02). With e = 3, ke declined from 3 x 10 –10 m2 to 1 x 10 –11 m2, a 96% reduction in 30 days (Figure 8).


Our values of [dz/dt]2mm hr-1 range from 61 to 176 mm hr-1 and are low compared to those reported in the literature. For example, McGurk and Kattelmann (1986) reported [dz/dt]2mm hr-1 values of about 360 mm hr-1in a Sierra Nevada snowpack in March 1984, declining to approximately 180mm hr-1 by May. They report even higher values in March 1985 (approximately 540 mm hr-1) with a decline by May to values similar to those of the previous year. Jordan (1983a) found [dz/dt]2mm hr-1 values between approximately 250 and 300 mm hr-1in late June/early July, at an alpine site in British Columbia. Colbeck and Davidson (1973) calculated [dz/dt]2mm hr-1 rates of approximately 288 and 576 mm hr-1 from drainage experiments with 0.6 and 1.8 m columns of homogeneous snow and with densities between 600 and 650 kg m-3. Colbeck (1977) and Colbeck and Anderson (1982), using data from a 55.7 m2 snow lysimeter at the Central Sierra Snow Laboratory collected in April 1954 (US Army, Corp of Engineers, 1956), calculated [dz/dt]2mm hr-1 values of approximately 685 and 755 mm hr-1. Similar values were reported in March 1973 from a snow lysimeter 0.7 m2 in area located in the Sleepers River Research Watershed, Vermont (Anderson et al., 1977). All of these studies were conducted in predominately maritime climates near either the west or east coast of North America.

Similarly, despite the large range in ke calculated in this study, there was only limited overlap between them and values calculated at other field sites. Many daily values reported here were over an order of magnitude smaller than the minimum values reported elsewhere. Denoth et al. (1979) report some of the highest ks values in the literature, ranging between 1.0 x 10-9 and 1.0 x 10-8 m2 for disaggregated wet snow with a porosity of approximately 0.48-0.49 and grain size > 1 mm. These values increased to 2.5 x 10–8 m2 for disaggregated, wet snow with a porosity of 0.52 and grain size of approximately 3 mm. As k s values represent a rather different snowpack property than k e values, as already discussed, the limited overlap in ranges, and lower k e values might be expected. The other reported permeability values used here were calculated in a way as to be more similar to k e as we have defined it, so direct comparison and altered notation is more appropriate. Somewhat lower values of ke equal to 1.0 x 10-10 m2 were reported by Colbeck and Davidson (1973) for wet snow with a porosity of 0.32 and density of 650 kg m3, with a slightly higher value of 3.0 x 10-10 m2 for a higher porosity of 0.36 and lower density of 620 kg m3. For a spring snowpack in the Sierra Nevada with a porosity of 0.47 and density of 480 kg m3, more similar to the snowpack properties found in the Colorado Front Range, Colbeck and Anderson (1982) found ke was 2.6 x 10-9 m2, and a higher value of 3.0 x 10-9 m2 for the Vermont snowpack, with a porosity of 0.55 and density of 410 kg m3. Jordan (1983a) does not determine a value of ke from eake1/ef*-1, which he found to approximately vary between 70 and 250, the equivalent of over a range of magnitude in ke. Calculating ke from the data presented by Jordan (1983b) and using the standard values adopted by Jordan (1983b) of e = 3 and f* = 0.43 (from density measurements), for eake1/ef*-1 = 130, ke = 3.29 x 10-10 m2. Over their two year study period, McGurk and Kattelmann (1986) report ke ranging from 3 x 10-10 m2 to 1.8 x 10-8 m2.

There is no obvious single reason for the generally lower values of [dz/dt]U and ke reported here compared with elsewhere. The lower values may be explained by a combination of factors, including physical properties of the snowpack, and characteristics of the continental, alpine climatology of the site. Compared with some investigations, the level of heterogeneity in our alpine snowpacks may have been much higher. For example, Colbeck and Davidson (1973) used repacked, homogenized snow columns in their experiments. Denoth et al. (1979) used glacier firn. Both of those experiments were likely to involve lower grain size variations and fewer ice lenses, resulting in higher wavespeeds and higher ke.

Another important difference is the synoptic-scale climatology at our continental site. During the spring and early summer, due to its high altitude there are rarely long periods of stable, warm conditions resulting in sustained high melt rates over many days. Instead, frequent cool periods, characterized by greatly reduced melt rates, are interspersed throughout the melt season, reducing the levels of saturation in the snowpack which may decrease effective permeability and wavespeeds. Antecedent conditions were not considered when selecting days for analysis, and it is possible that the snowpacks in this study often had lower levels of saturation than those found in other studies. This could result in lower wavespeeds overall, even though the actual surface meltwater fluxes used in this study were similar to those used by McGurk and Kattelmann (1986).

The 50% reduction in [dz/dt]2mm hr-1 values over the 1997 season observed at this continental, alpine site was similar to the decrease reported for two seasons by McGurk and Kattelmann (1986). The decrease in wavespeed (Figure 6) and the marked reduction in ke (Figure 8) over time may be caused by the break down of ice layers in the snowpack. Furbish (1988) suggested that ice layers in a snowpack can reduce transit times. Ice layers may act to increase wavespeeds because saturated flow over the impermeable ice layers has a much greater velocity than unsaturated flow. This leads to a rapid concentration of meltwater in a limited area of the snowpack, and once this water breaks through the ice layer, k w of the snowpack below will be markedly higher because k w increases exponentially with increasing effective saturation (Colbeck, 1978). This will result in preferential flowpaths, with local meltwater travel times much faster than in the surrounding snowpack. Preferential flowpaths are known to occur frequently in the snowpack on Niwot Ridge (Williams et al., 1999a), where their hydrological effects are a subject of on going investigation.

Wankiewicz (1979) attempted to quantify the effect that concentration of meltwaters caused by initial ponding over ice layers may have on flow. He suggested that preferential flowpaths would yield travel times faster by a factor F than those in the uniformly wetted case. He further proposed that F would be related to e as a power function, such that meltwater would move faster by F (e-1)/e for a snowpack with an initial water excess and and F if there is a water deficit. Once the snowpack becomes "ripe", the increase in vertical flow velocities which resulted along these preferential pathways may be sufficient to compensate for the delay resulting from horizontal flow above the ice layers (which would be relatively short if saturated flow occurred). The ice layers would have the overall effect of reducing transit times for a meltwater flux moving from the surface through the entire snowpack depth to the lysimeters at the snowpack base. Consequently, we would have calculated higher wavespeeds and ke values in our study.

As ice layers breakdown during the melt season, these preferential pathways are lost, resulting in more uniform meltwater movement through the snowpack, with decreased relative permeability due to reduced levels of saturation. This would cause longer transit times for the same snowpack depth, resulting in decreased ke. Our measurements of snow stratigraphy lend qualitative support to this argument. On 13 May the snowpack was characterized by four distinct ice layers. Over time these ice layers decayed and disappeared. This corresponds with the observed decrease in effective permeability.

The trend of increasing e over time is also consistent with accelerated movement of meltwater through preferential flow channels early in the melt season followed by a more uniform and slower wave front later in the melt season. Increasing e over time indicates that the distribution of pore-sizes increases over the same time period. This is opposite of what may be expected given that it is generally accepted that wet snow metamorphism leads to increasingly homogeneous, large, and rounded snow grains (e.g. Colbeck, 1986). However, if meltwater flow was predominately through preferential flow channels early in the melt season, wet snow metamorphism in these channels may have resulted in a relatively small range of pore sizes. The change to matrix flow later in the melt season may result in meltwater coming into contact with a greater range of pore sizes. Thus, the range of effective pore sizes that liquid water moving through the snowpack comes into contact with may increase even while the range of total pore sizes in the snowpack is decreasing, as the two would only be the same if the snowpack was saturated and so a two phase system.

An alternative suggestion is that the apparent trend in e values is an artifact of our experimental methodology. Peak melt rates, and thus the range of U values used in analysis, were slightly lower earlier in the melt season. As the relationship between [dz/dt]U - U is dependent on the range of U (Figure 7), and e values were calculated from these relationships, it is possible that this small change in the overall nature of the relationship results in an increase in e over the course of the melt season.

Other studies have reported the dependence of the [dz/dt]U - U relationship on the range of U values (e.g. McGurk and Kattelmann (1986)). However, there has been no attempt at explanation, or discussion of its implication for the use of simple gravity-flow models for describing meltwater movement through snowpacks. In our study there was no significant relationship between [dz/dt]U - U for surface meltwater fluxes below 1 mm hr-1. However, there was a clear relationship between these two variables when higher fluxes were included. For the [dz/dt]U - U relationships at higher ranges of U, the slope of the relationship tends to be higher, and when used to calculate e values, gives results more similar to those reported in other snow studies and values reported for soil textures with similar grain sizes as snow (Muelem, 1978). For example, in our study the regression line of the [dz/dt]U - U relationship for all wavespeeds calculated in 1997 from U values greater than 1 mm hr-1 was y = 68.457 x 0.608, giving e = ~ 3, the value used in most snow hydrology studies. As these e values derived from calculation and laboratory experimentation in other snow and soil hydrology studies generally assume gravity-flow, it appears that only at relatively high surface meltwater fluxes, >1 mm hr-1, can meltwater movement be accurately described by this simple model.

For low U values, capillary affects seem to have an important influence on meltwater movement through the snowpack. At lower surface meltwater fluxes, more pore-spaces are air-filled and the conducting portion of the snowpack's cross-sectional area is reduced relative to times of higher meltwater flux. As snow is such a coarse-textured medium, at low levels of saturation much of the liquid water available for flow will be held in capillary wedges in contact with the grains boundaries, and the large, empty pores are circumvented, causing increased tortuosity. Similarly, the relative importance of these different flow regimes will also be influenced by the concentration of meltwaters in certain areas of the snowpack along preferential flowpaths, and antecedent moisture conditions, reflecting recent surface melt history. These results can be interpreted to suggest that there is a gradual switch between meltwater movement dominated by gravity-flow at higher values of U, to movement dominated by capillary-flow at lower values of U. A 'clean break' in the slope of the regression line describing the [dz/dt]U - U relationship would not be expected, as presumably there is a range of U values over which both capillary-flow and gravity-flow are contributing significant amounts of runoff to the collecting lysimeters, and so the resulting curve would be some combination of the two.

This change between flow regimes may also occur diurnally. There may be a switch from capillary-flow and low lysimeter discharge early in the day, to gravity-flow and high lysimeter discharge when meltwater fluxes exceed 1 mm hr-1. These diurnal changes in flow regime may explain daily patterns of chemical concentration in meltwaters. For example, Bales et al. (1993) and Harrington and Bales (1998) reported that the conductance of snowpack meltwaters increased on the recession limb of the diurnal hydrograph and then decreased sharply on the rising limb. This linkage of chemistry and hydrology is consistent with capillary water at low U rates flowing around grain boundaries and removing impurities, while at high U rates there is gravity-flow which has been moving through snowpack pore-space with limited interaction with individual grains and hence lower conductance.


Surface meltwater fluxes and lysimeter discharge were used to calculate wavespeeds of vertical meltwater movement at a continental, alpine site. These wavespeeds, and equivalent permeabililties for the entire snowpack depth derived from them, were found to be considerably lower then those reported in the literature. Similar to other studies, they exhibited a marked decrease over the melt season, possibly due to the declining influence of ice layers. Although not a significant trend, e values seemed to increase over the melt season. This is contrary to the change expected given simple wet snow metamorphism, and is interpreted to suggest a change from meltwaters moving through predominantly preferential flow channels early in the melt season to a matrix-type flow as the melt season progresses. The nature of the relationship between U and [dz/dt]U , and thus ke and e values, were found to be dependent on the range of [dz/dt]U values used. At high U values the relationship is more similar to that predicted by a gravity-flow dominant regime, but breaks down for lower U values. At these lower values capillary-flow may become more important.

The switch between these two flow regimes will often occur diurnally, with a switch between capillary-flow and low lysimeter discharge during the morning, and gravity-flow and high lysimeter discharge later in the afternoon, this would also go some way to explaining daily patterns of chemical concentration in the meltwaters (Bales et al., 1993). Concentrations are high at low discharge as capillary water flows around grain boundaries removing impurities, and lower at high discharge, when flow consists largely of gravity-flow which has been moving through snowpack pore-space with limited interaction with individual grains. All the wavespeeds analyzed in this study were based on the falling limb of the diurnal hydrograph, when gradual drainage would result in a slow change from gravity- to capillary-flow. The switch in the opposite direction on the rising limb of the daily cycle, should be more pronounced, and occur at the meltwater wave front where there is a rapid increase in the level of saturation in the snowpack. Further work is required to detect these processes

Given that a switch between gravity-flow dominant and capillary-flow dominant regimes seems to occur at relatively high meltwater flux values over 1mm hr-1, and possibly as high as 2 mm hr-1, this has considerable implications for snow hydrology models which ignore capillary effects in their treatment of meltwater movement through snow. Some degree of error is introduced at lower meltwater fluxes, in which much of the seasonal melt will occur under natural conditions at many locations. These results suggest that current models of meltwater flow through snow need to incorporate capillary-flow processes when melt rates are less than 1 mm hr-1. However, modeling meltwater flow as gravitational-flow at higher melt rates appears to be correct.


Q lysimeter discharge (mm hr-1)

DQS cold content of the snowpack

DQM latent heat flux during melting or freezing

Q* net all wave radiation flux

QH sensible heat flux

QE latent heat flux

QG ground heat flux

QR heat advected by precipitation

S* effective water saturation of snow, expressed as a fraction of pore volume occupied by mobile liquid

U surface meltwater flux (mm hr-1)

[dz/dt]U velocity at which flux U moves through the snow (m s-1) or (mm hr-1)

a flow factor constant = 5.47 x 106 m-1 s-1 at 0oC

kw relative, or liquid water permeability (m2)

ke snowpack equivalent permeability (m2)

ks saturated permeability (m2)

f* effective porosity





Ambach, W. M., M. Blumthaler, and P. Kirchlechner, 1981 Applications of the gravity flow theory to the percolation of melt water through firn, Journal of Glaciology, 27, 67-75.

Anderson, E. A., R. Z. Whipkey, H. J. Greenan, and C. T. Machell, 1977. NOAA –ARS Cooperative Snow Research Project-Watershed hydro-climatology and data for years 1960-1974, report US Dept of Commerce, Silver Spring, Md.

Bales, R. C., R.E. Davis and M.W. Williams, 1993. Tracer release in melting snow: diurnal and seasonal patterns. Hydrological Processes, 7, 389-401.

Caine, N., 1992. Modulation of the diurnal streamflow response by the seasonal snowcover of an alpine basin, Journal of Hydrology, 137, 245-260.

Caine, N., 1995. Snowpack influences on geomorphic processes in Green Lakes Valley, Colorado Front Range, Geographical Journal, 161, 55-68.

Caine, N., 1996. Streamflow patterns in the alpine environment of North Boulder Creek, Colorado Front Range, Zeitschrift fur Geomorphologie, 104, 27-42.

Cline, D. W., 1997a. Snow surface energy exchanges and snowmelt at a continental, midlatitude Alpine site. Water Resources Research, 33(4), 689-701.

Cline, D. W., 1997b. Effect of seasonality of snow accumulation and melt on snow surface energy exchanges at a continental alpine site, Journal of Applied Meteorology, 36, 32-51.

Colbeck, S. C., 1972. A theory of water percolation in snow. Journal of Glaciology, 11, 369-385.

Colbeck, S. C., 1973. Effects of stratigraphic layers on water flow through snow. Cold Regions Research and Engineering Laboratory, Research Report 311.

Colbeck, S. C., 1974. The capillary effects on water percolation in homogeneous snow. Journal of Glaciology, 13, 85-97.

Colbeck, S. C., 1975. A theory of water flow through a layered snowpack. Water Resources Research, 11, 261-266.

Colbeck, S. C., 1977. Short-term forecasting of water runoff from snow and ice. Journal of Glaciology, 19, 571-588.

Colbeck, S. C., 1978. The physical aspects of water flow through snow in Chow, V. T. (ed.), Advances in Hydroscience, VII, 165-206.

Colbeck, S. C., 1979. Water flow through heterogeneous snow, Cold Regions Science and Technology, 1, 37-45.

Colbeck, S. C., 1986. Statistics of coarsening in water-saturated snow, Acta Metallurgica, 34, 347-352.

Colbeck S. C. and G. Davidson, 1973. Water percolation through homogeneous snow. In The role of snow and ice in hydrology: Proceedings of the Banff Symposium, IAHS, 242-256.

Colbeck, S. C. and E. A. Anderson, 1982. The permeability of a melting snow cover. Water Resources Research, 18, 904-908.

Denoth, A., W. Seidenbusch, M. Blumthaler, P. Kirchlechner, W. Ambach and S. C. Colbeck, 1979. Study of water drainage from columns of snow. Cold Regions Research and Engineering Laboratory, Research Report 79-1.

Furbish, D. J., 1988. The influence of ice layers on the travel time of meltwater flow through a snowpack, Arctic and Alpine Research, 20, 265-272.

Harrington, R. K. and R. C. Bales, 1998. Interannual, seasonal, and spatial patterns of meltwater and solute fluxes in a seasonal snowpack. Water Resources Research, 34 (4), 823-831.

Jordan, P., 1983a. Meltwater movement in deep snowpack 1. Field observations. Water Resources Research, 19 (4), 971-978.

Jordan, P., 1983b. Meltwater movement in deep snowpack 1. Simulation Model. Water Resources Research, 19 (4), 979-985.

Marks D., J. Dozier and R. E. Davis, 1992. Climate and energy exchange at the snow surface in the alpine region of the Sierra Nevada. 1. Meteorological measurements and monitoring. Water Resources Research 28 (11) 3029-3042.

Marsh, P., 1991. Water flux in melting snow covers. In: Advances in Porous Media (ed. by M. Y. Corapcioglu), vol. 1, 61-124. Elsevier Science Pub. Co., New York.

Marsh, P. and Woo, M. K., 1985. Water movement in natural heterogeneous snow covers. Water Resources Research, 21, 1710-1716.

McGurk, B. J. and Kattelmann, R. C., 1986. Waterflow rates, porosity and permeability in snowpacks in the Central Sierra Nevada. In Kane, D. L. (ed.), Cold Regions Hydrology Symposium, American Water Resources Association Tech. Publ. Ser. TPS-86-1, 359-366, American Water Resources Association, Bethesda, MD.

McGurk, B. J. and Kattelmann, R. C., 1988. Evidence of liquid water through snow from thick-section photography, in Proceedings of the International Snow Science Workshop, 137-139, Canadian Avalanche Association, Vancouver.

Richter-Menge J. A. and S. C. Colbeck, 1991. Recent progress in snow and ice research. Reviews of Geophysics, Supplement, April 1991, 218-226.


Schneebeli, M., 1998. Unsaturated water flow in snow: experiment and simulation, US Army Cold Regions Research and Engineering Laboratory Special Report 98-10.

U.S. Army. Corp of Engineers, 1956. Snow Hydrology: Report of the Cooperative Snow Investigations. Corps of Engineers, Portland, Oregon.

Wankiewicz, A., 1978. Water pressure in ripe snowpacks. Water Resources Research, 14, 593-600.

Wankiewicz, A., 1979. A review of water movement in snow. In Proceedings Modeling of Snow Cover Runoff (ed. Colbeck and Ray), 222-253. U.S. Army CRREL, Hanover, New Hampshire.

Williams, M. W. and J. M. Melack, 1989. Effects of spatial and temporal in snowmelt on nitrate and sulphate pulses in meltwaters within an alpine basin, Annals of Glaciology, 13, 285-289.

Williams, M. W., M Losleben, N. Caine, and D. Greenland, 1996a. Changes in climate and hydrochemical responses in a high-elevation catchment, Rocky Mountains, Limnology and Oceanography 41 (5), 939-946.

Williams, M. W., P. D. Brooks, A. Mosier, and K. A. Tonnessen, 1996b. Mineral nitrogen transformation in and under seasonal snow in a high-elevation catchment, Rocky Mountains, USA, Water Resources Research, 32, 3175-3185.

Williams, M. D., R. Sommerfeld, S. Massman, and M. Rikkers, 1999a. Correlation lengths of vertical flowpaths in melting snowpacks, Colorado Front Range, USA. Hydrological Processes, 13, 1807-1826.

Williams, M. W., D. Cline, M. Hartman, and T. Bardsley, 1999b. Data for snowmelt model development, calibration, and verification at an alpine site, Colorado Front Range, Water Resources Research, 35 (10) 3205-3209.


Table 1. Unsaturated hydrological properties of a continental, alpine snowpack, 1997.








[dz/dt]2 mm hr-1 2

(mm hr-1)
































































































































1constant (a) and exponent (b) in the regression equation, y = axb

2normalized wavespeed (calculated for U = 2 mm hr-1 from the regression equation)

xke is the equivalent permeability, calculated using the e values calculated from equation 3

bke is the equivalent permeability, calculated using the e equals 3





Figure 1. Location map of Niwot Ridge.

Figure 2. Snowpack depth and depth averaged temperature.

Figure 3. Hourly lysimeter discharge averaged from 10 minute values for the entire 1997 melt season.

Figure 4. (a) Surface meltwater flux (U) and lysimeter 2 discharge (Q) on 17 May, 1997. (b) Cumulative percentage curves of surface meltwater flux (U) and lysimeter 2 discharge (Q) on 17 May, 1997.

Figure 5. A log-log plot of wavespeed ([dz/dt]U) regressed on 8 surface meltwater flux (U) values, calculated from lysimeter 2 on 17 May, 1997.

Figure 6. The change in [dz/dt]2 mm hr-1 and best fit line from 21 May to 20 June, 1997. A simple linear regression shows a significant decline from 176 mm hr-1 to 61 mm hr-1 over the 30 day period (R2 = 0.54, p = 0.002).

Figure 7. A log-log plot of wavespeed ([dz/dt]U, mm hr-1) regresssed against surface meltwater flux (U). All available data in 1997 (n = 118) was split into U < 1 mm hr-1 (n = 34) and U > 1 mm hr-1 (n = 83)

Figure 8. The change in ke and best fit line from 21 may to 20 June, 1997, using e = 3 in equation 3. A simple linear regression shows a significant decrease at the 95% significance level (R2 = 0.40, p = 0.02).