List of symbols and units

Introduction

The one-dimensional, gravity theory of water percolation through snow was described in a previous paper (Reference ColbeckColbeck, 1972). The role of capillarity was omitted and only the driving force of gravity was considered. This simplification was justified by showing that the percolation processes were dominated by gravity everywhere except at the leading edges of melt waves and at very low flow rates. The theoretically constructed waves compared favorably with the measured waves (Reference ColbeckColbeck, 1972; Reference Colbeck and DavidsonColbeck and Davidson, in press) and it was concluded that the gravity-flow theory was capable of making accurate predictions of water percolation through snow as long as the snow was fairly homogeneous and the flow rates were not too small. For the purpose of increasing our understanding of more complicated situations, the basic theory of capillary effects in snow is developed in this paper. First the results of capillary-pressure experiments performed in the laboratory with kerosene and snow are presented. The governing equations are next derived and a numerical method presented for their solution. An analytical solution is also presented which is applicable to the region of the wave front and provides new insight from consideration of capillary effects at various positions along a typical melt-water wave. Finally, end effects are examined.

It is difficult to apply two-phase, Darcian concepts to snow because of the lack of experimental data on the permeability and capillarity of water-bearing snow. Indirect methods were used (Reference Colbeck and DavidsonColbeck and Davidson, in press) to find the relationship between relative permeability and water saturation for snow and the same is done here for the relationship between capillary pressure and water saturation. While it must be recognized that any quantitative conclusion drawn from this indirect evidence is tentative, the basic theory of capillary flow can be established.

The capillary-pressure curve

Water-bearing snow consists of rounded grains of uniform size with a small amount of liquid water. Snow has a uniform pore structure and is expected to have a relationship between capillary pressure and water saturation characteristic of materials with a similar structure. As an approximation to the snow-water system, capillary-pressure experiments were performed with snow and kerosene at a temperature of —10° C in a cold laboratory. The snow, which consisted of small grains (dia. ≈ 1 mm), was packed to a density of 0.56 Mg m⁻³ in the sample holder, saturated with kerosene and then drained through a semi-permeable barrier under careful capillary control (see, e.g. Reference ScheideggerScheidegger, 1957, p. 48). The result (sec Fig. 1) is fairly typical of materials with a uniform pore size. When the capillary pressure was increased from 320 to 540 N m⁻², 45% of the liquid drained. This rapid displacement occurs because capillary pressure is dependent on the effective pore radius and only a narrow range of pore sizes occurs in snow. The other important feature of this curve is the value of the irreducible liquid saturation, 0.07. This value is very close to that found by Reference Harris and MorrowHarris and Morrow (1964) in random packs of quidimensional spheres and is probably a good estimate of the value of S_wi for most types of snow.

Fig. 1. Capillary pressure as a function of liquid saturation. The experimental data points are shown.

Only a small part of the curve is significant with regard to water drainage in homogeneous snow because water saturations generally exist within a narrow range (about 0.10 to 0.20). In Figure 2 capillary pressure is plotted as a function of effective liquid saturation S* and in Figure 3 the equation

is shown to be a reasonable representation of the data. Equations (1) is integrated to give

which is shown on Figure 2 with the experimental curve. Equations (2) is an accurate representation of the experimental data over the range of saturations expected during water percolation in homogeneous snow.

The relationship between capillary pressure and water saturation depends on a number of parameters including contact angle, interfacial tension, particle shape, and pore-size distribution. While the experimental results must be cautiously interpreted, they agree qualitatively with the results which would be expected for water-bearing snow. Only the value of the coefficient in Equations (1) is likely to be changed when capillary pressure data becomes available for natural snow. The algebraic form of Equations (1) is convenient and when the governing equations are derived, analytical solutions are possible.

Fig. 2. Capillary pressure as a function of effective liquid saturation. The solid line was calculated from Figure 1 and the dashed line represents Equations (2)

Theory

Using two-phase, Darcian flow concepts it has been shown that (Reference ColbeckColbeck, 1972)

In the development of this equation it was argued that:

since water saturations are generally less than 0.20 (Reference Colbeck and DavidsonColbeck and Davidson, in press), and, (c) that the total volume flow (air plus water) must be zero.

Normally there is a hysteresis in the relationship between capillary pressure and water saturation because of the “trapping” which characteristically occurs during imbibition (Reference ScheideggerScheidegger, 1957, p. 49). Since the range of saturations is small in simple drainage, hysteresis is neglected and capillary pressure P_c is taken as a single-valued function of effective water saturation S*. Equations (3) becomes

The generally assumed form of the relationship between relative permeability and water saturation is

and the continuity equation is

where Φ_e is the “effective porosity” which includes only the pore space which is available for flow. Equations (6), (7) and (8) are combined to give the basic equation describing water percolation in homogeneous snow:

where n is about 3 (Reference Colbeck and DavidsonColbeck and Davidson, in press). The diffusive effects of capillarity have introduced second-order terms into Equations (9) which, if capillarity is negligible (∂P _c/∂S* = 0), simplifies to the first-order equation which describes the percolation under the influence of gravity alone.

Fig. 3. ∂P _c/∂S* as a function of S*.

Equations (9) is equivalent to the Fokker-Pianck equation in diffusion theory which has received wide use in the theory of infiltration into soils (e.g. Reference PhilipPhilip, 1969). For three-dimensional flow (taking n = 3),

Where

This equation is difficult to solve although solution for either the gravity or the capillary case is relatively easy. Also, the application of Equations (10) to snow is difficult because of the complicated boundary and initial conditions which occur. In the case where

And

accurate solutions have been obtained for Equations (10) (Reference PhilipPhilip, 1969). While these solutions could be applied to snow, these simple criteria do not approach the complicated nature of the diurnal melt waves which characterize percolation through snow. Equations (10) could be linearized by holding k _w

and constant. While this simplification destroys the dynamic interplay between capillary and gravity forces, it may be a useful approach for describing the three-dimensional character of the flow around deeply buried ice layers or the firn-ice transition. Also these criteria may be useful in describing the drainage from glacial firn during the winter months when little flow occurs across the upper surface.

Numerical Solution

In the case of gravity flow through homogeneous snow, the first-order equation was solved by Reference Colbeck and DavidsonColbeck and Davidson (in press) using the method of characteristics as described by Reference Sheldon, Sheldon and CardwellSheldon and others (1959). The second-order equation developed here can be solved using a variation of that method (see Reference Lighthill and WhithamLighthill and Whitham, 1955). Whereas saturation is invariant along the characteristics in the case of gravity flow, in capillary flow the value of saturation changes along each characteristic. Therefore the characteristics can change shape in the region of a wave front and no discontinuity occurs.

When effective saturation S* varies with space z and time t,

where D/Dz means differentiation along a characteristic. The speed of propagation of any value of S* is

then

Following Reference Lighthill and WhithamLighthill and Whitham (1955), c is assumed to vary only with S* along a characteristic, thus where c is the slope of the curve relating flow u _w to concentration S*,

Upon combining Equations (9), (16) and (17) to eliminate

,

This equation describes the change of S* along a characteristic and when sufficient information has been specified, the solution can be constructed using numerical techniques.

As a characteristic approaches the region of a wave front, ∂S*/∂Z and ∂² S*/∂Z ² become sufficiently large to cause a significant change in DS*/D _z and therefore a change in shape of the characteristic lines.

Combining Equations (1) and (18), the spatial rate of change of S* along a characteristic is given by (taking n = 3)

Under normal conditions of flow, DS*/D_z is less than 2 X 10^–5 m⁻¹ (ΔS* is less than 10⁻⁴) between wave fronts. The slope of the characteristics is given by

where the first term is due to gravity and the second to capillarity. As shown later, under normal conditions of flow the second term is very much smaller than the first term except in the region of a wave front.

The wave front

The complete solution can be constructed from Equations (19) and (20) but analytical methods are better suited for estimating the effects of capillarity on the wave front. The actual shape of the wave front is calculated analytically over the entire region of the front except within the neighborhoods of the two points where the wave front joins the trailing edges. Thus the diffusive effect of capillarity is estimated.

Following the method described by Reference Morel-Seytoux and De WiestMorel-Seytoux (1969, p. 507), a new system of coordinates is defined which advances with the shock front. The shock front is located at (t) relative to the fixed coordinates z and any position Z in the new coordinates is given by

Thus

It follows that

and

In this coordinate system, Equations (10) becomes,

The advantage of working in this coordinate system is that, within the interior of the wave front, S* changes very rapidly with depth but very slowly with time. This is not true, however, in the matching regions between the leading and trailing edges. In these regions the solution breaks down since

is finite while is zero. Therefore the matching condition must be estimated. Within the interior region of the wave front,

At depths below 0.5 m, the construction of the characteristics shows that the shock-front speed (d£/dt) is nearly constant. In fact,

where S is the daily average of S* at the surface. Holding d/dt constant permits direct integration of Equations (26) whence

When

is known from experimental results, Equations (28) can be integrated to give the distribution of saturation within the interior of the wave front.

By substituting Equations (1) into Equations (28) when n is 3,

Integrating between the shock front and any position Z,

where

and S is the value of S* at the shock front. To illustrate Equations (30), a particular case is examined. Taking Φ_e = 0.45, k = 3 X 10^–10 m² and a truncated sine wave for the boundary condition where u_max = 1.20 X 10⁻⁶ m _s–, the gravity-flow solution was constructed using the method of characteristics (see Reference ColbeckColbeck, 1972, for examples). The values of saturation at 4.85 X 10₄ s after the start of melting are shown as a function of depth on Figure 4. Also from this construction it was determined that d/dt ≈ 0.465 X 10_-4 m s_-1 for all depths greater than 0.5 m. Using these particular values in Equations (30),

where the value of S must be found by material balance between the shock front and capillary solutions. The calculated wave front is shown in Figure 5.

Fig. 4. The saturation distribution 48 500 s after the start of melting an day 3. This profile was calculated using the gravity-flow theory.

In this example the tendency for the leading edge of a wave to form into a shock front is very strong because the top of the wave travels about four times as fast as the bottom. While capillarity does cause some departure from a shock front, the diffusive effect which capillarity has on the leading edge is not very pronounced. Wave fronts observed in homogeneous snow (e.g. Reference Colbeck and DavidsonColbeck and Davidson, in press) always have a much more “diffused” appearance than the one shown in Figure 5, This may be explained by the phenomenon known as “fingering”, i.e. diffusion of the wave front resulting from inhomogeneous flow at the local scale. Application of a stability criterion for two-phase flow in homogeneous porous media (Reference Bear, Bear, Zaslavsky and IrmayBear and others, 1968, p. 283) shows that perturbations should not grow. Fingering could result from other causes, however, such as the inhomogeneous nature of the snow itself or a break down of the assumption that air pressure is hydrostatic.

Fig. 5. A melt-water wave with a calculated wave front and values of R along the wane.

While it was possible to deduce the form of k_v(S*) from data on the percolation of water through snow, it does not seem likely that P_e(S*) can be found in a similar manner. This is unfortunate in view of the problems associated with handling water-bearing snow samples during experimental studies. The information available so far indicates that capillarity is of very little importance when the flow is essentially one-dimensional. Further insight into the relative importance of capillarity can be gained by referring back to Darcy's law. From Equations (6) and (7) where n = 3 for snow (Reference Colbeck and DavidsonColbeck and Davidson, in press)

While the gravity force is constant, the capillary force varies with position and time. The ratio of capillary to gravity forces is equal to the ratio of the two terms in Equations (20). Taking

the ratios calculated for various positions along the melt-water wave are shown in Figure 5:

1. At the leading edge,

   

   (35)
   Thus the gravity and capillary forces are nearly equal over the interior region of the wave front except at the matching points where gravity dominates.

2. At intermediate points along the trailing edge,

   

   (36)
   thus capillary forces account for less than 1% of the total force.

3. It is now determined at what flow rate (along the trailing edge) capillarity becomes large. R increases as S* approaches zero although

   
   becomes quite small at the same time. In general,
   

   (37)
   At a flow rate of 10⁻⁷ m s^–1, R is about 0.02; at flow rates of 10 ^–8 m s⁻¹ R is about 0.10; and at 10^–9 m s⁻¹, R is about 0.5. The lowest flow rates observed in homogeneous snow undergoing water flux due to normal surface melting were greater than 10⁻⁷ m s⁻¹ hence capillarity is not important at the lower limit of normal fluxes. However, when surface fluxes stop, the role of capillarity increases as drainage proceeds, until ultimately R approaches infinity as S* approaches zero.

4. The surface effect of capillarity can now be examined. When S* decreases to 0.01 at the surface (see Fig. 5), R is about 23 indicating that capillary effects control surface fluxes during periods when the surface saturations are low. This “surface effect” is probably not as serious as the “end effect” at the outflow end since only very small volumes of mobile water are involved and the surface effect is limited to a region within 0.25 m of the surface. This surface effect causes a skewing of the surface flux similar to that caused by the percolation process itself. However, the slowly moving melt-water wave associated with this surface effect would be overtaken by the next day's wave front before reaching a depth of even 1 m (see Fig. 5).

End Effects

End effects are the increased saturation (above the saturation associated with flow) due to capillary action at a discontinuity of the porous medium. These effects can be seen in water-bearing snow when an interface exists between the snow and an underlying air space (such as an undercut snow pit). Under such circumstances the snow appears to be saturated for about 20 mm above the interface with the saturation rapidly decreasing beyond that distance. If the amount of water trapped by the end effect varies in time, a significant distortion of the travelling melt-water wave could result. End effects occur in snow lysimeters and, since we are dependent upon the results of lysimeter studies to measure the movement of liquid water in snow masses, any difference between the wave movement through the lysimeter and through the in situ snow mass must be detected. Also, the presence of the end effect can confuse tracer studies since there is a time delay between the arrival of a water particle at the upper end of the end-effect region and its discharge at the interface. Note, however, that as long as the saturation profile associated with the end effect remains fixed, the end effect has no influence on travelling waves since a particle of water will be released at the interface as soon as another particle of water arrives at the top of the region of the end effect.

Although Equations (1) is an accurate representation of the experimental data for values of S* below 0.6, it is a poor approximation for the larger values which must be included here. Therefore to gain an estimate of the influence of the end effect on lysimeter response, take

It is shown later that this approximation is sufficient.

Given Equations (38), Equations (33) can be integrated if u _w is held constant. In this manner the saturation profiles are constructed first for various values of flux and then the variation in the profile with time is established. This method relies on the assumption that transient effects are small but, at least for the case of the trailing edge where changes occur very slowly (transit time ≈ 22 h), the assumption is probably valid. Rearranging and integrating Equations (33),

where the snow is assumed to be saturated at the interface. Combining Equations (38) and (39)

Taking as an example φ_e = 0.45 and k = 3 X 10⁻¹⁰ m², Equations (40) becomes

where

The background saturation associated with each value of flux is derived from the gravity-flow theory,

The two extreme values of water flux associated with fair-weather melting (2 X 10⁻⁶ and 10⁻⁷ m s⁻¹) are chosen because the largest change in the saturation profile should occur between these two values. The calculated profiles are shown in Figure 6. The sudden change in the slope at S* = 0.6 occurs because of the approximation made for

(Equations (38)). Because these two profiles are essentially identical for values of S* above 0.30, any approximation made for above 0.6 would have worked. Only the difference between these profiles is significant.

It is now possible to deduce the changing end effect while the trailing edge is passing through a lysimeter. During this time the profile must change from the upper limit to the lower limit. For a lysimeter with radius = 0.113 m, the change in volume of water associated with the passage of the trailing edge is 130 X 10_-6 m³. Thus, while the water flux decreases from 2 X 10⁻⁶ to 10⁻⁷ m s⁻¹, 130 X 10⁻⁶ m³ of water are added to the end effect. This however, is not water which would have to be diverted from the flow and placed in storage because this water is already in place. In fact 9.5 x 10⁻⁶ m³ of water will be released during the time of passage of the trailing edge. With a total volume flux of about 2 000 X 10⁻⁶ m₃ during this period, the amount of water released from the end effect is insignificant.

The changes in end effect associated with the passage of the wave front must balance those just calculated for the trailing edge, i.e. about 9.5 X 10_-6 m_-3 of the water arriving with the leading edge will be stored as end effect. In spite of this, the leading edge will arrive earlier than it would in the absence of an end effect because of the water which is already in place. Thus arrival should occur about

h (a 400 s) faster than otherwise expected. Also, because saturation gradually increases in the direction of wave propagation, the wave front will be diffused by the end-effect phenomenon.

Fig. 6. The end-effect regions associated with the limiting values of u_w expected during fair-weather melting. The end effect includes only the water trapped at the interface by capillary action.

Conclusions

The theory of water percolation through snow is expanded to include the effects of capillarity and, when combined with experimental data, further insight into the nature of water percolation is obtained. Apparently, for all flow rates occurring during periods of strong diurnal melting, the flow processes are dominated by gravity and unless flow rates as low as 10_-8 ms_-1 are reached, capillarity can be ignored. While capillarity is important at the wave front, the shock front approximation is sufficient for most purposes. The possible diffusion of waves by the lysimeter is a serious problem. This should be tested by using lysimeters which neutralize the end effect by reducing the water pressure in the discharge line.